R Émi

Chunmei Su

NUMERICAL STUDY OF THE LOGARITHMIC SCHR ÖDINGER EQUATION WITH REPULSIVE HARMONIC POTENTIAL

Keywords: 2010 Mathematics Subject Classification. 35B05, 35B40, 35Q55, 65M15, 81Q05 Nonlinear Schrödinger equation, logarithmic nonlinearity, repulsive harmonic potential, dispersion, splitting methods, error estimates. RC

We consider the nonlinear Schrödinger equation with a logarithmic nonlinearity and a repulsive harmonic potential. Depending on the parameters of the equation, the solution may or may not be dispersive. When dispersion occurs, it does with an exponential rate in time. To control this, we change the unknown function through a generalized lens transform. This approach neutralizes the possible boundary effects, and could be used in the case of the nonlinear Schrödinger equation without potential. We then employ standard splitting methods on the new equation via a nonuniform grid, after the logarithmic nonlinearity has been regularized. We also discuss the case of a power nonlinearity and give some results concerning the error estimates of the first-order Lie-Trotter splitting method for both cases of nonlinearities. Finally extensive numerical experiments are reported to investigate the dynamics of the equations.

Introduction

We consider the initial value problem of a time-dependent nonlinear Schrödinger equation with potential V and logarithmic nonlinearity:

(1.1)

   i∂ t u + 1 2 ∆u = V u + λu ln |u| 2 , x ∈ R d , t > 0, u(0, x) = u 0 (x),
where u 0 is the initial data.

In the absence of the logarithmic nonlinearity (λ = 0), there is a large amount of literature on the so-called repulsive potential setting [START_REF] Barceló | Some dispersive estimates for Schrödinger equations with repulsive potentials[END_REF][START_REF] Reed | Methods of modern mathematical physics. IV Analysis of Operators[END_REF], i.e., (∂ r V ) + (r) := sup

x∈R d :|x|=r (∂ r V ) + (x) = 0.
Particularly, solutions of the Schrödinger equation with repulsive 1/r n -type longrange potentials play an essential role in many physical processes such as collisions of similar atoms in a radiation field [START_REF] Burke | Multichannel cold collisions: Simple dependences on energy and magnetic field[END_REF][START_REF] Gallagher | Exoergic collisions of cold Na * -Na[END_REF][START_REF] Orzel | Time-resolved studies of ultracold ionizing collisions[END_REF], molecular spectra converging to thresholds where two fragments can interact via resonant dipole-dipole interactions [START_REF] Abraham | Photoassociative spectroscopy of long-range states of ultracold 6 Li 2 and 7 Li 2[END_REF][START_REF] Miller | Photoassociation spectrum of ultracold Rb atoms[END_REF][START_REF] Napolitano | Line shapes of high resolution photoassociation spectra of optically cooled atoms[END_REF], atom-electron, atom-ion interactions [START_REF] Gao | Repulsive 1/r 3 interaction[END_REF] and atom-surface interactions [START_REF] Stwalley | Simple long-range model and scaling relations for the binding of isotopic hydrogen atoms to isotopic helium surfaces[END_REF]. Compared to the free-potential equation, the repulsive potential creates acceleration of the field [START_REF] Antoine | Absorbing boundary conditions for the one-dimensional Schrödinger equation with an exterior repulsive potential[END_REF][START_REF] Lorin | A numerical Maxwell-Schrödinger model for intense laser-matter interaction and propagation[END_REF]. In the case of a repulsive harmonic potential, V (x) = -|x| 2 , the acceleration is exponential in time, as recalled below.

In the absence of potential (V = 0), the logarithmic Schrödinger equation (logNLS) has been adopted in many physical models [START_REF] Avdeenkov | Quantum Bose liquids with logarithmic nonlinearity: Self-sustainability and emergence of spatial extent[END_REF][START_REF] Buljan | Incoherent white light solitons in logarithmically saturable noninstantaneous nonlinear media[END_REF][START_REF] Hansson | Propagation of partially coherent solitons in saturable logarithmic media: A comparative analysis[END_REF][START_REF] Hefter | Application of the nonlinear Schrödinger equation with a logarithmic inhomogeneous term to nuclear physics[END_REF][START_REF] Krolikowski | Unified model for partially coherent solitons in logarithmically nonlinear media[END_REF][START_REF] Martino | Logarithmic Schrödinger-like equation as a model for magma transport[END_REF][START_REF] Zloshchastiev | Logarithmic nonlinearity in theories of quantum gravity: Origin of time and observational consequences[END_REF] since it was introduced in [START_REF] Bia Lynicki-Birula | Nonlinear wave mechanics[END_REF]. For instance, as proposed in [START_REF] Bouharia | Stability of logarithmic Bose-Einstein condensate in harmonic trap[END_REF][START_REF] Zloshchastiev | Spontaneous symmetry breaking and mass generation as built-in phenomena in logarithmic nonlinear quantum theory[END_REF], the logarithmic model may generalize the Gross-Pitaevskii equation, used in the case of two-body interaction, to the case of multi-body interaction. A particular feature of the logarithmic nonlinearity is that it leads to very special solitary waves, called Gaussons in [START_REF] Bia Lynicki-Birula | Nonlinear wave mechanics[END_REF][START_REF] Bia Lynicki-Birula | Gaussons: Solitons of the logarithmic Schrödinger equation, Special issue on solitons in physics[END_REF] when λ < 0. These solitary waves are orbitally stable [START_REF] Ardila | Orbital stability of Gausson solutions to logarithmic Schrödinger equations[END_REF][START_REF] Cazenave | Stable solutions of the logarithmic Schrödinger equation[END_REF]. Furthermore, for λ < 0, no solution is dispersive ( [26, Proposition 4.3]), while for λ > 0, every solution is dispersive with an enhanced rate compared to the usual rate of the free Schrödinger equation, and the modulus of the solution converges to a universal Gaussian profile [START_REF] Carles | Universal dynamics for the defocusing logarithmic Schrödinger equation[END_REF].

When it comes to the logarithmic model with potential, a harmonic trapping potential was considered in [START_REF] Bouharia | Stability of logarithmic Bose-Einstein condensate in harmonic trap[END_REF] to describe the logarithmic Bose-Einstein Condensation:

(1.2)

i∂ t u + 1 2 ∆u = ω 2 2 |x| 2 u + λu ln |u| 2 , x ∈ R d .
Due to the presence of the potential, stationary solutions (generalized Gaussons) are available and orbitally stable in both cases λ < 0 [START_REF] Ardila | Logarithmic Bose-Einstein condensates with harmonic potential[END_REF][START_REF] Bouharia | Stability of logarithmic Bose-Einstein condensate in harmonic trap[END_REF] and λ > 0 [START_REF] Carles | Orbital stability of Gausson solutions to logarithmic Schrödinger equations[END_REF]. For the logNLS with repulsive potential,

(1.3) i∂ t u + 1 2 ∆u = - ω 2 2 |x| 2 u + λu ln |u| 2 , x ∈ R d ,
it was shown in [START_REF] Carles | Orbital stability of Gausson solutions to logarithmic Schrödinger equations[END_REF]Proposition 1.3] that for λ ∈ R and any

u 0 ∈ Σ := H 1 ∩ F(H 1 ) = f ∈ H 1 (R d ), x → |x|f (x) ∈ L 2 (R d ) , f Σ := f L 2 + ∇f L 2 + xf L 2 ,
there exists a unique solution u ∈ L ∞ loc (R; Σ) ∩ C(R; L 2 (R d )) to (1.3). Furthermore, when λ > 0, the solution shares the same exponential dispersive rate as in the linear case λ = 0, and no solitary wave or universal dynamics exists in this case. In the nondispersive case (λ < 0), the situation is different. It was shown in [START_REF] Zhang | Bound states for logarithmic Schrödinger equations with potentials unbounded below[END_REF] that (1.3) admits at least one positive bound state under a suitable range of the coefficients. Recently, we proved that if -λ > ω > 0, there exist two positive stationary Gaussian solutions, which are orbitally unstable [START_REF] Carles | Nonuniqueness and nonlinear instability of Gaussons under repulsive harmonic potential[END_REF]. We emphasize the fact that so far, (1.3) does not seem to be related to a physical model. However, it raises difficulties for the mathematical and numerical analysis, which may help gain insight for other problems.

Along the numerical part, there have been few studies for the model with repulsive potential or logarithmic nonlinearity. For the Schrödinger equation with repulsive potential (1.3) and without logarithmic interaction (λ = 0), the strong dispersive effects made it almost impossible to simulate the dynamics on a truncated domain with naive homogeneous or periodic boundary conditions [START_REF] Antoine | Absorbing boundary conditions for the one-dimensional Schrödinger equation with an exterior repulsive potential[END_REF]. To overcome this difficulty, some types of artificial or absorbing boundary conditions were proposed for one-dimensional Schrödinger equation with a general variable repulsive potential [START_REF] Antoine | Absorbing boundary conditions for the one-dimensional Schrödinger equation with an exterior repulsive potential[END_REF] and for two-dimensional Schrödinger equation with a time and space varying exterior potential [START_REF] Antoine | Absorbing boundary conditions for the two-dimensional Schrödinger equation with an exterior potential part I: construction and a priori estimates[END_REF][START_REF] Antoine | Absorbing boundary conditions for the two-dimensional Schrödinger equation with an exterior potential[END_REF]. On the other hand, for the logNLS without potential, the singularity of the logarithmic nonlinearity also makes it very challenging to design and analyze numerical schemes. To overcome the singularity at the origin, some numerical methods were proposed and analyzed for the logarithmic Schrödinger equation ((1.3) with ω = 0) based on a global nonlinearity regularization model [START_REF] Bao | Error estimates of a regularized finite difference method for the logarithmic Schrödinger equation[END_REF][START_REF] Bao | Regularized numerical methods for the logarithmic Schrödinger equation[END_REF] or a local energy regularization approximation [START_REF] Bao | Error estimates of local energy regularization for the logarithmic Schrödinger equation[END_REF]. Note that even though it is rather natural to regularize the nonlinearity at least in the case of splitting methods (otherwise the solution of the ODE is singular), it may not be necessary to do so, as proved in [START_REF] Paraschis | On the convergence of the Crank-Nicolson method for the logarithmic Schrödinger equation[END_REF], in the case of the Crank-Nicolson method.

The model (1.3) cumulates several difficulties in terms of computation and numerical analysis. Most importantly, the strong dispersion caused by the repulsive harmonic potential requires to be handled very carefully as simulations are usually performed on a bounded domain. In the case of a logarithmic nonlinearity, the solution to (1.3) may or may not be dispersive, as recalled above. Moreover, the singularity of the logarithm at the origin makes the use of splitting methods delicate. In this paper, we propose a numerical strategy to face these two features, and extend it to the case where the logarithmic nonlinearity is replaced by a power nonlinearity. Our numerical methods are based on a generalized lens transform. This formulation is able to extract the main dispersion as well as oscillation and transforms the original equation into an equivalent one with non-dispersive solutions. This enables standard numerical methods to work successfully. The use of such transform in the context of numerical simulations appears to be one of the main novelties of this paper, and has possible extensions discussed in Remark 3.1 below. The other technique of our methods is to use a nonuniform temporal grid for the derived equivalent equation according to some designed rule consistent with the transform in the first step. This is the basis to establish the error estimates and numerical experiments show that it is superior than that by using a uniform grid directly. As pointed out in Remark 3.1, the approach presented here can be adapted to the case when one starts from a nonlinear Schrödinger equation without potential.

The rest of the article is organized as follows. Section 2 is devoted to recalling some properties, particularly the dispersion of the logNLS with repulsive potential (1.3). We introduce the generalized lens transform and present the splitting methods with a series of non-equidistant time steps in Section 3. This approach is extended to the nonlinear Schrödinger equation with repulsive potential and power nonlinearity in Section 4. Some error estimates of the time discretization are presented in Section 5. Extensive numerical experiments are displayed in Section 6 to show the accuracy and feasibility of the methods and investigate the dispersion properties of the equations with various parameters.

Some features of the logNLS under repulsive potential

In this section, we recall some properties of the dynamics of the logNLS under repulsive potential. A first unusual property associated to this logarithmic nonlinearity is that the size of the initial data plays no role (apart from a purely time dependent oscillation): if u solves (1.1), then for all k ∈ C, so does ku(t, x)e -itλ ln |k| 2 .

In the case of (1.3) (as well as (1.

2)), the potential V (x) is the sum of potentials depending only on one space variable, V (x) = d j=1

V j (x j ). The second property we emphasize is a tensorization property, which initially motivated the introduction of the logarithmic nonlinearity [START_REF] Bia Lynicki-Birula | Nonlinear wave mechanics[END_REF]: if the initial data is a tensor product,

u 0 (x) = d j=1 u 0j (x j ),
then the solution is given by

u(t, x) = d j=1 u j (t, x j ),
where each u j solves a one-dimensional equation,

i∂ t u j + 1 2 ∂ 2 xj u j = V j (x j )u j + λ ln |u j | 2 u j , u j (0, x j ) = u 0j .
Like the first property recalled above, this is reminiscent of linear Schrödinger equations, but effects caused by the logarithmic nonlinearity actually modify the dynamics, as we will see.

The third property, which is extremely convenient for numerical simulations, is that initial Gaussian data leads to a solution which remains Gaussian for all time. In view of the second property, we consider the case d = 1, and we seek the Gaussian solution with the form

(2.1) u(t, x) = b(t)e -a(t)x 2 /2 ,
particularly the initial data u 0 (x) = b(0)e -a(0)x 2 /2 (a(0) = α -iβ, α > 0) is Gaussian, for the equation (1.3). It can be easily checked that

(2.2) i ḃ = 1 2 ab + λb ln |b| 2 , i ȧ = a 2 + 2λRe a + ω 2 ,
where ḃ is the derivative with respect to time. Writing a as

(2.3) a = 1 µ 2 -i μ µ , µ ∈ R + , leads to (2.4) μ = 2λ µ + 1 µ 3 + ω 2 µ, µ(0) = 1 √ α , μ(0) = β √ α .
Plugging (2.3) into (2.2) and integrating in time yields

(2.5) b(t) = b(0)e iθ(t) µ(0)/µ(t), θ(t) ∈ R.
Multiplying (2.4) by μ and integrating, we get

(2.6) ( μ) 2 = C 0 + 4λ ln µ - 1 µ 2 + ω 2 µ 2 ,
where C 0 = μ(0) 2 -4λ ln µ(0)+ 1 µ(0) 2 -ω 2 µ(0) 2 is related to the initial data. Noticing that F (q) = C 0 + 4λ ln q -1 q 2 + ω 2 q 2 → -∞ when q → 0, this implies that µ is bounded away from zero, i.e., ∃δ > 0, µ(t) δ, ∀t 0. Combining (2.1), (2.3) and (2.5), one arrives at

|u(t, x)| = |b(0)| µ(t)α 1/2 e -x 2 2µ(t) 2 .
Hence it suffices to study the dynamics of µ(t), i.e., the ODE (2.4). Concerning this, we have:

Proposition 2.1 (Propagation of Gaussian data, [START_REF] Carles | Nonuniqueness and nonlinear instability of Gaussons under repulsive harmonic potential[END_REF]). Let d = 1, λ < 0 < ω.

1. If λ < -ω, then (2.4) has exactly two stationary solutions, µ ± = 1/ k ± with k ± = -λ± √ λ 2 -ω 2 , which correspond to the positive stationary solutions of (1.3):

φ ± (x) = exp -k± 4λ -k±x 2 2
and generate a continuous family of solitary waves,

u ±,ν (t, x) = φ ±,ν (x)e iνt , φ ±,ν (x) = e -ν 2λ φ ± (x), ν ∈ R.
The other solutions to (2.4) are either periodic or unbounded, which correspond to time-periodic or dispersive Gaussian solutions to (1.3), respectively.

2. If λ = -ω, then (2.4) has exactly one stationary solution, µ 0 = 1/ √ ω. All the other solutions are unbounded. In other words, any Gaussian solution to (1.3) which is not of the form e iθ e (2ν+ω)/(4ω) e iνt e -ωx 2 /2 , ν ∈ R, θ ∈ [0, 2π], is dispersive. 3. If λ > -ω, then every solution to (2.4) is unbounded and e ωt µ(t) e ωt . This implies every Gaussian solution to (1.3) disperses exponentially in time.

While for λ > 0, it can be easily observed from (2.2) that there are no stationary Gaussian solutions. Furthermore, it was shown in [START_REF] Carles | Orbital stability of Gausson solutions to logarithmic Schrödinger equations[END_REF]Proposition 1.8] that for general initial data u 0 ∈ Σ, the solution of (1.3) can be scaled as

u(t, x) = 1 ν(t) 1/2 v t, x ν(t) exp i ν(t) ν(t) x 2 2 ,
where v is bounded and non-dispersive, and

ν(t) ∼ ν ∞ e ωt , ν(t) ∼ ων ∞ e ωt , t → ∞.

Splitting numerical methods based on a transformation

It can be seen from the above section that for λ < 0 and Gaussian initial data, the solution of (1.3) disperses very fast, except some special solutions, e.g., standing waves or time-periodic solutions. While for λ > 0 and general initial data, each solution disperses exponentially in time. On the one hand, we are curious whether the solution is dispersive in the case λ < 0, for general initial data. One the other hand, possible quick dispersion makes commonly used naive truncation by homogeneous Dirichlet or periodic boundary conditions infeasible in practical computation. In this section we introduce a transformation which extracts the main dispersion such that the simple truncation technique works well. Then a classical time splitting discretization is presented.

For simplicity of notation, we consider d = 1 and extensions to higher dimensions are straightforward, provided that the potential is isotropic (like in (1.3)), as generalized lens transforms are available then (see [START_REF] Carles | Nonlinear Schrödinger equations with repulsive harmonic potential and applications[END_REF]). This fact is independent of the tensorization property related to the logarithmic nonlinearity pointed out in Section 2. The lens transform seems to have appeared in [START_REF] Niederer | The maximal kinematical invariance groups of the harmonic oscillator[END_REF], where the author investigated the largest group of coordinate transformations leaving invariant the Schrödinger equation with harmonic potential. It is a general fact that for isotropic quadratic potentials, possibly time dependent, similar formulas are available, see [START_REF] Carles | Nonlinear Schr dinger equation with time dependent potential[END_REF]. In the case of the (time independent) repulsive harmonic potential, the formula appears in [START_REF] Carles | Nonlinear Schrödinger equations with repulsive harmonic potential and applications[END_REF]: u is changed to v via the identity

(3.1) u(t, x) = 1 cosh(ωt) v (s, y) e i ω 2 x 2 tanh(ωt) , s = tanh(ωt) ω , y = x cosh(ωt)
.

Plugging this formula into the equation (1.3) leads to

0 = i∂ t u + 1 2 ∂ xx u + ω 2 2 |x| 2 u -λu ln |u| 2 = e i ω 2 x 2 tanh(ωt) cosh(ωt) i∂ s v sech 2 (ωt) + 1 2 sech 2 (ωt)∂ yy v -λv ln |u| 2 = e i ω 2 x 2 tanh(ωt) cosh 5/2 (ωt) i∂ s v + 1 2 ∂ yy v -λ cosh 2 (ωt)v ln |v| 2 -ln(cosh(ωt)) .
Noticing that cosh 2 (ωt) =

1 1-tanh 2 (ωt) = 1 1-ω 2 s 2 , we see that v solves the equation    i∂ s v + 1 2 ∂ yy v = λv 1 -ω 2 s 2 ln |v| 2 + 1 2 ln(1 -ω 2 s 2 ) , (s, y) ∈ (0, 1/ω) × R, v(0, y) = u(0, y) = u 0 (y).
We emphasize that the time interval t ∈ [0, ∞) has been transformed into s ∈ [0, 1/ω): on such bounded time interval, dispersive effects will not be too strong, also because we have filtered out the maximum (exponential) dispersive effects by changing the space variable. On the other hand, the equation is now nonautonomous, as a time dependent factor has appeared in front of the nonlinearity.

In the case of splitting methods, this will require to compute various integrals in the time variable, when solving the ODE part. The last term above involves a purely time-dependent potential, and can be absorbed by a gauge transform. Set

(3.2) v(s, y) = κ(s, y)e -iλg(s) , g(s) = 1 2 s 0 ln(1 -ω 2 p 2 ) 1 -ω 2 p 2 dp. Then κ satisfies the equation (3.3)    i∂ s κ + 1 2 ∂ yy κ = λκ 1 -ω 2 s 2 ln |κ| 2 , (s, y) ∈ (0, 1/ω) × R, κ(0, y) = u(0, y) = u 0 (y).
Hence investigating the large-time behavior of u(t, x) reduces to studying the asymptotic dynamics of κ when s approaches 1/ω. Noticing (3.1), for each dispersive solution u(t, x), the coefficient in the exponential function in space is of order e -ωt , hence it requires large computational bounded domain increasing exponentially in time. Moreover, due to the term exp i ω 2 x 2 tanh(ωt) , which is highly oscillatory in space, it requires very tiny mesh size in practical computation and brings difficulties for dynamics.

To simulate the dynamics of (3.3), firstly we regularize (3.3) by introducing a small parameter 0 < ε 1 and truncate the problem on a bounded computational domain Ω = (a, b) with periodic boundary condition:

(3.4)    i∂ s κ ε + 1 2 ∂ yy κ ε = 2λκ ε 1 -ω 2 s 2 ln(|κ ε | + ε), (s, y) ∈ (0, 1/ω) × Ω, κ ε (0, y) = u(0, y) = u 0 (y).
We use the splitting method to solve the equation (3.4), based on the splitting

∂ s φ = A(φ) + B(s, φ),
where

A(φ) = i 2 ∂ yy φ, B(s, φ) = - 2iλφ 1 -ω 2 s 2 ln(|φ| + ε),
and the solutions of the subproblems

(3.5) ∂ s φ(s n + s, y) = A(φ), y ∈ Ω, s > 0, φ(s n , y) = φ n (y), (3.6) ∂ s z(s n + s, y) = B(s n + s, z), y ∈ Ω, s > 0, z(s n , y) = z n (y).
The associated evolution operators are given by

φ(s n + s, •) = Φ s A (φ n ) = exp is 2 ∂ yy φ n , z(s n + s, •) = Φ s B (s n , z n ) = z n exp -2iλ ln(|z n | + ε) s 0 1 1 -ω 2 (s n + p) 2 dp = z n (|z n | + ε) -iλ ω ln( (1-ωsn )(1+ω(sn+s)) (1+ωsn)(1-ω(sn+s)) ) .
For fixed T > 0 and the number of time steps N ∈ N * , we denote the time step τ = T /N > 0, t n = nτ , s n = tanh(ωtn) ω for n ∈ N * and set

(3.7) δ n = s n+1 -s n .
We consider the Lie-Trotter splitting for solving (3.4) by using a series of nonequidistant time steps {δ n }:

(3.8) κ n+1 ε = Φ L (κ n ε ) = Φ δn A (Φ δn B (s n , κ n ε )
), κ 0 ε = u 0 , and the Strang splitting scheme:

(3.9) κ n+1 ε = Φ S (κ n ε ) = Φ δn/2 A (Φ δn B (s n , Φ δn/2 A (κ n ε ))), κ 0 ε = u 0 .
With the value of κ n ε , we can approximate u(t n , x) via (3.1) by setting

(3.10) u n ε (x) = 1 cosh(ωt n ) exp i ω 2 x 2 tanh(ωt n ) -iλg(s n ) κ n ε x cosh(ωt n ) ,
where g is defined as (3.2). The above scheme (3.10) together with (3.8) or (3.9) displays a time integrator (semi-discretization) for solving (1.3). In practical computation, we combine it with the Fourier pseudospectral method for spatial discretization [START_REF] Shen | Spectral Methods: Algorithms, Analysis and Applications[END_REF]. The schemes are explicit and efficient thanks to the fast Fourier transform (FFT) for calculating (3.8) or (3.9) and the nonuniform fast Fourier transform (NUFFT) [START_REF] Greengard | Accelerating the nonuniform fast Fourier transform[END_REF][START_REF] Jiang | Fast and accurate evaluation of nonlocal Coulomb and dipole-dipole interactions via the nonuniform FFT[END_REF] for computing (3.10). That is, we get the value of κ n ε on a fixed uniform mesh, or equivalently the Fourier coefficients in phase space by standard FFT. To approximate the value of u n ε on a uniform mesh, we need to calculate the value of κ n ε at the points x cosh(ωtn) by recalling (3.10), which varies with respect to time. This can be realised by NUFFT. For the details, we refer to [START_REF] Greengard | Accelerating the nonuniform fast Fourier transform[END_REF][START_REF] Jiang | Fast and accurate evaluation of nonlocal Coulomb and dipole-dipole interactions via the nonuniform FFT[END_REF]. Remark 3.1 (Case without potential: lens transform). In the case of a nonlinear Schrödinger equation without potential,

(3.11) i∂ t u + 1 2 ∆u = f |u| 2 u, (t, x) ∈ (0, ∞) × R d ,
one may introduce the lens transform

v(s, y) = 1 (cos(ωs)) d/2 u(t, x)e -i ω 2 |y| 2 tan(ωs) , t = tan(ωs) ω , x = y cos(ωs)
.

Note that the relation between u and v is somehow reversed compared to (3.1). The study of (3.11) is equivalent to the study of

i∂ s v + 1 2 ∆v = ω 2 2 |y| 2 v + 1 cos 2 (ωs) f |v| 2 cos d (ωs) v, (s, y) ∈ (0, π 2ω ) × R d .
Like above, the time interval is now bounded, and the presence of the (confining) harmonic potential makes it possible to avoid boundary effects. This provides an alternative approach to compute the scattering operator associated to (3.11), compared to [START_REF] Carles | Numerical aspects of nonlinear Schrödinger equations in the presence of caustics[END_REF]. Note that in the case of an L 2 -critical nonlinearity, f (|z| 2 ) = λ|z| 4/d , the nonlinear term reduces to λ|v| 4/d v, and the equation in v is autonomous [START_REF] Carles | Critical nonlinear Schrödinger equations with and without harmonic potential[END_REF].

Power nonlinearity case

For a comparison, next we consider the nonlinear Schrödinger equation with repulsive potential and power nonlinearity. For simplicity, we consider the onedimensional case only (note that with a power nonlinearity, the tensorization property recalled in Section 2 is lost):

(4.1)    i∂ t u + 1 2 ∂ xx u = -ω 2 x 2 2 u + λ|u| 2σ u, x ∈ R, u(0, x) = u 0 (x),
with σ > 0. Note that in the case of a power nonlinearity, initial Gaussian data do not propagate as Gaussians (this can be seen by seeking solutions under the form of a Gaussian function). The Cauchy problem and scattering theory have been investigated in [START_REF] Carles | Nonlinear Schrödinger equations with repulsive harmonic potential and applications[END_REF]. Among other results, it is proven that in the case λ > 0 (defocusing case, considered in the simulations in Section 6), every power nonlinearity is short range, in the sense that the (nonlinear) solution behaves like a solution to the linear equation ((4.1) with λ = 0): for every

u 0 ∈ Σ, (4.2) ∃u + ∈ Σ, e -i t 2 (∂xx+ω 2 x 2 ) u(t) -u + Σ = O e -σω
e i t 2 (∂xx+ω 2 x 2 ) u + = e -iπ 4 ω sinh(ωt) F u + e i ω 2 coth(ωt)|•| 2 ωx sinh(ωt) e iω coth(ωt) x 2 2 ∼ t→+∞ e -iπ 4 ω sinh(ωt) F u + e i ω 2 |•| 2 ωx sinh(ωt) e iω coth(ωt) x 2 2 ,
where we normalize the Fourier transform as

F(f )(ξ) = f (ξ) = 1 √ 2π R e -ixξ f (x)dx.
While for the focusing case λ < 0, finite-time blow-up might occur under some assumptions on λ, ω, σ and u 0 (e.g. [19, Theorem 1.1]), that is, there exists T > 0 such that

(4.3) lim t→T ∇u(t) L 2 = +∞.
We study the dispersion properties of (4.1) numerically. With the same formulation (3.1), we get the equation for v as (4.4)

   i∂ s v + 1 2 ∂ yy v = λ 1 -ω 2 s 2 σ 2 -1 |v| 2σ v, (s, y) ∈ (0, 1/ω) × R, v(0, y) = u(0, y) = u 0 (y).
Denote Φ s B by the evolution operator of the nonlinear subproblem

(4.5) ∂ s z(s n + s, y) = -iλ(1 -ω 2 (s n + s) 2 ) σ 2 -1 |z| 2σ z, s > 0, y ∈ R, z(s n , y) = z n (y).
Then it can be written explicitly as

z(s n + s, y) = Φ s B (s n , z n ) = z n exp -iλ|z n | 2σ sn+s sn (1 -ω 2 ζ 2 ) σ 2 -1 dζ .
Using the same notations as before, i.e., s n = tanh(ωtn) ω with t n = nτ , and

δ n = s n+1 -s n , we get (4.6) z(s n+1 , y) = Φ δn B (s n , z n ) = z n exp -iλ|z n | 2σ tn+1 tn sech σ (ωζ)dζ .
Combining this with the evolution operator of the linear subequation Φ s A (3.5), we can get the first-order Lie-Trotter splitting and the second-order Strang splitting approximations as (4.7)

v n+1 = Φ L (v n ) = Φ δn A (Φ δn B (s n , v n )), v 0 = u 0 ; and (4.8) v n+1 = Φ S (v n ) = Φ δn/2 A (Φ δn B (s n , Φ δn/2 A (v n ))), v 0 = u 0 ,
respectively. While the approximation u n+1 can be recovered by (4.9)

u n+1 (x) = 1 cosh(ωt n+1 ) exp i ω 2 x 2 tanh(ωt n+1 ) v n+1 x cosh(ωt n+1 )
.

In practical computation, we combine the splitting integrator with Fourier pseudospectral discretization in space and a highly accurate quadrature rule for approximating the integral in (4.6).

On error estimates for the splitting methods

In both cases, (1.3) and (4.1), we have used the generalized lens transform (3.1) to turn the original equation into a nonlinear Schrödinger equation without potential, with a nonautonomous nonlinearity, (3.3) and (4.4), respectively. It is possible to rely on error estimates which have been established for splitting methods applied to such equations. Three differences must be taken into account though:

• A time dependent factor has appeared in front of the nonlinearity (expect in the very special L 2 -critical case σ = 2/d in (4.1)).

• The time step is uniform in terms of the t variable, hence not in terms of the s variable. • The factor 1 1-ω 2 s 2 is singular, and not integrable, as s → 1/ω. In the case of nonlinear Schrödinger equations without potential ((3.3) and (4.4) without the time dependent factor in front of the nonlinearity), error estimates for Lie-Trotter time splitting schemes are available with initial data in H 1 : see [START_REF] Bao | Regularized numerical methods for the logarithmic Schrödinger equation[END_REF] for the case of a (regularized) logarithmic nonlinearity, and [START_REF] Choi | On the splitting method for the nonlinear Schrödinger equation with initial data in H 1[END_REF][START_REF] Ignat | A splitting method for the nonlinear Schrödinger equation[END_REF] for the power nonlinearity. Throughout this section, we assume that the numerical solution is given by a Lie-Trotter splitting scheme, as opposed to the Strang splitting scheme considered so far: indeed error estimates for Strang splitting scheme in the presence of a logarithmic nonlinearity are far less satisfactory, see [START_REF] Bao | Regularized numerical methods for the logarithmic Schrödinger equation[END_REF]Remark 4].

In view of (3.1), applying the gradient ∇ y to the solution of (3.3) and (4.4) amounts to considering the action of the vector-field J on the initial unknown u, where J is defined by

J(t) = ωx sinh(ωt) + i cosh(ωt)∇ = e -itH (i∇) e -itH = i cosh(ωt)e iω tanh(ωt) |x| 2 2 ∇ e -iω tanh(ωt) |x| 2 2
• ,

where

H = - 1 2 ∆ -ω 2 |x| 2 2 .
This vector-field plays a central role in the proofs in [START_REF] Carles | Nonlinear Schrödinger equations with repulsive harmonic potential and applications[END_REF].

The splitting considered for v (or, equivalently, κ), is equivalent to a splitting for u, thanks to the definition of s n ,

s n = tanh(ωt n ) ω , t n = nτ.
Indeed, the case λ = 0 shows that (3.1) maps a solution u of the linear Schrödinger equation with potential to a solution v of the linear Schrödinger equation without potential, with u |t=0 = v |s=0 . We readily check that the same holds regarding the ordinary differential equations. The last formulation for J shows that

J(t)u L 2 (R d ) = ∇v(s) L 2 (R d )
, and we see that it is equivalent to consider splitting methods,

• For v with an analysis in H 1 , that is, involving v and ∇v, with e i s 2 ∆ as the linear operator;

• For u with an analysis involving u and Ju, with e -itH as the linear operator. 5.1. Logarithmic nonlinearity. In [START_REF] Bao | Regularized numerical methods for the logarithmic Schrödinger equation[END_REF], an error estimate is proven with initial data in H 1 (in the 1D case -data in H 2 if d = 2 or 3), on any bounded time interval: the constants given by the proof grow at least exponentially in time. By resuming the same proof step by step, following either of the two strategies described above, we can prove that for any T > 0, there exists ε 0 > 0 such that when 0 < ε ε 0 and 0 nτ T , we have:

u ε (t n ) -u n ε L 2 (R) = κ ε (s n ) -κ n ε L 2 (R) C (T, u 0 Σ ) ln(ε -1 )τ 1/2 ,
where C(•, •) is independent of ε ∈ (0, ε 0 ], and u ε is related to κ ε via the same formula as the one relating u to κ.

To compare the numerical solution with the exact solution u, we then have to compare u with u ε , or, equivalently in L 2 , κ with κ ε . We first note that [8, Lemma 2.3] remains unchanged in the presence of a time dependent factor in front of the logarithm, but we have to stick to time intervals where this factor remains bounded (or equivalently, bounded time intervals for u): for S ∈ [0, 1/ω), there exists C = C(S) such that for all s ∈ [0, S],

d ds κ ε (s) -κ(s) 2 L 2 (R) C κ ε (s) -κ(s) 2 L 2 (R) + ε κ ε (s) -κ(s) 2 L 1 (R) .
Note that the transform (3.1) does not preserve the L 1 -norm in space (L 2 is the only Lebesgue norm which is preserved), but for bounded t ∈ [0, T ] (that is, for s ∈ [0, S] and S = tanh(ωT ) ω ∈ [0, 1/ω)), there exists C(T ) such that for all t ∈ [0, T ],

1 C(T ) κ ε (s) -κ(s) L 1 (R) u ε (t) -u(t) L 1 (R) C(T ) κ ε (s) -κ(s) L 1 (R) .
Modifying e.g. the proof of [8, Lemma 2.6], we readily prove: there exists C such that for any v ∈ Σ,

(5.1) v L 1 (R) C v 1/2 L 2 (R) v 1/2
Σ . In view of [21, Proposition 1.3], we infer that there exists C(T ) such that

κ ε (s) -κ(s) Σ C(T ), ∀s ∈ [0, S], S = tanh(ωT ) ω ,
and so

d ds κ ε (s) -κ(s) 2 L 2 (R) C(T ) κ ε (s) -κ(s) 2 L 2 (R) + ε κ ε (s) -κ(s) L 2 (R) .
Using Gronwall Lemma, we conclude:

Proposition 5.1. Suppose d = 1, and consider the Lie-Trotter splitting method (3.8) combined with (3.10). Let T > 0 and u 0 ∈ Σ. There exists ε 0 > 0 such that when 0 < ε ε 0 and 0 nτ T , we have

u(t n ) -u n ε L 2 (R) = κ(s n ) -κ n ε L 2 (R) C (T, u 0 Σ ) ln(ε -1 )τ 1/2 + ε , where C(•, •) is independent of ε ∈ (0, ε 0 ].

Power nonlinearity.

In the case of a power nonlinearity, it is not necessary to regularize the nonlinearity, but for σ < 1/2, the nonlinearity in (4.1) is not C 2 : working with an H 2 regularity in space may be delicate (one way to overcome the lack of regularity might be to consider one time derivative, and use the equation to infer H 2 regularity in space, see [START_REF] Cazenave | Semilinear Schrödinger equations[END_REF]). In [START_REF] Choi | On the splitting method for the nonlinear Schrödinger equation with initial data in H 1[END_REF][START_REF] Ignat | A splitting method for the nonlinear Schrödinger equation[END_REF], discrete Strichartz inequalities were used in order to decrease the required regularity, to H 1 (still in the case of Lie-Trotter method). However, as noted in [START_REF] Ignat | Dispersive properties of numerical schemes for nonlinear Schrödinger equations[END_REF][START_REF] Ignat | Numerical dispersive schemes for the nonlinear Schrödinger equation[END_REF], in the absence of potential, it is necessary to introduce a frequency cut-off in the free propagator e i t 2 ∆ , and consider instead e i t 2 ∆ Π τ , where Π τ is a Fourier multiplier of symbol χ(τ 1/2 ξ), where χ is (sufficiently) smooth and compactly supported: this corresponds to a modification of the operator A as defined in Section 3. This implies that if one wants to prove error estimates based on discrete Strichartz inequalities like in [START_REF] Choi | On the splitting method for the nonlinear Schrödinger equation with initial data in H 1[END_REF][START_REF] Ignat | A splitting method for the nonlinear Schrödinger equation[END_REF], in the case of (4.1), one should consider the same frequency cut-off when working with the unknown v. We note that the frequency cut-off amounts to imposing the frequency localization |ξ| τ -1/2 .

In the case of (4.4), the time step is not uniform, but

|s n+1 -s n | = tanh(ω(n + 1)τ ) ω - tanh(ω nτ ) ω τ,
and so |ξ| τ -1/2 implies |ξ| |s n+1 -s n | -1/2 uniformly in n.

For s ∈ [0, S] and S ∈ [0, 1/ω), the time dependent factor in front of the nonlinearity in (4.4) is harmless, and the proofs of [START_REF] Choi | On the splitting method for the nonlinear Schrödinger equation with initial data in H 1[END_REF] can be repeated, in order to obtain:

Proposition 5.2. Let d = 1, σ > 0, u 0 ∈ Σ and T > 0. In (4.7), replace the operator Φ s A with Φ s à = Φ s A Π τ , where Π τ φ(ξ) = χ(τ 1/2 ξ) φ(ξ), and χ ∈ C ∞ (R) is a cut-off function supported in [-2, 2] such that χ ≡ 1 on [-1, 1]. There exists C = C(σ, T, u 0 Σ ) such that max 0 nτ T u n -u(t n ) L 2 (R) = max 0 sn S v n -v(s n ) L 2 (R) Cτ 1/2 ,
where S = tanh(ωT ) ω . Remark 5.1. In view of the analysis performed in [START_REF] Carles | Scattering and uniform in time error estimates for splitting method in NLS[END_REF], it is likely that in the case λ > 0, the above estimate is uniform in time, that is, C can be chosen independent of T . Indeed, when working with the unknown u, the vector-field J provides some exponential decay in time of various Lebesgue norms in space, as exploited in [START_REF] Carles | Nonlinear Schrödinger equations with repulsive harmonic potential and applications[END_REF] to prove scattering results. The proof in [START_REF] Carles | Scattering and uniform in time error estimates for splitting method in NLS[END_REF] being already quite technical, we do not explore the details of the argument here.

Remark 5.2. In the case σ > 1/2, the nonlinearity in (4.1) is C 2 , and considering initial data in H 2 ∩F(H 2 ), the above error rate τ 1/2 can be improved to τ , following the approach of [START_REF] Besse | Order estimates in time of splitting methods for the nonlinear Schrödinger equation[END_REF][START_REF] Lubich | On splitting methods for Schrödinger-Poisson and cubic nonlinear Schrödinger equations[END_REF].

Numerical Results

In this section, we first test the order of accuracy of the proposed method (3.10) combined with the Lie-Trotter splitting (3.8) or the Strang splitting (3.9). Then we apply the Strang-splitting method to investigate some long time dynamics of the nonlinear Schrödinger equation with repulsive potential and logarithmic nonlinearity (1.3) (or equivalently (3.4)) or power nonlinearity (4.1) (or equivalently (4.5)). Example 1. Here, we set λ = -3, ω = 2 in the equation (1.3) and ε = 10 -15 in the regularized model (3.4). Choose the initial data as

u 0 (x) = Ae -αx 2 /2 with α = -λ - √ λ 2 -ω 2 = 3 - √ 5.
Then according to Proposition 2.1, it generates the solitary wave as u(t, x) = Ae iνt e -αx 2 /2 , with ν = -(2λ ln(A) + α/2). Set A = 2, then the exact solution is given by (6.1) u(t, x) = 2e -αx 2 /2 e -i(α/2-6 ln 2)t .

We set Ω = (-L, L) with L = 10 in (3.4). To quantify the numerical error, we define the error function as

e(t n ) = u n ε (x) -u(t n , x),
where u n ε is obtained by (3.10) combined with (3.8) or (3.9). In practical computation, we compute until T = 2.5 by using different values of N (the number of time steps) and a fixed mesh size h = 1/2 9 for spectral spatial discretization, which is small enough for neglecting the spatial error introduced. The time integral in g(s) is approximated by Simpson rule with very fine mesh such that the error introduced by numerical integration is ignorable.

Before presenting the numerical results, we recall the temporal grid introduced in Section 3. We notice that the computation is performed by utilizing a series of non-equidistant time steps δ n defined by (3.7), which results as a series of numerical solutions κ n ε as an approximation for κ on a nonuniform grid {s n = tanh(ωtn) ω } and a series of solutions u n via (3.10) as an approximation for u on a uniform grid {t n = nT N }. We denote the corresponding methods as Lie I ((3.8)), and Strang I ((3.9)), respectively, or Lie-t and Strang-t, which means the temporal grid is uniform in variable t. On the other hand, another intuitive approach is to solve the equation (3.4) directly by a standard temporal discretization: we set the step size δ n ≡ δ := s N N with s N = tanh(ωt N ) ω . Then we get a series of solutions κ n ε as an approximation for κ on a uniform grid {s n = nδ} and a series of solutions u n via (3.10) as an approximation for u on a nonuniform grid {t n = artanh(ωsn) ω }. For convenience afterwards, we denote the corresponding methods as Lie II ((3.8)) and Strang II ((3.9)), or Lie-s and Strang-s, which means the temporal grid is uniform in variable s. It can be clearly seen that via both cases of temporal grid, the Lie-Trotter splitting method converges at the first order in time. The right plot shows the evolution of the error when the step number N is fixed as N = 25000. We observe that the error increases exponentially in time for both cases and the error of Lie I increases a bit slower than that of Lie II, which suggests to choose Lie I for long time simulation. Figure 6.2 shows the error e(2.5) L 2 for the methods Strang I and Strang II. It can be clearly observed that the splitting scheme converges at the second order in time and the error of Strang I is far less than that of Strang II. Moreover, for the method Strang II, it leads to a correct solution and is convergent quadratically only when the step number N is large enough N N T . Figure 6.3 depicts the evolution of the error obtained via Strang I and Strang II with fixed step number N = 25000. We are surprised to find that the error of Strang I method ncreases much more slowly than that of Strang II, the latter of which increases exponentially in time by noticing that the longitudinal axis represents the error (left plot in Figure 6.3) and the logarithm of error (right plot in Figure 6.3), respectively. The excellent convergence behavior of Strang I is different from that of Lie I and we can benefit a lot in practical computation, especially for long time simulation. We postpone the analytical study of this method to a future work. -6 -3 0 3 6

x -2 -1 0 1 2 Im(u(t, x)) t = 1 t = 2 t = 3 -1.5 -0.5 0.5 1.5 y -2 -1 0 1 2
Im( κε(s=tanh(ωt)/ω,y) For the solitary wave, it is easily checked that direct classical methods for u work well since dispersion does not occur. On the contrary, the transform (3.1) (and (3.2)) makes the "essential" support of the solution v (and κ) shrink exponentially with respect to time, as is shown in Figure 6.4. This requires that the mesh size h has to be set tinier and tinier as time evolves in order to capture the correct solution, which introduces huge costs if one keeps the computational domain fixed. However, this difficulty can be overcome by shrinking the computational region correspondingly in time, which enables the computational cost to be comparable with the direct splitting methods for (1.3) (regularized version).

cosh 1/2 (ωt) ) t = 1 t = 2 t = 3 -0.02 0 0.02
This suggests that a reasonable way to proceed may be to first simulate (1.3) directly: if the solution is not dispersive, then nothing specific is needed (apart from the regularization of the logarithm). If the solution is dispersive (which is always the case when λ > 0), then boundary effects become strong, and it is more efficient to consider (3.4) and use Strang I method via (3.9)-(3.10) on a nonuniform grid {s n = tanh(ωtn) ω }.

Example 2. We set λ = -3 and ω = 2 in the equation (1.3). We consider the following several cases: (i). u 0 (x) = 2e -αx 2 /2 for α = 2; (ii). u 0 (x) = 2e -αx 2 /2 for α = 1/2; (iii). u 0 (x) = sech(x 2 /2).

According to Proposition 2.1, for λ = -3, ω = 2, there exist two stationary solutions. The ODE (2.4)

μ = 1 µ 2 -k - 1 µ 2 -k + µ, k ± = -λ ± λ 2 -ω 2 ,
produces different trajectories due to different initial data for µ, which is illustrated in Figure 6.5 displaying the phase portrait for the equation (2.4). For Case (i), α = 2 corresponds to µ(0) = 1/ √ α = 1/ √ 2, μ(0) = 0, which lies in the region producing time periodic solution µ(t), as seen from Figure 6.5. This is confirmed by the dynamics shown in Figure 6.6. We remark here that Figure 6.6 is obtained by the standard Strang splitting method for (1.3) directly to avoid necessary treatments including adaptive mesh refinement and computational domain cut, since there is no dispersion for u. Figure 6.7 shows the dynamics of κ ε (s, y) for Case (ii) initial data by using Strang I method till T = 4, or equivalently s = tanh(4ω)/ω with time step τ = 0.0001 and mesh size h = 1/2 10 . It can be observed that κ ε , or accordingly κ, is well localized in space and |κ(s, y)| stays almost invariant after some time (cf. the top left one in Figure 6.7), however, the argument of κ varies in time (the top right one in Figure 6.7). Moreover, κ oscillates in space. The bottom part in Figure 6.7 shows the profile of u(t, x) obtained by (3.10), which shows the dispersion of u. The right one displays Re(u(4, x)), which shows that u is highly oscillatory in space. The oscillation is more and more drastic when observed further away from the origin, which agrees with the transformation (3.10), due to the term exp i ω 2 x 2 tanh(ωt n ) . Because of the rapid expansion and drastic oscillation in space, it is a disaster if one solves (1.3) directly. Figure 6.8 shows the profile of κ(s, y) for Case (iii) initial data. Different from Case (ii) where the modulus of κ keeps invariant after some time, here the modulus keeps increasing with respect to time, which is illustrated in the top graphics. The wave is split into two symmetric wave packets and a wave in between. As time evolves, the width of the two wave packets gets narrower and narrower. With the formulation (3.10), we get the dynamics of u(t, x) in the bottom in Figure 6.8. We observe that the initial wave is split into two symmetric breathers which spread out quickly (exponentially) in space with modulus varying periodically in time and a dispersive wave in between. Figure 6.9 displays the profile of Re(κ) and Re(u) at t = 3, which suggests the feasibility of solving u via the formulation (3.10) and the virtual difficulty of solving u directly. Example 3. We set λ = -2 and ω = 2 in the equation (1.3), in which case there exists a continuous family of solitary wave when u 0 (x) = Ae -αx 2 /2 for α = 2 (or, equivalently µ(0) = 1/ √ 2), according to Proposition 2.1, that can be illustrated from the phase portrait for the ODE (2.4) (cf. Figure 6.10). We consider the following three cases for the initial data u 0 : (i). u 0 (x) = 2e -αx 2 /2 for α = 3; (ii). u 0 (x) = 2e -αx 2 /2 for α = 1/2; (iii). u 0 (x) = sech(x 2 /2). for T = 5 by using the method Strang I with τ = 0.0002 and h = 1/2 12 on a computational domain Ω = (-10, 10) with periodic boundary condition. We observe that for all three cases, the modulus |κ| remain invariant after some time while the variation of the argument introduces the oscillation of κ in space. This suggests that |u| decreases exponentially in time, with rate e -ωt/2 in view of (3.10). Moreover, the argument of κ introduces more and more oscillation in space as time involves (cf. the right column in Figure 6.11).

Example 4. We set λ = 1 and ω = 2 in the equation (1.3), in which case there exists no stationary solution for equation (2.4). We consider the following two cases for the initial data u 0 : Figure 6.12 shows the dynamics of κ until s = tanh(ωT ) ω for T = 5 with τ = 0.0002 and h = 1/2 12 on a computational domain Ω = (-20, 20) with periodic boundary condition. We observe that the behavior is similar to that in Example 3, the modulus can attain the invariant state soon and more and more mild oscillations are created as time evolves. With the formulation (3.1), the dispersion and oscillation properties of u can be clearly revealed. Next we display some numerical experiments for dynamics of the solution v for the equation (4.4), or accordingly the dynamics of u (through (4.9)) for the power nonlinearity case (4.1). We still use Strang I method for temporal discretization with τ = 0.0002 and h = 1/2 10 as the spatial mesh size on a computational domain Ω = (-20, 20) with periodic boundary condition. Figures 6.13 and 6.14 display the dynamics of v for Cases (i) and (ii) initial data, respectively. We observe that for both cases and all chosen σ = 1, 2, 3, both the modulus |v| and the argument varies very slowly (almost invariant) after some time, which is different from the case with logarithmic nonlinearity where the argument varies quickly with respect to time. For Case (i) initial data, the "essential" support of the (almost) invariant state gets larger (cf. left column in Figure 6.13) and more and more oscillation is created in space (cf. right column in Figure 6.13) when σ increases. While for Case (ii) initial data, one can hardly tell any remarkable difference for various choices of σ, as shown in Figure 6.14.

(i). u 0 (x) = sech(x 2 /2) sin(x); (ii). u 0 (x) = 2e -(x-3) 2 /4 + 2e -(x+3) 2 .
In the above examples, the numerical solution of (4.4) stays almost invariant after some time. This is understandable by noticing that the exponent in (4.6) gets very small as t n increases, and due to the rapid convergence toward the linear dynamics, As was shown in [START_REF] Carles | Nonlinear Schrödinger equations with repulsive harmonic potential and applications[END_REF], finite-time blow-up might occur when the parameters λ, ω, σ and u 0 satisfy some condition. Here we choose σ = 1, 2, 3, respectively, and solve the equation (4.5) by Strang I method. For σ = 1, the dynamics of v is shown in Figure 6.15. Similar results are obtained as in the case λ > 0, which suggests that in this case, the solution in u has the same dispersion as before, i.e., it disperses exponentially in time at the rate e ωt/2 . While for σ = 3, the dynamics of v is displayed in Figure 6.16, from which we observe that the solution concentrates around the origin quickly. This is in agreement with the virial computation leading to [19, Theorem 1.1, case 3]: in the present case we check that 1 2 ∇u 0 2

L 2 + λ σ + 1 u 0 2σ+2 L 2σ+2 < - ω 2 2 xu 0 2 L 2 -ω Im u 0 x • ∇ x u 0 ,
as every term can be computed in the Gaussian case, and so the solution blows up in the future and in the past. The right plot shows the dynamics of ∇v, which suggests that blow-up (4.3) occurs at around T = 0.02445. The situation is similar for σ = 2, which is omitted here for brevity. In one word, the proposed method based on the formulation (3.1) is able to capture the correct dynamics for both dispersive solutions and finite-time blow-up solutions. It is superior than that by solving the original equation (4.1) directly when the finite-time blow-up does not occur so soon.

Conclusion

We proposed some time splitting methods for the nonlinear Schrödinger equation with a logarithmic nonlinearity and a repulsive potential based on the generalized lens transform. This transformation can capture the main dispersion and oscillation in the solution. It neutralizes the possible boundary effects and enables the classical splitting methods to work smoothly for the equivalent formulation. This approach was extended to the case with a power nonlinearity. Error estimates for the semidiscrete Lie-Trotter splitting method were established for the equation with both types of nonlinearities. Finally we investigated the dynamics and revealed the dispersion property of the mentioned two types of equations with different parameters by employing the proposed methods.

Figure 6 . 1 .

 61 Figure 6.1. Temporal errors of the Lie-Trotter splitting (3.8) and (3.10) for the solitary wave (6.1).

Figure 6 .

 6 Figure 6.1 displays the temporal errors of the Lie-Trotter splitting method (3.8) & (3.10) for the solitary wave. The left plot displays the errors at t N = T = 2.5 under different values of N .It can be clearly seen that via both cases of temporal grid, the Lie-Trotter splitting method converges at the first order in time. The right plot shows the evolution of the error when the step number N is fixed as N = 25000. We observe that the error increases exponentially in time for both cases and the error of Lie I increases a bit slower than that of Lie II, which suggests to choose Lie I for long time simulation. Figure6.2 shows the error e(2.5) L 2 for the methods Strang I and Strang II. It can be clearly observed that the splitting scheme converges at the second order in time and the error of Strang I is far less than that of Strang II. Moreover, for the method Strang II, it leads to a correct solution and is convergent quadratically only when the step number N is large enough N N T . Figure6.3 depicts the evolution of the error obtained via Strang I and Strang II with fixed step number N = 25000.

Figure 6 . 2 .

 62 Figure 6.2. Temporal errors of the Strang splitting scheme (3.9)-(3.10) for the solitary wave solution (6.1) under different step numbers.

L 2 )Figure 6 . 3 .

 263 Figure 6.3. Error evolution by using Strang I (left) and Strang II (right) for the solitary wave (6.1).

Figure 6 . 4 .

 64 Figure 6.4. Dynamics of u(t, x) (left) and κ(s, y) (right) for the solitary wave (6.1).

Figure 6 . 5 .

 65 Figure 6.5. Phase portraits for (2.4) with ω = 2 and λ = -3.

Figure 6 . 6 .

 66 Figure 6.6. Dynamics of u for Case (i) initial data in Example 2.

Figure 6 . 7 .

 67 Figure 6.7. Dynamics of κ (top) and u (bottom) for Case (ii) initial data in Example 2.

Figure 6 . 8 .Figure 6 . 9 .

 6869 Figure 6.8. Dynamics of κ(s, y) (top) and u(t, x) (bottom) for Case (iii) in Example 2.

Figure 6 . 10 .Figure 6 . 11 .

 610611 Figure 6.10. Phase portraits for the ODE (2.4) with ω = -λ = 2.

Figure 6 .

 6 Figure 6.11 shows the dynamics of κ until s = tanh(ωT ) ω

5 Figure 6 . 12 .

 5612 Figure 6.12. Dynamics of κ for Cases (i)-(ii) (up to down) in Example 4.

Figure 6 . 16 .

 616 Figure 6.16. Dynamics of v in the equation (4.5) with σ = 3 in Example 6.

Example 5. We set λ = 1, ω = 2 in the equation (4.1) with different σ and initial data u 0 : (i). u 0 (x) = 2e -x 2 ; (ii). u 0 (x) = sech(x 2 /2) sin(x).

-5 0 5 Re(v(s = tanh(ωt) ω , y)) recalled in (4.2), as well as the linear behavior, which involves, asymptotically as t → +∞, the same rescaling and phase shift as in (3.1).

Example 6. We set λ = -1, ω = 2 in the equation (4.1) with initial data u 0 (x) = 2e -x 2 and different σ.


Jean-Yves Delenne
email: jean-yves.delenne@inrae.fr

Luisa Orozco

Philippe Sornay

Farhang Radjai

DEM Modeling of particle breakage inside rotating drums

Keywords: Discrete Element Method, granular media, particle breakage, rotating drum, upscaling

In this work, we used the Contact Dynamics -Discrete Element Method (CD-DEM) to investigate the effect of material and system parameters on the grinding process in a 2D rotating drum. To model breakable particles we implemented the Bonded Cell Method (BCM) [1], in which the particles are discretized into bonded polygonal cells (figure 1a). A debonding criterion consistent with the classical framework of fracture mechanics (both in terms of yield stress and fracture energy) was employed [2]. We used a smooth drum without grinding media in which particle breakage is a consequence of granular flow. In this self-grinding or autogenous process (figure 1b), each particle breaks into fragments composed of unbreakable primary cells with different shapes and sizes depending on grinding time, surface energy, rotation speed and other mechanical properties. For an extensive parametric study [3] we performed long lasting simulations with large number of particles and cells in order to get meaningful statistics of fracture events. We varied system parameters such as drum size, rotation speed, filling degree and initial particle shape. The effect of each parameter on the granular flow and evolution of grinding in terms of the mean particle size and specific surface of the material was quantified. We show that the specific surface (defined as the sum of the surface areas of all particles divided by their total weight) increases almost linearly with time up to a transition point to a nonlinear regime where many unbreakable fragments are generated, and thus the probabil-ity of breakage declines. For all values of system parameters, this point corresponds to the same amount of specific surface equal to slightly more than half the maximum specific surface that can be generated in the simulations. This transition was used to define a characteristic time associated to the grinding efficiency. For all system parameters, when the times are scaled by this characteristic time, the rate of particle breakage collapses on a master curve. Finally, we show that the rate of particle breakage can be expressed as a linear function of a general scaling parameter that incorporates all system parameters. This scaling behavior provides a framework for the upscaling of drum grinding process from laboratory to industrial scale.

Figure 1 .

 1 Figure 1. a) Voronoï tessellation applied to polygonal breakable particles. Each cell is presented in a different color. b) Snapshots of a rotating drum simulation. The color is proportional to damage, defined as the number of cells detached from a particle, represented on color scale from bright green for intact particles to black for highly-damaged particles.

  ________________[1] D.-H. Nguyen, E. Azéma, F. Radjai, and P. Sornay, Physical Review E, 2014, 90, 012202.[2] L.F. Orozco, J.Y. Delenne, P. Sornay, F. Radjai, Int. J. SolidsStruct, 2019, 166, 47-56. [3] L. Orozco, D. Nguyen, J.-Y. Delenne, P.Sornay, and F. Radjai, Powder Technology, 2019, 362, p157.

17th European Symposium on Comminution & Classification (ESCC 2022) 27-29 June 2022, Toulouse, France


Christophe Hédon

Frédéric Schnell

Philippe Sosner

Frédéric Chagué

Iris Schuster

Marc Julia

Alexandre Duparc

Jean-Michel Guy

Nicolas Molinari

Lionel Michaux

Cardiac screening before returning to elite sport after SARS-CoV-2 infection

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.


A Jacques
email: alicia.jacques@inrae.fr

D Duclos

C Danchin-Burge

M J Mercat

M Tixier-Boichard

G Restoux

IDI: an index to assess the usefulness of cryopreserved collections for the management of animal genetic diversity

Confronted with the erosion of genetic diversity in domestic animal populations, many countries such as France have created gene banks to limit this loss. However, it is crucial to evaluate the potential of these collections before using them for genetic diversity management. Thus, we propose a new Index of Diversity Impact, IDI, to evaluate and compare the diversity potential of collections. IDI is negative when the use of collections decreases the overall relatedness of the current population and positive when it increases it. IDI was used to evaluate the cryopreserved collections of the French National Cryobank of 14 breeds of three livestock species, pig, sheep and cattle. We obtained significantly negative IDIs for three breeds that differ in terms of collection strategy and population management. The IDI appears to be a relevant tool to help managers and breeders to quickly identify promising collections to be used in future breeding programmes.

Introduction

Genetic diversity is essential for the sustainability of livestock populations by allowing for continuous response to selection and for adaptation to new conditions. However, selection has led to an erosion of genetic diversity over time, which can limit further genetic gain. In addition, for small populations, such as local breeds, the effect of genetic drift is even more important and leads to accelerated depletion of genetic diversity within the breed. To address this risk, many countries have invested in conserving genetic diversity through the creation of gene banks. These ex situ conservation entities have been storing genetic material for several years, including reproductive material such as semen. However, these genetic resources are still underused, even though they could have a potential interest for managing breed diversity. In their study, Leroy et al. (2011) showed that the use of frozen semen from cryobanked bulls can be a valid choice for introducing genetic diversity, provided that attention is paid to gene diffusion, or when there is a change in breeding goals. Two other studies showed that the use of cryopreserved semen from ancient bulls could introduce genetic diversity into a local breed under selection (Eynard et al., 2018) but also in an international breed [START_REF] Doekes | [END_REF]. France is a country with a rich livestock biodiversity, with for instance 47, 57 and 12 French breeds officially recognized in 2015 for cattle, sheep and pigs, respectively, but a majority of local ones are threatened [START_REF] Verrier | [END_REF]. French breeds greatly differ in terms of population size and management programme, with, for example, the Prim'Holstein, first French dairy breed under selection, having nearly 3,400,000 females born in the period 2016-2019 versus the local Tarentaise dairy breed, having only 17,266 females in the same period [START_REF] Danchin-Burge | Indicateurs de variabilité génétique -races bovines -Edition[END_REF]. This raises the question of genetic erosion and motivated the creation of the French National Cryobank, FNC, in 1999, in order to preserve the reproductive material from French breeds of domestic animals. Currently, the FNC brings together reproductive resources of 21 species and nearly 8,000 donors with stored material over the last 20 years. However, until today, the use of these available resources to manage the genetic diversity of current populations is very scarce. Thus, this study aims at characterizing the collections of the FNC and evaluating their potential use for the management of breed diversity in future breeding programmes.

Materials & Methods

Data for individuals in the National Cryobank. Donor animal data were obtained from the FNC database (www.cryobanque.org). We focused on three species of livestock that are likely to use frozen semen in their management: cattle, sheep and pig. The national identifier (unique identification in a national pedigree database), the breed and the year of birth of donor animals were recorded.

Breeds, pedigree data and definition of cohorts.

A total of 14 breeds were chosen to cover a large range of situations encountered in France, encompassing two pig breeds, the Pietrain (PIE) and the Cul Noir du Limousin (CNL), five sheep breeds, the Dairy Lacaune (LAC), the Manech Tête Noire (MTN), the Manech Tête Rousse (MTR), the Basco-Béarnaise (BAS) and the Mouton Vendéen (VEN), and seven cattle breeds, the Prim'Holstein (HOL), the Montbéliarde (MON), the Abondance (ABO), the Froment du Léon (FRO), the Tarentaise (TAR), the Saosnoise (SAO) and the Maraîchine (MAR). We defined contemporary male and female cohorts, M and F respectively, as individuals with at least one known parent, born between 2018 and 2020 for pigs, between 2017 and 2020 for sheep, and between 2015 and 2020 for cattle, according to the species' generation interval. These contemporary cohorts are assumed to be currently alive and available for breeding. A cohort composed of all the cryopreserved sires present in the FNC for each breed, C, was also defined. For the pig breed (PIE), 13 boars with semen cryopreserved by the breeding companies were added to our study. The pedigrees of all individuals in the three cohorts were extracted from the national database and limited to the 25 last generations. The pedigree quality of the contemporary male and female cohorts were evaluated with the overall equivalent generations as the weighted average of the generation equivalents (EqG) for each year using PEDIG [START_REF] Boichard | Proc. of the 7 th WCGALP[END_REF].

Computation of kinship and definition of an Index of Diversity Impact.

The calculations of the average kinship between the three cohorts were carried out using the colleau module of PEDIG [START_REF] Colleau | [END_REF]. For each breed, we calculated the average kinship between all the sires present in the cryobank (cohort C) and the contemporary female population (cohort F) (𝝓𝝓 𝑪𝑪-𝑭𝑭 ������� ); and between the contemporary male and female cohorts (𝝓𝝓 𝑴𝑴-𝑭𝑭 ������� ). The average kinship within cohort C between cryobank sires (𝝓𝝓 𝑪𝑪 ���� ) was also calculated. Then, we defined an Index of Diversity Impact (IDI) for cryobank individuals. This indicator is an estimation of the expected impact of the use of cryobank individuals in random mating with contemporary females on the level of diversity of the next generation under the assumption that only the males in the cryobank would be utilized. It is calculated following the Equation (1).

𝑰𝑰𝑰𝑰𝑰𝑰 = ( 𝝓𝝓 𝑪𝑪-𝑭𝑭 ������� -𝝓𝝓 𝑴𝑴-𝑭𝑭 ������� ) 𝝓𝝓 𝑴𝑴-𝑭𝑭 ������� � (1)
It is expressed as a percentage of loss (negative) or gain (positive) in mean kinship relative to the current kinship of the population. The confidence intervals (CI) at 95 % were defined from the individual IDI of each cryobank individual and applying Equation (2).

𝑪𝑪𝑰𝑰 = � 𝑰𝑰𝑰𝑰𝑰𝑰 ����� ± 𝟏𝟏, 𝟗𝟗𝟗𝟗 𝒔𝒔𝒔𝒔(𝑰𝑰𝑰𝑰𝑰𝑰) √𝒏𝒏 � (2)
With sd(IDI) the standard deviation observed within the cryobank population of n donors. All statistical analyses were conducted using R, version 4.0.2 (R Core Team, 2019).

Results

Composition of the cryobank's collections

For the three species, a total of 84 breeds are represented in the FNC, with 21 breeds for cattle, 15 for pigs and 48 for sheep. This corresponds to 977 bulls, 451 boars and 1,177 rams with frozen semen. The average number of donors per breed is 46.5 for cattle (min=1, max=210, se=13.4), 30.1 for pig (min=10, max=76, se=4.4) and 24.5 for sheep (min=1, max=302, se=6.4). Donor birth years ranged from 1959 to 2013, 1979 to 2018 and 1974 to 2019 for cattle, pig and sheep respectively. This shows that material from ancient sires was retrieved by the FNC after 1999. The range of birth years of donors varied according to species and breed (Table 1).

Genetic diversity and IDI for the 14 breeds analysed with pedigree data

The number of donors in the cryobank for the 14 selected breeds ranged from 10 bulls in the SAO to 302 rams in the LAC. The size of the contemporary cohorts was very variable according to breed, with the number of females ranging from 534 cows for the FRO to 5,069,731 for the HOL.

Table 1. 4 -16.3 [-20.8;-11.9] 1 Number of individuals with at less one parent known in the contemporary cohorts 2 Number of Equivalents generations in the pedigree computed for the contemporary cohorts 3 Confidence intervals at 95 %, computed from the individual IDI of each cryobank individual The quality of the pedigrees ranged from 3.1 to 19.1 EqG on the female pathway and from 3.1 to 19.3 EqG on the male pathway (Table 1). The mean within-cryobank kinship values ranged from 0.8 % to 24.3 %. The average kinship values between cryobank sires and contemporary females (cryo-F) ranged from 3.1% to 24.9%, and the average kinship values between contemporary males and females (M-F) ranged from 1.9% to 36.1%. The IDI values varied between -31% for the CNL and 66.5% for the SAO. Three breeds exhibiting different situations have a significantly negative IDI: a local pig breed with quite ancient sires in the cryobank (CNL), a moderately selected sheep breed with a large number of males in the cryobank as compared to contemporary males (VEN), and a highly selected cattle breed (HOL) for which the cryobank favored donors with extreme breeding values and includes old material, thus, increasing the genetic distance with the current population. The 11 remaining breeds have a positive or null IDI.

Discussion

With more than 8000 donor animals, the FNC needs to know the potential value of its collections for the management of genetic diversity, all the more that ancient sires are represented in the collections and that living populations can suffer from selection and drift. For sheep, the collections include donors born fairly continuously over time since 1985. For pigs, local breeds have ancient donors whereas breeds under selection have quite recent donors. For cattle, few donors are born after the 2010s following the switch to genomic selection, which considerably shortened the reproductive carrier of males and decreased the quantity of material available for long term storage. For the breeds studied, the number of donors in the collection is not correlated with the average intra-breed relatedness, which rather depends on the quality of the pedigree of the individuals in collections and the breed management programmes. The IDI reflects the average impact of randomly mating sires in collections with contemporary females relative to expectations for the use of contemporary males only; and makes it possible to rank the usefulness of collections resulting from different sampling strategies. A positive or null IDI reveals a representative collection of the current population, either by the presence of the strongest contributors to the breed (12 AI bulls in collection out of the 13 used in the FRO), or by recent collections likely more related to current females (example of PIE or BAS). However, a positive index does not definitively mean that the use of cryobanked sires is not useful for genetic diversity management, but that on average there will be fewer individuals in the collection that are genetically distant from the current population. In these cases, large confidence intervals suggest that the use of some well-chosen cryobank sires could produce matings with a lower level of relatedness than with contemporary sires, even if the IDI is positive, as for FRO or MON.

In conclusion, the IDI summarizes the impact, on genetic diversity, of using cryopreserved individuals relative to the current level of genetic diversity of the population; it is thus a convenient tool for managers to rapidly screen collections. In our study, IDI values were obtained from pedigree data but could also be calculated from genomic data, especially to be applicable to local breeds without pedigrees. The next step is to analyse options for individual mating plans between cryobank sires and contemporary females.

Table of genetic diversity measures for the 14 breeds studied by pedigree data.

		Cryobank sires	Contemporary cohort size 1	No. of EqG 2	Average Kinship (in %)	Index of Diversity Impact (in %)
	Breed No.	Birth year range	M	F	M	F	C	C-F M-F	IDI	[CI] 3
	SAO 10 1991;2011 3084 3089	3.2 3.3	00.8	03.1 01.9	66.5	[25.5;107.5]
	FRO 12 1978;2006 536	534	3.1 3.1	10.1	12.8 11.3	13.5	[-08.2;35.2]
	TAR 75 1994;2006 25781 26030	6.5 6.6	10.6	09.3 08.6	08.5	[03.1;14.0]
	ABO 66 1995;2006 470	71698	6.7 6.8	15.1	12.9 12.0	07.2	[01.1;13.4]
	MON 142 1994;2013 3658 1201388 9.4 9.5	12.1	10.8 10.8	00.0	[-03.0;03.0]
	MAR 13 1989;2004 3195 3368	6.7 6.5	05.7	07.6 08.0 -04.7 [-28.5;19.2]
	HOL 190 1993;2013 10776 5069731 10.4 10.0 10.8	08.9 10.1 -12.6 [-14.7;-10.6]
	PIE	27 2016;2017 49073 48503	19.3 19.1 12.8	12.7 12.4	02.3	[00.4;04.1]
	CNL 32 1976;2013 1243 1367	15.5 15.5 24.3	24.9 36.1 -31.0 [-35.4;-26.6]
	MTN 26 2007;2016 1209 11977	5.9 4.5	11.3	06.0 04.1	44.1	[34.1;54.1]
	MTR 48 2008;2016 6546 104184 9.7 8.2	06.8	05.5 05.1	07.2	[-01.3;15.7]
	LAC 302 1980;2018 19206 708459 13.3 10.8 05.7	03.5 03.4	03.3	[00.3;06.2]
	BAS 32 2009;2016 2418 32547	8.4 6.1	10.1	06.8 06.6	03.2	[-04.9;11.3]
	VEN 69 1989;2012 306	33855	13.7 10.3 05.2	03.7 04.


M Tixier-Boichard
email: michele.tixier-boichard@inrae.fr

N Peynot

D Duclos

S Weigend

L Monteagudo

A Martinez Martinez

J V Delgado

R Gonzalez Prendes

R Crooijmans

G Restoux

Mapping genetic diversity in European gene banks: preliminary results on chickens for the validation of IMAGE001 array

Gene banks are a component of a national strategy for the preservation of genetic diversity. Gene bank managers need to have a global and comparable picture of the diversity in their collections in order to rationalize them. Facing a diversity of molecular tools is a difficulty. The IMAGE H2020 project aimed at developing a low cost 60k SNP array to facilitate the mapping of diversity in gene banks. The first test of this array was performed for chicken with 204 samples from 18 local breeds and nine experimental lines provided by Germany, France and Spain. The MAF across population was 0.34, showing that this tool is useful over a range of populations. The principal component analysis and the Neighbor-joining tree showed that local breeds did not cluster according to country and were generally homogenous. Comparison with on-farm populations remains to be done to assess the value of the gene bank collections.

Introduction

Genetic diversity is essential for the sustainability of animal breeding. Countries are investing in long term conservation of genetic diversity with the creation of gene banks. FAO promotes this strategy in the Global Plan of Action for Animal Genetic Resources. In Europe, the European Regional Focal Point is gathering national coordinators of animal genetic resources from all countries and includes a permanent working group on ex situ cryoconservation. Until now, sampling strategies remain national and may vary according to countries, resulting in a variety of situations in terms of species and breeds represented. The ERFP is currently promoting the set up of a European network of Animal Gene Banks (EUGENA) to increase collaboration between national gene banks in Europe. An economic analysis of gene banking was recently performed in the frame of the IMAGE H2020 project (Grant # 677353). Results showed that the cost of preservation of a given number of breeds could be decreased by 20% by a better coordination across European countries (de Oliveira Silva et al., 2018). Assessing the value of gene bank collections for breeding programmes require data on withinbreed diversity as well as between-breeds, since both levels are important for the management of farmed animal species. Whereas pedigree data have been used since a long time to characterize within-breed diversity, for in situ populations as well as ex situ collections, only molecular data make it possible to map the diversity at all levels within a species. Mapping diversity has made significant progress thanks to the development of genotyping arrays in all farmed species, whereas whole genome sequencing can provide a deeper analysis of the status of collections. However, sequencing costs are still relatively high for gene bank managers, who generally lack the capacity to analyse such data. Thus, genotyping arrays are still a relevant tool to provide a global description of the diversity present in gene bank collections and to make comparisons between gene banks as well as between gene bank collections and on farm populations. However, several SNP arrays for chickens coexist, either private (around 60K) or public (around 580K), at a rather high cost (from 90 to 150 € per animal), with a low % of markers shared between the arrays. This is complicated for gene bank managers who need a global and comparable picture of the diversity in their collections. To facilitate a common approach across gene banks, the IMAGE H2020 project has developed a multi-species 60k array, with the aim to provide 10k markers for each of six species, i.e. cattle, sheep, goat, horse, pig and chicken, at a price of 20 US$ per sample maximum for both the array and the genotyping, thus excluding the cost of DNA extraction. A subset of gene bank collections from different countries was genotyped to validate the IMAGE001 array. In this paper, we show how the first version of IMAGE001 can map diversity of chicken collections from gene banks in Germany, France and Spain.

Materials & Methods

Gene bank samples.

A total of 211 chicken DNA samples were provided for 27 breeds from the gene banks of three countries (Table 1). Depending on the breed, the DNA had been extracted from frozen semen or from blood. Quality and concentration were determined with spectrophotometry and at least 1.5 µg of DNA was provided to the IMAGE partner at Wageningen University who was coordinating the test. Sampling represented only local breeds for Germany and Spain, with 10 to 20 animals per breed, whereas experimental lines were also provided by France, with a lower number of animals (2-5) per population (Table 1).

SNP selection.

SNP selection was performed based on 1) overlap with existing arrays with a high allele frequency across populations; 2) location in genes affecting phenotypic traits; 3) mtDNA; 4) Ancestral SNPs as compared to wild species of jungle fowl; 5) MHC region; 6) random genes located into QTL regions. A total of 9,306 SNPs were initially selected to be tested for the chicken, covering 32 autosomes (9071 markers) and the sex chromosomes (235 markers).

Genotyping.

Genotyping was performed on an Affymetrix platform by Eurofins® biotech company.

Analysis.

Analysis was performed with plink software. The Minor Allele Frequency (MAF) was calculated per chromosome and per breed. Autosomal markers were used to determine the genetic structure with a Principal Component Analysis and a neighbour-joining tree.

Results

Genotyping

Altogether 7743 autosomal markers passed the calling rate and were included in the analysis. Seven samples did not provide satisfactory results, leading to a total of 204 individual genotypes. The mean MAF per chromosome was above 0.30 for most of the autosomes, exceptions being chr 16 (0.03), chr 25 (0.27), chr 27 (0.20) chr 31 (0.02) and chr 33 (0.04). The overall MAF across all populations was 0.34. At population level, the median MAF was above 0.25 except for three congenic lines (LB4, LB13, LB19), which were less variable.

Principal component analysis

The first four principal components (PC) explained only 16% of the total variance. PC1 explained 7%, populations were scattered across a gradient from Langshan breed to congenic lines, Westfaelische Totleger and Ostfriesische_Moewen breeds. PC2 isolated the Gallina_del_Sobrarbe from other breeds, PC3 isolated the Deutsche Lachshuhn and PC4 isolated the Andalusian.

Neighbor-joining tree

The NJ tree (Figure 1) featured 4 main branches: branch A included only one Spanish breed (Blue Andalusian) and one German breed (Augsburger); branch B included three German breeds (Rheinlander, Westfaelische Totleger and Ostfriesische Moewen) and the four congenic lines which are all of White Leghorn origin; branch C included three German breeds (Krüper, Deutsche Sperber and Bergische Schlotterkamm); branch D included the remaining German breeds, most animals of the Gallina del Sobrarbe, the French local breeds and the experimental lines of a non White Leghorn origin. The divergent lines such as Fat/Lean lines and DPF+/DPFclustered together.

Discussion

First of all, the high MAF value observed across all populations demonstrates the ability of the array to be used over a large range of populations. The total number of markers was quite reduced and another set of markers will be added to reach the target of 10k markers. Yet, a few chromosomes, (16,31,33) are likely to remain underrepresented in the array. The populations did not particularly cluster according to country, which shows the interest of sharing information across gene banks to assess their complementarity. This finding, however, may result from the unbalanced distribution of breeds according to country in the data set. There was only one branch made up of breeds from a single country (Germany) which is consistent with the fact that these breeds are typical Northwestern European breeds and quite different from Mediterranean breeds such as the White Leghorn. Enlarging the sampling to Northern Europe would make possible to check this hypothesis. More information is needed on the history of the two Spanish breeds to understand their separate clustering pattern in the NJ tree.

Besides, the history of breeds from different countries which are located within the same branch of the NJ tree also needs to be further explored. Each population appeared in a single branch of the NJ tree, with a few exceptions, such as the Gallina del Sobrarbe, for which most animals appeared mainly in branch D but a few were in branch B. Local chicken breeds managed by a breeders'club generally appear to be homogenous, unless the breed is distributed across farmers who do not exchange their animals, thus leading to genetic fragmentation. In conclusion, using a single SNP array for mapping diversity in gene bank collections is feasible and the next step is to include data from on-farm populations in order to map the diversity within and outside gene banks.

Figure 1 .

 1 Figure 1. Neighbor-Joining tree of the populations. (A, B, C, D are discussed in the text)

  De Oliveira Silva R., Ahmadi B.V, Hiemstra S.J., and Moran D. (2018) J. Anim. Breed. 136(1) : 63 -73. https://doi.org/10.1111/jbg.12368 Acknowledgments : this work was funded by H2020 IMAGE grant #677353.

Table 1 -

 1 List of chicken populations included into the genotyping test for IMAGE001 array

	Type of Population	Breeds	Samples Country
		Congenic line_LB13 (MHC)	2	France
		Congenic line LB21 (MHC)	2	France
		Congenic line LB4 (MHC)	4	France
		Congenic line_LB19 (MHC)	3	France
	Experimental lines	DPF-(low duration of fertility line)	4	France
		DPF+ (high duration of fertility line)	5	France
		DWNA (dwarf naked neck layer)	5	France
		Fat line(high abdominal fatness)	5	France
		Lean line Low abdominal fatness)	5	France
		High quality Bresse chicken (LB99)	5	France
		Barbezieux	3	France
		Gasconne	3	France
		Gauloise_Doree	12	France
		Blue_Andalusian	10	Spain
		Deutsches Reichshuhn	10	Germany
		Krüper	10	Germany
		Deutsches Lachshuhn	10	Germany
	Local breeds	Langshan Ostfriesische Moewen	10 10	Germany Germany
		Rheinlaender	10	Germany
		Sachsenhuhn	10	Germany
		Sundheimer	10	Germany
		Westfaelische Totleger	10	Germany
		Augsburger	11	Germany
		Bergische Schlotterkamm	11	Germany
		Deutsche Sperber	11	Germany
		Gallina_del_Sobrarbe	20	Spain


,130) Z PD3m Liquidity

∆𝐸𝐷𝐹 1 𝑍 -𝑠𝑐𝑜𝑟𝑒 = [𝐸(𝑟𝑜𝑎) + 𝑐𝑎𝑟]/𝜎(𝑟𝑜𝑎) 𝜎(𝑟𝑜𝑎) 𝑐𝑎𝑟 + 𝑟𝑜𝑎 ≤ 0 𝑟𝑜𝑎 ≤ -𝑐𝑎𝑟 𝑐𝑎𝑟 < 𝜋 𝜋 (obs=18

https://www.newyorkfed.org/research/policy/rstar


A Guintard
email: ana.guintard@inrae.fr

L Bourguignon

G Restoux

V Ducrocq

Accounting for heterogeneity of data in implementing genomic selection models applicable to crossbred cattle in India

Indian dairy production takes place in very diverse environments. The objective of this study was to identify genotype by environment interactions in the context of establishment of a reference population for genomic selection. Data from 9,265 cows in six states were analysed. The high variability in the distribution of breed types and milk yield across India, as well as the wide range of climatic zones existing in the country are highlighted. This heterogeneity must be taken into account, defining, for example, groups of districts (administrative regions within states), rather than considering India as one homogeneous country, in order to better predict the genomic merit of individuals in their local environment.

Introduction

The dairy sector plays a major role in India, contributing 5.3% of the gross domestic product, and representing an essential source of animal protein in a country where one-third of the population is vegetarian (Belhekar and Dash, 2016). Improving productivity at the animal level is necessary to meet the demand and to ensure self-sufficiency [START_REF] Chawla | Milk and Dairy Products in India -Production, Consumption and Exports[END_REF]. Genomic selection is a way to increase the milk production in Indian smallholder farming systems (Al Kalaldeh et al., 2021). To implement it, BAIF Development Research Foundation (https://baif.org.in/) has been collecting performances and genotypes for the past four years to set up a reference population. However, the heterogeneity of collected data reflects a huge diversity of production systems over India (across states, climatic zones, herd size, breeds, etc.). Thus, the aim of this study was to analyse phenotypic data from BAIF's database in order to characterise the heterogeneity of performances to account for in the establishment of a reference population for genomic selection.

Material and Methods

Data used for this study originates from the BAIF database and encompass phenotypic data recorded in the states of Andhra Pradesh, Bihar, Gujarat, Jharkhand, Maharashtra, Odisha, Punjab, Rajasthan and Uttar Pradesh between August 2016 and November 2021. Animals were crossbred cattle ranging from pure indigenous Bos indicus to pure exotic Bos taurus (Holstein or Jersey) cows. The database, contains information on daily milk yield, day in milk, lactation number, last calving date, breed type, the type of recording (by an enumerator or by the farmer himself) the farm, the CDC (Cattle Development Centre, which is a local centre providing services such as artificial insemination for eight to ten villages), the district and the state.

In order to make the dataset suitable for analysis, we filtered the data as follows: (1) Andhra Pradesh and Gujarat states were removed because data were available only on pure zebus (breeding tracks), and were not representative of the diversity of breeds in these states, (2) data from cows with fewer than four records per lactation and with a lactation lasting less than 200 days were removed, (3) data were retained only if they were recorded between eight and 340 days after calving, (4) data from cows with a coefficient of variation in daily milk yield greater than 0.8 on a given lactation or a daily milk yield greater than 25kg were removed, being more than three times greater than the average daily milk yield (Al Kalaldeh et al., 2021), ( 5) data from cows with a lactation rank greater than seven were deleted, and (6) for each CDC and breed, data were deleted if they were recorded on fewer than ten cows.

Weather data (from worldweatheronline.com), temperature, humidity and precipitation were also considered at different levels of season which were defined as: winter (January to February), pre-monsoon (March to May), monsoon (June to September) and post-monsoon (October to December).

After data filtering and integration of weather data (not available for all the districts), 211,545 records from 9,265 cows in six states (Bihar, Jharkhand, Maharashtra, Odisha, Punjab and Uttar Pradesh), 26 districts, 75 CDC and 4,168 farms were retained for the analysis. The number of cows per herd ranged from one to 47 with an average of 2.2. Overall, 44% of the herds were composed of a single cow while 85% of the herds had three cows or less. The number of cows per CDC varied from ten to 456.

The daily milk yield was transformed into a 305-day production and corrected, using a linear mixed model including the effects of breed type, calving season, lactation rank and type of recording, all treated as fixed.

A recursive partitioning classification using decision trees allowing for clustering of individuals in homogenous sub-populations according to explanatory variables (Breiman et al., 1984), was performed on corrected milk production to identify regions and subsequent characteristics with contrasted production levels.

Statistical analyses were done using RStudio (version 3.5.0) and the ggplot2 package. Recursive partitioning classification was performed using the rpart package.

Results and discussion

Breed distribution.

To understand the variability of milk yield, we first looked at the distribution of uncorrected milk production according to breed types between and within states (Figure 1). The percentage of each breed between states varies a lot, ranging from places where indigenous zebus and nondescriptive cows are an important part of the herd (Bihar, Odisha and Uttar Pradesh), to states that orientated their breeding towards high levels of exotic breeds (Maharashtra, Jharkhand).

Holstein crosses usually show the highest milk production in most of the states. However, we can identify some scaling effect (a significantly higher production in Maharashtra than in other states) and re-ranking effect (in Odisha, Holstein crosses have, on average, a lower milk yield than indigenous breeds).

Figure 1 shows a high variability in the distribution of breed types across the six states. Corrected milk yield was significantly variable among districts (ANOVA, F=170.5, df =25, P<0.05)

Milk production under diverse climatic conditions.

The recursive portioning classification separated districts into four climatic clusters (figure 2), with the following predicted average corrected milk productions: 2,999 kg cluster A: Ahmednagar, sd= 649.4 kg, cv=0.2), 2,750 kg (cluster B: Faizabad, Hazaribag, Khunti, Lohardaga, Pune, Ranchi, Ropar and Saharanpur, sd=767.3 kg, cv=0.3), 2,459 kg (cluster C: Ara, Bahraich, Chhapra, Gonda, Hardoi, Hoshiarpur, Kannauj and Meerut, sd=759.4 kg, cv=0.4), and 2,237 kg (cluster D: Allahabad, Bargarh, Bolangir, Buxar, Etawah, Lucknow, Sambalpur, Unnao and Varanasi, sd=701.8 kg, cv=0.4).

Maharashtra shows the highest corrected milk production, especially the district of Ahmednagar (cluster A). This district, situated in a rain shadow region, has a tough climate, with very low rainfall, even during monsoon. The higher average milk production in Maharashtra, is more likely due to a longer technical breeding support of BAIF to local farmers in this state, where BAIF influence started more than 30 years ago. The higher milk production results from years of breeding with exotic breeds, and particularly Holstein cattle. This illustrates that genotype by environment interactions do not only depend on the state, and that beyond climatic conditions, other factors greatly impact production. Particular attention should be paid to the farming system, and further analyses of agro-economic data (major crops, size of the farm, family income, economic status and feeding system of the cows) should help to identify significant determinants of milk yield in India.

Figure 2 also shows that homogeneous groups of weather are not always located within the same state. Indeed, all dots of the same cluster share similar climatic conditions, and we can see that districts of the same cluster can be found in different locations. Mountains and rivers (especially the Ganges river) play an important role in the local weather, resulting in various sub-groups of climates, influencing crops and feedstuff in the area. Considering these regional variations by classifying districts into homogenous climatic zones could help to account for the wide heterogeneity of weathers across India in a genomic evaluation.

Conclusion

India, given its size and its particular milk production, mostly ensured by small and marginal farmers, shows extreme variability in terms of production, breeds and climates, resulting in strong genotype by environment interactions. It, thus, appears crucial to account for them, along with agro-economic characteristics, when establishing a reference population in order to get an appropriate estimation of the genomic breeding value of individuals. In particular, defining a few groups of homogeneous environments across India, rather than considering it as a uniform country, could help to conduct more accurate genomic evaluation in the field.

Figure 1 :

 1 Figure 1: Milk Yield production (MY 305) and percentage of each breed type within each state.