Christophe Besse

AND Rémi Carles
email: remi.carles@math.cnrs.fr

Florian Méhats
email: florian.mehats@univ-rennes1.fr

An asymptotic preserving

est destinée au dépôt

Introduction and main results

1.1. Motivation. We consider the cubic, defocusing, nonlinear Schrödinger equation (NLS)

iε∂ t u ε + ε 2 2 ∆u ε = |u ε | 2 u ε , (t, x) ∈ R + × R d , (1.1)
with WKB type initial data u ε (0, x) = a 0 (x)e iφ 0 (x)/ε , (

the functions a 0 and φ 0 being real-valued. The aim of this article is to construct an Asymptotic Preserving (AP) numerical scheme for this equation in the semiclassical limit ε → 0. We seek a scheme which provides an approximation of the solution u ε with an accuracy, for fixed numerical parameters ∆t, ∆x, that is not degraded as the scaled Planck constant ε goes to zero. In other terms, such a scheme is consistent with (1.1) for all fixed ε > 0, and when ε → 0 converges to a consistent approximation of the limit equation, which is the compressible, isentropic Euler equation (1.4). The difficulty here is that, when ε is small, the solution u ε becomes highly oscillatory with respect to the time and space variables, and converges to its limit only in a weak sense. In order to follow these oscillations without reducing ∆t and ∆x at the size of ε, which may be computationally demanding, our construction relies on a fluid reformulation of (1.1), well adapted to the semiclassical limit.

In kinetic theory, plasma physics and radiative transfer, the interest in AP schemes has considerably grown in the recent years since the pioneer works on the subject [START_REF] Coron | Numerical passage from kinetic to fluid equations[END_REF][START_REF] Gabetta | Relaxation schemes for nonlinear kinetic equations[END_REF][START_REF] Jin | Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations[END_REF][START_REF] Klar | An asymptotic-induced scheme for nonstationary transport equations in the diffusive limit[END_REF]. Among a long list of works on this topic, we can mention [START_REF] Jin | Uniformly accurate diffusive relaxation schemes for multiscale transport equations[END_REF][START_REF] Gosse | Asymptotic-preserving & well-balanced schemes for radiative transfer and the Rosseland approximation[END_REF][START_REF] Buet | Asymptotic preserving and positive schemes for radiation hydrodynamics[END_REF][START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF][START_REF] Bennoune | Uniformly stable numerical schemes for the Boltzmann equation preserving the compressible Navier-Stokes asymptotics[END_REF][START_REF] Crispel | An asymptotic preserving scheme for the twofluid Euler-Poisson model in the quasineutral limit[END_REF][START_REF] Filbet | A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources[END_REF][START_REF] Carrillo | Numerical schemes of diffusion asymptotics and moment closures for kinetic equations[END_REF][START_REF] Degond | Asymptoticpreserving particle-in-cell method for the Vlasov-Poisson system near quasineutrality[END_REF][START_REF] Lafitte | Asymptotic-preserving projective integration schemes for kinetic equations in the diffusion limit[END_REF][START_REF] Degond | An asymptotic-preserving method for highly anisotropic elliptic equations based on a micro-macro decomposition[END_REF]. For the Schrödinger equation, fewer works exist. In the stationary case, one can cite [START_REF] Arnold | WKB-based schemes for the oscillatory 1D Schrödinger equation in the semiclassical limit[END_REF] which proposes a WKB-type transformation in order to filter out the oscillations in space. In the time-dependent case, the usual numerical methods for solving the nonlinear Schrödinger equationfinite-difference schemes [START_REF] Wu | Dufort-Frankel-type methods for linear and nonlinear Schrödinger equations[END_REF][START_REF] Delfour | Finite-difference solutions of a nonlinear Schrödinger equation[END_REF][START_REF] Akrivis | Finite difference discretization of the cubic Schrödinger equation[END_REF][START_REF] Karakashian | On optimal order error estimates for the nonlinear Schrödinger equation[END_REF], time-splitting methods [START_REF] Pathria | Pseudo-spectral solution of nonlinear Schrödinger equations[END_REF][START_REF] Weideman | Split-step methods for the solution of the nonlinear Schrödinger equation[END_REF][START_REF] Besse | Order estimates in time of splitting methods for the nonlinear Schrödinger equation[END_REF], relaxation scheme [START_REF] Besse | A relaxation scheme for the nonlinear Schrödinger equation[END_REF] or even symplectic methods [START_REF] Sanz-Serna | Conservative and nonconservative schemes for the solution of the nonlinear Schrödinger equation[END_REF] -are designed in order to guarantee the convergence of the discrete wave function u ε when ε > 0 is fixed. As it is analyzed for instance in [START_REF] Markowich | Numerical approximation of quadratic observables of Schrödinger-type equations in the semi-classical limit[END_REF][START_REF] Markowich | A Wigner-measure analysis of the Dufort-Frankel scheme for the Schrödinger equation[END_REF] for the case of finite-differences or in [START_REF] Klein | Fourth order time-stepping for low dispersion Korteweg-de Vries and nonlinear Schrödinger equations[END_REF] for Runge-Kutta schemes, these methods suffer from the oscillations in the semiclassical regime if the time and space steps are not shrunk. Note however that, among these methods, the status of splitting methods is particular. Indeed, in the linear case, it is shown in [START_REF] Bao | On time-splitting spectral approximations for the schrödinger equation in the semiclassical regime[END_REF] (see also the recent review article [START_REF] Jin | Mathematical and computational methods for semiclassical Schrödinger equations[END_REF]) by Wigner techniques that, if the space step ∆x has to follow the parameter ε when it becomes small, the time step can be chosen independently of ε (before the appearance of caustics) if we are only interested in getting a good approximation of the quadratic observables associated to u ε . Still, none of these methods is Asymptotic Preserving in the sense defined above.

The quantum fluid reformulation of the Schrödinger equation provides a natural framework for the construction of AP numerical methods. The Madelung transform [START_REF] Madelung | Quanten theorie in hydrodynamischer form[END_REF] is a polar decomposition of the solution of (1.1), written as u ε (t, x) = a ε (t, x)e iφ ε (t,x)/ε , where the amplitude a ε and the phase φ ε are real-valued. Inserting this ansatz into (1.1) yields the system of quantum hydrodynamics (QHD), or Madelung equations, satisfied by ρ ε = (a ε ) 2 and v ε = ∇φ ε :

     ∂ t ρ ε + div (ρ ε v ε ) = 0, ρ ε |t=0 = (a 0 ) 2 , ∂ t v + v ε • ∇v ε + ∇ρ ε = ε 2 2 ∇ ∆ √ ρ ε √ ρ ε , v ε |t=0 = ∇φ 0 . (1.3)
This system takes the form of the compressible Euler equation (with the pressure law p(ρ) = ρ 2 /2), with an additional term testifying for the quantum character of the equation, the so-called Bohm potential ε 2 2 ∇ ∆ √ ρ ε √ ρ ε . Two comments are in order. First, the term in ε does not appear any more as a singular perturbation and, in the limit ε → 0, the quantum pressure disappears, leading formally to the compressible, isentropic Euler equation

∂ t ρ + div (ρv) = 0, ρ |t=0 = (a 0 ) 2 , ∂ t v + v • ∇v + ∇ρ = 0, v |t=0 = ∇φ 0 , (1.4)
which is the semiclassical limit of our problem (actually valid before the appearance of shocks). Second, the fluid-like form of this system enables to consider many numerical methods originated from computational fluid mechanics.

For the linear Schrödinger equation, the Madelung formulation is at the heart of the method of quantum trajectories (or Bohmian dynamics), see [START_REF] Wyatt | Quantum Dynamics with Trajectories: Introduction to Quantum Hydrodynamics[END_REF], which are based on a Lagrangian resolution of this system. In Eulerian coordinates, the work [START_REF] Degond | An asymptotic preserving scheme for the Schrödinger equation in the semiclassical limit[END_REF] also exploits this formulation to construct Asymptotic-Preserving schemes for (1.3). Unfortunately, the main drawback of the QHD system is the form of the quantum potential, which becomes singular as the density ρ ε vanishes (see [START_REF] Carles | [END_REF] for a recent survey). Hence, these methods do not provide AP schemes in the presence of vacuum. Another point of view was recently developed in [START_REF] Carles | Numerical aspects of the nonlinear Schrödinger equation in the semiclassical limit in a supercritical regime[END_REF] for the nonlinear equation (1.1), taking advantage of another fluid reformulation of the equation -due to Grenier [START_REF] Grenier | Semiclassical limit of the nonlinear Schrödinger equation in small time[END_REF] -which tolerates the presence of vacuum. In the next Subsection 1.2, we present in detail this formulation relying on a complex amplitude version of Madelung transform. However, as we will see, the drawback of this reformulation is that this system may develop singularities even for fixed ε > 0, whereas the solution of NLS remains smooth. Then, in the last Subsection 1.3 of this introduction, we present a modification in order to remedy this problem and preserve the regularity of the solution. The new QHD model that we propose contains a viscosity term in the equation for the velocity, but is still equivalent to the original equation (1.1).

1.2. Backgrounds on NLS and its semiclassical limit. In this section, we recall some known results on the Cauchy problem for NLS (1.1) and on its semiclassical limit as ε → 0. For details on this subject, we refer to the textbooks [START_REF] Cazenave | Semilinear Schrödinger equations[END_REF] and [START_REF] Carles | Semi-classical analysis for nonlinear Schrödinger equations[END_REF].

The following result is standard in the case ε = 1, and can easily be deduced when ε ∈ (0, 1], by scaling arguments.

Proposition 1.1. Let d 3, and u ε 0 ∈ H 1 (R d ). (i) There exists a unique solution u ε ∈ C(R; H 1 (R d ))∩L 8/d loc (R; W 1,4 (R d )) to (1.1), such that u ε |t=0 = u ε 0 .
Moreover, the following conservations hold:

Mass: d dt u ε (t) 2 L 2 = 0. Momentum: d dt Im R d u ε (t, x)∇u ε (t, x)dx = 0.
Energy:

d dt ε∇u ε (t) 2 L 2 + u ε (t) 4 L 4 = 0. (ii) If in addition u ε 0 ∈ H k (R d ) for k ∈ N, k 2, then u ε ∈ C(R; H k (R d )). Note that, if u ε = a ε e iφ ε /ε with φ ε ∈ R, the above three conservation laws become d dt a ε (t) 2 L 2 = 0. (1.5) d dt Im R d a ε (t, x)∇a ε (t, x)dx + 1 ε R d |a ε (t, x)| 2 ∇φ ε (t, x)dx = 0. d dt ε∇a ε (t) + ia ε (t)∇φ ε (t) 2 L 2 + a ε (t) 4 L 4 = 0. (1.6)
Let us now describe the limit ε → 0 for (1.1). The toolbox for studying semiclassical Schrödinger equations contains a variety of methods, depending on the quantities we are interested in. If one is only interested in a description of the dynamics of quadratic observables, such as the mass, current or energy densities, the Wigner transform is well adapted and has been applied to linear Schrödinger equations and Schrödinger-Poisson systems [START_REF] Gérard | Homogenization limits and Wigner transforms[END_REF][START_REF] Lions | Sur les mesures de Wigner[END_REF][START_REF]Wigner measure and semiclassical limit of Schrödinger-Poisson equations[END_REF]; on the other hand, it is not adapted to the study of NLS [START_REF] Carles | On the time evolution of Wigner measures for Schrödinger equations[END_REF]. Still for a description of quadratic observables, the modulated energy method [START_REF] Brenier | Convergence of the Vlasov-Poisson system to the incompressible Euler equations[END_REF][START_REF] Zhang | Semiclassical limit of nonlinear Schrödinger equation[END_REF][START_REF] Lin | Semiclassical limit of the Gross-Pitaevskii equation in an exterior domain[END_REF][START_REF] Alazard | Loss of regularity for super-critical nonlinear Schrödinger equations[END_REF] enables to treat the nonlinear equation (1.1). Yet, for a pointwise description of the wavefunction u ε , WKB techniques are more convenient. We refer to [START_REF] Carles | Semi-classical analysis for nonlinear Schrödinger equations[END_REF] for a presentation of WKB analysis of Schrödinger equations.

In the case of (1.1), also referred to as supercritical geometric optics in this context, the justification of the WKB expansion was done by Grenier [START_REF] Grenier | Semiclassical limit of the nonlinear Schrödinger equation in small time[END_REF]. Let us briefly present his idea. As we said above, the main inconvenience of the QHD system (1.3) is the singularity of the Bohm potential at the zeroes of the function ρ ε (i.e. at the zeroes of the wavefunction u ε ). Looking again for solutions u ε under the form

u ε (t, x) = a ε (t, x)e iφ ε (t,x)/ε , (1.7)
but allowing now a ε to take complex values, Grenier proposed to define a ε and φ ε as the solutions of the system

     ∂ t a ε + ∇φ ε • ∇a ε + 1 2 a ε ∆φ ε = i ε 2 ∆a ε , a ε |t=0 = a 0 , ∂ t φ ε + 1 2 |∇φ ε | 2 + |a ε | 2 = 0, φ ε |t=0 = φ 0 .
(1.8)

This system, which is formally equivalent to (1.1), can also be expressed in terms of the unknowns a ε and v ε = ∇φ ε . Differentiating with respect to x the second equation in (1.8), we find, using the fact that v ε is irrotational:

   ∂ t a ε + v ε • ∇a ε + 1 2 a ε div v ε = i ε 2 ∆a ε , a ε |t=0 = a 0 , ∂ t v ε + v ε • ∇v ε + ∇|a ε | 2 = 0, v ε |t=0 = ∇φ 0 .
(1.9)

The important remark made in [START_REF] Grenier | Semiclassical limit of the nonlinear Schrödinger equation in small time[END_REF] is to notice that the above system is hyperbolic symmetric, perturbed by a skew-symmetric term. In the limit case ε = 0, and if a 0 is real-valued and nonnegative, (1.9) is the Euler equation (1.4) written in symmetric form (see e.g. [START_REF] Makino | Sur la solution à support compact de l'équation d'Euler compressible[END_REF][START_REF] Chemin | Dynamique des gaz à masse totale finie[END_REF]):

   ∂ t a + v • ∇a + 1 2 a div v = 0, a |t=0 = a 0 , ∂ t v + v • ∇v + ∇|a| 2 = 0, v |t=0 = ∇φ 0 . (1.10)
This remark enables to obtain, in small time (but independent of ε ∈ [0, 1]), uniform estimates of the solution (a ε , v ε ) of (1.9) in Sobolev norms and prove the convergence of these quantities to the solution (a, v) of (1.10), as long as this solution exists and is smooth. From the numerical point of view, the regular structure of (1.9), in particular when a ε vanishes or has very small values, is an advantage compared to (1.3), and is at the basis of the schemes constructed in [START_REF] Carles | Numerical aspects of the nonlinear Schrödinger equation in the semiclassical limit in a supercritical regime[END_REF].

1.3. New model and main results. A drawback of (1.9) is that large time existence results (for fixed ε > 0) do not seem available. In particular, it is not known whether (1.9) still has a smooth solution past the time where the solution to (1.10) has ceased to be smooth (such a time necessarily exists if φ 0 and a 0 are compactly supported, regardless of their size, from [START_REF] Makino | Sur la solution à support compact de l'équation d'Euler compressible[END_REF] -see also [START_REF] Xin | Blowup of smooth solutions of the compressible Navier-Stokes equation with compact density[END_REF]). A practical consequence of this fact is that a numerical method based on (1.9), such as the one proposed in [START_REF] Carles | Numerical aspects of the nonlinear Schrödinger equation in the semiclassical limit in a supercritical regime[END_REF], may not approximate correctly the original equation (1.1) on arbitrary time intervals, for fixed ε > 0.

To overcome this issue, we use the degree of freedom given by (1.7) in a manner which is slightly different from the approach introduced by E. Grenier. From (1.7), we have

iε∂ t u ε + ε 2 2 ∆u ε -|u ε | 2 u ε = -∂ t φ ε + 1 2 |∇φ ε | 2 + |a ε | 2 a ε e iφ ε /ε + iε ∂ t a ε + ∇φ ε • ∇a ε + 1 2 a ε ∆φ ε -i ε 2 ∆a ε e iφ ε /ε .
Introduce the viscous term ε 2 a ε e iφ ε /ε ∆φ ε , and reorder the terms so as to consider the new system

     ∂ t a ε + ∇φ ε • ∇a ε + 1 2 a ε ∆φ ε = i ε 2 ∆a ε -iεa ε ∆φ ε , a ε |t=0 = a 0 , ∂ t φ ε + 1 2 |∇φ ε | 2 + |a ε | 2 = ε 2 ∆φ ε , φ ε |t=0 = φ 0 .
(1.11)

Proceeding like above, we get the following new system:

   ∂ t a ε + v ε • ∇a ε + 1 2 a ε div v ε = i ε 2 ∆a ε -iεa ε div v ε , a ε |t=0 = a ε 0 , ∂ t v ε + v ε • ∇v ε + 2 Re (a ε ∇a ε ) = ε 2 ∆v ε , v ε |t=0 = ∇φ 0 .
(1.12)

We recover physical observables, respectively particle, current and energy densities, thanks to the following formula

ρ ε = |a ε | 2 j ε = ε Im(a ε ∇a ε ) + ρ ε v ε , e ε = |ε∇a ε + ia ε v ε | 2 + |a ε | 4 .
(1.13) Remark 1.2 (Linear case). In the case of a linear Schrödinger equation

iε∂ t u ε + ε 2 2 ∆u ε = V u ε ,
we will see that the introduction of this artificial diffusion is even more striking than in the nonlinear framework; see Subsection 4.1.

Remark 1.3 (Irreversibility). Our choice for introducing the viscous term in (1.11) makes the system irreversible (in time), while (1.1) is reversible. To solve the Schrödinger equation backward in time, it suffices to change the signs to

     ∂ t a ε + ∇φ ε • ∇a ε + 1 2 a ε ∆φ ε = i ε 2 ∆a ε + iεa ε ∆φ ε , a ε |t=0 = a 0 , ∂ t φ ε + 1 2 |∇φ ε | 2 + |a ε | 2 = -ε 2 ∆φ ε , φ ε |t=0 = φ 0 .
Let us now present our main theoretical results on this system. Our first result concerns the local in time well-posedness of the Cauchy problem in any dimension, for fixed values of ε 0.

Theorem 1.4. Let (a 0 , φ 0 ) ∈ H s (R d ) × H s+1 (R d ), with s > d/2 + 1.
Then the following holds. (i) The system (1.12) admits a unique maximal solution

(a ε , v ε ) ∈ C([0, T ε max ); H s × H s ). (ii) There exists T > 0 independent of ε ∈ [0, 1] such that T ε max T . Moreover, the L ∞ ([0, T ]; H s × H s ) norm of (a ε , v ε ) is bounded uniformly in ε ∈ [0, 1]. (iii) Defining φ ε ∈ C([0, T ε max ); H s+1 ) by φ ε (t, x) = φ 0 (x) - t 0 1 2 |v ε (s, x)| 2 + |a ε (s, x)| 2 -ε 2 div v ε (s, x) ds, (1.14)
then the function

u ε = a ε e iφ ε /ε ∈ C([0, T ε max ); H s ) is the unique solution to (1.1) satisfying (1.2).
Remark 1.5. In fact, the proof of this theorem provides a criteria of global existence for the solution to (1.12). Indeed, we shall prove that the lifespan T ε max is independent of s > d/2 + 1 and that

T ε max < ∞ =⇒ T ε max 0 (a ε , v ε )(t) W 1,∞ + a ε (t) 2 L ∞ dt = ∞. (1.15)
Our second result states the convergence of the solution to (1.12) towards the solution of the Euler equation, as ε → 0.

Notation. For two positive numbers α ε and β ε , the notation α ε β ε means that there exists C > 0 independent of ε such that for all ε ∈ (0, 1],

α ε Cβ ε . Theorem 1.6. Let s > d/2 + 3, and (a 0 , φ 0 ) ∈ H s × H s+1 . Let T > 0 such that the Euler equation (1.4) has a unique solution (ρ, v) ∈ C([0, T ]; H s × H s ). Then (1.10) has a unique solution (a, v) ∈ C([0, T ]; H s × H s ). Moreover, for ε > 0 sufficiently small, we have T ε max T and (a ε , v ε ) -(a, v) L ∞ ([0,T ];H s-2 ×H s-2 ) ε.
Theorems 1.4 and 1.6 are not different from the corresponding results obtained by Grenier in [START_REF] Grenier | Semiclassical limit of the nonlinear Schrödinger equation in small time[END_REF] on his system (1.9). We now state a global existence result for fixed ε > 0 in dimension one, that is not available on (1.9). This result takes advantage of the viscous term in the second equation of (1.12).

Theorem 1.7. Suppose d = 1 and let (a 0 , φ 0 ) ∈ H s × H s+1 , with s 2. Then for ε > 0 fixed, the solution to (1.12) is global in time, in the sense that

T ε max = ∞ in Theorem 1.4.
This article is organized as follows. Section 2 is devoted to the proofs of our three theorems: the local existence result in Subsection 2.1, the semiclassical limit in Subsection 2.2 and the global existence result in Subsection 2.3. Section 3 is of different nature and concerns numerics. We describe a second-order Asymptotic Preserving numerical scheme for NLS, based on the reformulation (1.12), and we give the results of numerical experiments, in dimensions 1 and 2. We show that our scheme is indeed AP before the formation of singularities in the Euler equation, which means that the time and space steps ∆t and ∆x can be taken independently of ε. After the formation of singularities, we show that, in order to get a good approximation of the solution, we need to take ∆t and ∆x of order O(ε). At the end of the article, in Section 4, we sketch two extensions of our results: we discuss on the case of the linear Schrödinger equation, then we discuss on other nonlinearities.

Proofs of the main theorems

This section is devoted to the proofs of our main results stated in Theorems 1.4, 1.6 and 1.7.

Local existence.

In this paragraph, we show the local existence of a smooth maximal solution to (1.12).

Proof of Theorem 1.4. Items (i) and (ii): local existence and uniform estimates. Following [START_REF] Grenier | Semiclassical limit of the nonlinear Schrödinger equation in small time[END_REF], introduce the vector-valued unknown function

U ε =        Re a ε Im a ε v ε 1 . . . v ε d        =        a ε 1 a ε 2 v ε 1 . . . v ε d        ∈ R d+2 .
In terms of this unknown function, (1.12) reads

∂ t U ε + d j=1 A j (U ε )∂ j U ε = ε 2 LU ε + ε 2 DU ε + ε d j=1 B j (U ε )∂ j U ε , (2.1)
where the (k, ) 1 k, d+2 elements of the matrices A j (U ) ∈ M d+2 (R) are given by

             A j k,k (U ) = U j+2 for k = 1, . . . , d + 2 A j 1,j+2 (U ) = U 1 /2, A j 2,j+2 (U ) = U 2 /2 A j j+2,1 (U ) = 2U 1 , A j j+2,2 (U ) = 2U 2 A j k, (U ) = 0 otherwise. The linear operator L is given by L =   0 -∆ 0 . . . 0 ∆ 0 0 . . . 0 0 0 0 d×d   .
A first important remark is that even though L is a differential operator of order two, it causes no loss of regularity in the energy estimates, since it is skew-symmetric. Also, (2.1) is hyperbolic symmetric (or symmetrizable). Indeed, let

S = I 2 0 0 1 4 I d . (2.2)
This matrix is symmetric and positive, and SA j is symmetric for all k

SA j (U ) ∈ S d+2 (R), ∀U ∈ R d+2 .
The diffusive term D is given by

D =   0 0 0 . . . 0 0 0 0 . . . 0 0 0 I d×d ∆   .
Finally, the matrices B j , defined analogously to A j , are given by

B j 1,j+2 (U ) = U 2 , B j 2,j+2 (U ) = -U 1 B j
k, (U ) = 0 otherwise. When the right-hand side of (2.1) is zero, local existence of a unique solution in H s with s > d/2 + 1 is a consequence of standard quasilinear analysis for hyperbolic symmetric systems; see e.g. [START_REF] Taylor | Partial differential equations[END_REF]. The important point in energy estimates consists in computing

d dt Λ s (SU ε ), Λ s U ε , Λ = (1 -∆) 1/2 ,
and taking advantage of symmetry to obtain, thanks to tame estimates,

d dt Λ s (SU ε ), Λ s U ε C U ε W 1,∞ Λ s (SU ε ), Λ s U ε
, for some constant C depending only on A and d.

In the case of the complete system (2.1), we first recall that as noticed in [START_REF] Grenier | Semiclassical limit of the nonlinear Schrödinger equation in small time[END_REF], the term L does not alter the above energy estimate, since SL is skew-symmetric, so we have

d dt Λ s (SU ε ), Λ s U ε C U ε W 1,∞ Λ s (SU ε ), Λ s U ε + 2ε 2 SΛ s DU ε , Λ s U ε + 2ε d j=1 SΛ s B j (U ε )∂ j U ε , Λ s U ε .
In terms of (v ε , a ε ), the first new term reads

2ε 2 SΛ s DU ε , Λ s U ε = ε 2 2 d k=1 ∆ (Λ s v ε k ) Λ s v ε k = - ε 2 2 |∇Λ s v ε | 2 .
For the second new term, we have

2ε SΛ s B j (U ε )∂ j U ε , Λ s U ε = ε 2 Λ s a ε 2 ∂ j v ε j Λ s a ε 1 - ε 2 Λ s a ε 1 ∂ j v ε j Λ s a ε 2 .
By Kato-Ponce commutator estimates [START_REF] Kato | Commutator estimates and the Euler and Navier-Stokes equations[END_REF],

Λ s a ε k ∂ j v ε j -a ε k Λ s ∂ j v ε j L 2 ∇a ε k L ∞ v ε j H s + a ε k H s ∇v ε L ∞ .
Then using Cauchy-Schwarz and Young inequalities,

a ε k Λ s ∂ j v ε j Λ s a ε ∇Λ s v ε L 2 a ε k Λ s a ε L 2 δε ∇Λ s v ε 2 L 2 + 1 4δε a ε k Λ s a ε 2 L 2 .
Choosing δ > 0 sufficiently small and independent of ε, we infer

d dt Λ s (SU ε ), Λ s U ε + ε 2 4 ∇Λ s v ε 2 L 2 U ε W 1,∞ U ε 2 H s + a ε 2 L ∞ a ε 2 H s . (2.
3) The local existence result then follows easily by adapting the standard arguments from [START_REF] Taylor | Partial differential equations[END_REF]. In view of the energy estimate (2.3), we also infer (1.15). The lifespan T ε max of the maximal solution to (1.12) must be expected to actually depend on ε. For instance, if φ 0 , a 0 ∈ C ∞ 0 (R d ), then from [START_REF] Makino | Sur la solution à support compact de l'équation d'Euler compressible[END_REF], T 0 max < ∞, regardless of the dimension d. On the other hand, we prove further (proof of Theorem 1.7

) that if d = 1, T ε max = ∞ for all ε > 0. Item (iii): u ε is solution of the original NLS problem. Define φ ε by (1.14). Let us prove that (a ε , φ ε ) ∈ C([0, T ε max ); H s × H s+1 ) solves the system (1.11). Since (a ε , v ε ) ∈ C([0, T ε max ); H s ×H s ), we readily have φ ε ∈ C([0, T ε max ); H s-1
). Moreover, v ε is irrotational: indeed, curl v ε satisfies a homogeneous equation (see e.g. [4, p. 291]), it is zero at time t = 0, so curl v ε ≡ 0. Set Ω ε = Dv ε -∇v ε , where Dv ε stands for the Jacobian matrix of v ε , and ∇v ε stands for its transposed matrix. It solves

∂ t Ω ε + v ε • ∇Ω ε + Ω ε • Dv ε + ∇v ε • Ω ε = ε 2 ∆Ω ε . As a consequence, ∇|v ε | 2 = 2v ε • ∇v ε .
We then deduce from (1.14) and from the second equation of (1.12)

∂ t (∇φ ε -v ε ) = ∇∂ t φ ε -∂ t v ε = 0, so we infer from the initial condition v ε |t=0 = ∇φ 0 = ∇φ ε |t=0 that v ε = ∇φ ε . There- fore, ∇φ ε ∈ C([0, T ε max ); H s ), hence φ ε ∈ C([0, T ε max ); H s+1 ).
Replacing v ε with ∇φ ε in (1.14) shows that φ ε solves the first equation in (1.11), and the second equation in (1.11) is obviously satisfied. Now, define u ε by (1.7). The property

u ε ∈ C([0, T ε ]; H s ) is straightforward, since H s (R d ) is an algebra. It is clear that u ε solves (1.1) and (1.2). Since C([0, T ε ]; H s (R d )) ⊂ C([0, T ε ]; H 1 (R d )) ∩ L 8/d ([0, T ε ]; W 1,4 (R d )),
uniqueness stems from Proposition 1.1.

Convergence towards the Euler equation.

In this subsection, we study the semiclassical limit of (1.12).

Proof of Theorem 1.6. The proof proceeds as for [START_REF] Grenier | Semiclassical limit of the nonlinear Schrödinger equation in small time[END_REF]Theorem 1.2]. First, Theorem 1.4 yields a solution (a, v) ∈ C([0, T 0 max ); H s × H s ). Necessarily, T 0 max > T (where T is an existence time for the Euler equation), for if we had T 0 max T , then by uniqueness for (1.4), (|a|

2 , v) = (ρ, v) ∈ L ∞ ([0, T ]; H s × H s ). The equation for a in (1.10) then shows that a ∈ L ∞ ([0, T ]; H s-1 ), hence (a, v) ∈ L ∞ ([0, T ]; H s-1 × H s ), which is impossible if T 0 max
T : this yields the first part of the theorem. To prove the error estimate, introduce the vector-valued unknown U , associated to (a, v). We know that U ∈ C([0, T ]; H s ). Consider the error W ε = U ε -U . It solves

∂ t W ε + d j=1 A j (U ε )∂ j U ε -A j (U )∂ j U = ε 2 LU ε + ε 2 DU ε + ε d j=1 B j (U ε )∂ j U ε .
Rewrite this equation as:

∂ t W ε + d j=1 A j (U ε )∂ j W ε = - d j=1 A j (U ε ) -A j (U ) ∂ j U + ε 2 LW ε + ε 2 LU + ε 2 DW ε + ε 2 DU + ε d j=1 B j (U ε )∂ j U ε .
The operator on the left hand side is the same operator as in (2.1). The term LW ε is not present in the energy estimates, since it is skew-symmetric. The term εLU is considered as a source term: it is of order ε, uniformly in L ∞ ([0, T ]; H s-2 ). Next, the first term on the right hand side is a semi-linear perturbation:

A j (U ε ) -A j (U ) ∂ j U H s-2 A j (U ε ) -A j (U ) H s-2 U H s-1 A j (W ε + U ) -A j (U ) H s-2 C ( W ε L ∞ , U H s-1 ) W ε H s-2 ,
where we have used tame estimates. Finally, we know that

W ε is bounded in L ∞ ([0, T ] × R d ),
as the difference of two bounded terms. The terms involving D and B are treated in a similar way, by adapting the estimates used in the proof of Theorem 1.4. We end up with

d dt SΛ s-2 W ε , Λ s-2 W ε + ε 2 4 ∇Λ s-2 W ε 2 L 2 ε W ε H s-2 + W ε 2 H s-2 ε 2 + SΛ s-2 W ε , Λ s-2 W ε ,
and the result follows from Gronwall lemma.

2.3.

Global existence result in one dimension. In this section, we consider (1.12) with ε > 0 fixed, in dimension d = 1. To prove a global existence result, we shall need finer estimates than (2.3), in particular for low order derivatives.

Proof of Theorem 1.7. For s = 0 and s = 1, (1.12) can be refined, and this will be useful to prove a global existence result.

L 2 estimates. Multiplying the second equation in (1.12) by a ε , integrating in space, and considering the real part yields:

d dt a ε (t) 2 L 2 = 0.
We have simply recovered the conservation of mass (1.5).

Multiplying now the first equation in (1.12) by v ε and integrating in space, we take advantage of the dimensional one:

1 2

d dt v ε (t) 2 L 2 + ε 2 ∂ x v ε (t) 2 L 2 |a ε (t, x)| 2 |∂ x v ε (t, x)|dx. (2.4)
Ḣ1 estimates. Differentiating (1.12) in space, multiplying then the first equation by ∂ x v ε , the second by ∂ x a ε , and integrating, we find (thanks to the symmetrizer):

d dt ∂ x v ε (t) 2 L 2 + 4 ∂ x a ε (t) 2 L 2 + ε 2 ∂ 2 x v ε (t) 2 L 2 ∂ x v ε (t) L ∞ + a ε (t) L ∞ + a ε (t) 2 L ∞ v ε (t) 2 H 1 + a ε (t) 2 H 1 . (2.5)
Now, we have the tools to prove Theorem 1.7. In view of Remark 1.5, it suffices to prove that for any fixed ε > 0,

v ε , ∂ x v ε , a ε , ∂ x a ε , |a ε | 2 ∈ L 1 loc [0, ∞); L ∞ (R d ) .
Since ε > 0 is fixed, we shall omit it in the notations, and fixed it equal to 1 in (1.12). Using Item (iii) of Theorem 1.4, and in view of (1.6), we have the a priori estimate a(t) L 4 1, ∀t 0. The refined L 2 estimate (2.4) for v and Cauchy-Schwarz inequality yield

1 2 d dt v(t) 2 L 2 + ∂ x v(t) 2 L 2 a(t) 2 L 4 ∂ x v(t) L 2 ∂ x v(t) L 2 .
Using Young's inequality

∂ x v(t) L 2 δ ∂ x v(t) 2 L 2 + 1 4δ ,
we infer, for δ > 0 sufficiently small,

d dt v(t) 2 L 2 + ∂ x v(t) 2 L 2 1.
Therefore,

v(t) 2 L 2 + t 0 ∂ x v(τ ) 2 L 2 dτ 1 + t.
The Gagliardo-Nirenberg inequality yields

t 0 v(τ ) L ∞ dτ t 0 v(τ ) 1/2 L 2 ∂ x v(τ ) 1/2 L 2 dτ (1 + t) 1/4 t 0 ∂ x v(τ ) 1/2 L 2 dτ (1 + t) 1/4 t 0 1 + ∂ x v(τ ) 2 L 2 dτ (1 + t) 5/4 ,
where we have used the Young inequality. In view of (1.6), we infer

∂ x a(t) L 2 1 + a(t)v(t) L 2 1 + v(t) L ∞ ,
and the Gagliardo-Nirenberg inequality yields

a(t) L ∞ a(t) 1/2 L 2 ∂ x a(t) 1/2 L 2 1 + v(t) 1/2 L ∞ . We infer t 0 a(τ ) L ∞ + a(τ ) 2 L ∞ dτ (1 + t) 5/4 .

It only remains to prove that

∂ x v, ∂ x a ∈ L 1 loc ([0, ∞); L ∞ (R)). The Ḣ1 estimate (2.5) yields t 0 ∂ 2 x v(τ ) 2 L 2 dτ 1 + t 0 ∂ x v(τ ) L ∞ + 1 + a(τ ) 2 L ∞ a(τ ) 2 H 1 + v(τ ) 2 H 1 dτ 1 + t 0 ∂ x v(τ ) L ∞ + 1 + a(τ ) 2 L ∞ 1 + v(τ ) 2 L ∞ + v(τ ) 2 H 1 dτ 1 + t 0 ∂ x v(τ ) L ∞ + 1 + a(τ ) 2 L ∞ 1 + v(τ ) 2 H 1 dτ,
from Sobolev embedding. Using an integration by parts and the Young inequality, we find

∂ x v(τ ) 2 L 2 v(τ ) L 2 ∂ 2 x v(τ ) L 2 δ ∂ x v(τ ) L ∞ + 1 + a(τ ) 2 L ∞ ∂ 2 x v(τ ) 2 L 2 + ∂ x v(τ ) L ∞ + 1 + a(τ ) 2 L ∞ 4δ v(τ ) 2 L 2 .
Taking δ > 0 sufficiently small, we infer:

t 0 ∂ 2 x v(τ ) 2 L 2 dτ 1 + (1 + t) t 0 ∂ x v(τ ) L ∞ + 1 + a(τ ) 2 L ∞ dτ + (1 + t) t 0 ∂ x v(τ ) L ∞ + 1 + a(τ ) 2 L ∞ 2 dτ 1 + t 2 + (1 + t) t 0 ∂ x v(τ ) 2 L ∞ dτ + (1 + t) t 0 a(τ ) 4 L ∞ dτ.
We also have

t 0 a(τ ) 4 L ∞ dτ t + t 0 v(τ ) 2 L ∞ dτ t + t 0 v(τ ) L 2 ∂ x v(τ ) L 2 dτ t + √ 1 + t t 0 ∂ x v(τ ) L 2 dτ (1 + t) 3/2 ,
and so

t 0 ∂ 2 x v(τ ) 2 L 2 dτ (1 + t) 5/2 + (1 + t) t 0 ∂ x v(τ ) 2 L ∞ dτ.
Now Gagliardo-Nirenberg and Cauchy-Schwarz inequalities yield

t 0 ∂ x v(τ ) 2 L ∞ dτ t 0 ∂ x v(τ ) L 2 ∂ 2 x v(τ ) L 2 dτ t 0 ∂ x v(τ ) 2 L 2 dτ 1/2 t 0 ∂ 2 x v(τ ) 2 L 2 dτ 1/2 √ 1 + t (1 + t) 5/2 + (1 + t) t 0 ∂ x v(τ ) 2 L ∞ dτ 1/2 (1 + t) 7/4 + (1 + t) t 0 ∂ x v(τ ) 2 L ∞ dτ 1/2 .
The Young inequality implies

t 0 ∂ x v(τ ) 2 L ∞ dτ (1 + t) 2 , hence ∂ x v ∈ L 2 loc ([0, ∞); L ∞ ) ⊂ L 1 loc ([0, ∞); L ∞ ). Since ∂ x a ∈ L 1 loc ([0, ∞); L 2 ), Gagliardo-Nirenberg shows that it suffices to prove that ∂ 2 x a ∈ L 1 loc ([0, ∞); L 2 ) to conclude that ∂ x a ∈ L 1 loc ([0, ∞); L ∞ ).
For that, we use again Item (iii) of Theorem 1.4, from which we now that ae iφ ∈ L 1 loc ([0, ∞); H 2 ). Differentiating ae iφ twice, we have, since

∂ x φ = v, ∂ 2 x a(t) L 2 ae iφ H 2 + v∂ x a L 2 + a∂ x v L 2 + v 2 a L 2 ae iφ H 2 + v L ∞ ∂ x a L 2 + a L ∞ ∂ x v L 2 + v 2 L ∞ .
From the above estimate, we infer ∂ 2 x a ∈ L 1 loc ([0, ∞); L 2 ), which completes the proof of the theorem.

Remark 2.1 (Higher dimension). Even though we have obviously taken advantage of the one-dimensional setting to use the embedding H 1 ⊂ L ∞ , the true reason why Theorem 1.7 is limited to d = 1 lies elsewhere. The refined

L 2 estimate (2.4) becomes, if d 2, 1 2 d dt v ε (t) 2 L 2 + ε 2 ∇v ε (t) 2 L 2 |a ε (t, x)| 2 + |v ε (t, x)| 2 |∇v ε (t, x)|dx.

Note the appearance of the new term |v

ε (t, x)| 2 |∇v ε (t, x)|dx: it correspond to the fact that v ε • (v ε • ∇v ε ) is zero if d = 1, but not if d 2 (in general
). This aspect seems to be the only real limitation to extend Theorem 1.7 to d 2.

Numerics

In this section, we define a second order numerical scheme for (1.12). Since ε is no longer a singular perturbation parameter in this reformulation of NLS, this scheme is naturally Asymptotic Preserving for our original problem in the semiclassical limit ε → 0, as long as the solution of the Euler equation remains smooth. Since we solve our problem in dimension 1 on a bounded interval and in dimension 2 on a rectangle, we add periodic boundary conditions to (1.1) or to (1.12), both formulations remaining equivalent.

3.1. The AP numerical scheme. The nonlinear system to be solved reads

     ∂ t a + div(av) + (iε - 1 2 )a div v = i ε 2 ∆a ; a |t=0 = a 0 , ∂ t v + ∇ |v| 2 2 + |a| 2 = ε 2 ∆v ; v |t=0 = v 0 . (3.1)
The semi-discretization in time is realized through a Strang splitting scheme. Let us denote by ∆t n = t n -t n-1 the variable time step, such that t n = n k=1 ∆t k . On each time step [t n , t n+1 ], we split (3.1) into two subsystems and apply the second-order Strang splitting algorithm.

Step 1 for t ∈ [t n , t n + ∆t n+1 2 ]:

∂ t a 1 = i ε 2 ∆a 1 ; a 1|t=tn = a |t=tn , ∂ t v 1 = ε 2 ∆v 1 ; v 1|t=tn = v |t=tn .
Step 2 for t ∈ [t n , t n + ∆t n+1 ]:

       ∂ t a 2 + div(a 2 v 2 ) + (iε - 1 2 )a 2 div v 2 = 0 ; a 2|t=tn = a 1|t=tn+ ∆t n+1 2 , ∂ t v 2 + ∇ x |v 2 | 2 2 + |a 2 | 2 = 0 ; v 2|t=tn = v 1|t=tn+ ∆t n+1 2 .
Step 3 for t ∈ [t n , t n + ∆t n+1 2 ]:

∂ t a 3 = i ε 2 ∆a 3 ; a 3|t=tn = a 2|t=t n+1 , ∂ t v 3 = ε 2 ∆v 3 ; v 3|t=tn = v 2|t=t n+1 .
Then, (a 3 , v 3 )

|t=tn+ ∆t n+1 2
is an approximation of (a, v) |t=t n+1 , solution to (3.1).

This splitting scheme enables to decouple the fluid part (step 2) from the parabolic and Schrödinger parts (step 1 and 3), which allows a standard treatment of each step. Note that for ε = 0, we recover exactly the Euler equations (1.10). We denote by (a n,ε , v n,ε ) the numerical solution at time t n . In the spirit of [START_REF] Holden | Operator splitting for partial differential equations with Burgers nonlinearity[END_REF], the following result can be proven. Proposition 3.1. Under the assumptions of Theorem 1.4, there exist ε 0 > 0 and C, c 0 independent of ε ∈ [0, ε 0 ] such that for all n ∈ N such that t n ∈ [0, T ], where T is as in Theorem 1.4 (ii), for all ∆t n ∈ (0, c 0 ],

|a n,ε | 2 -ρ ε (t n ) L 1 (R d )∩L ∞ (R d ) C(max n ∆t n ) 2 , ε Im(a n,ε ∇a n,ε ) + |a n,ε | 2 v n,ε -j ε (t n ) L 1 (R d )∩L ∞ (R d ) C(max n ∆t n ) 2 .
We assume, without loss of generality, that we are in a periodic framework for the space variable. This allows to solve steps 1 and 3 in a spectral way by using fast Fourier transform, whereas for step 2, we use a Lax-Wendroff scheme (with directional splitting in dimension 2). The extension to Dirichlet or Neumann boundary conditions would require to consider sine or cosine transforms in place of FFT. The total scheme is consistent with (3.1), and second-order in time and space. Since the Lax-Wendroff scheme is explicit, a Courant-Friedrich-Lewy condition is necessary for stability:

∆t n CFL ∆x max j |v n j | + |a n j |
.

where v n j denotes the approximation to v(t n , x j ), x j being a node of the mesh discretization of the computational domain Ω, and ∆x is the space discretization parameter. In all our numerical experiments, we take CFL = 0.8.

Numerical experiments in dimension 1.

In order to analyze the behavior of the numerical solutions (3.1), we proceed like in [START_REF]Numerical study of time-splitting spectral discretizations of nonlinear Schrödinger equations in the semiclassical regimes[END_REF], and compare the behavior of numerical physical observables (1.13) computed by our AP numerical scheme to a reference solution. Since exact solutions to (1.1) with initial data (1.2) are not available, we numerically generated one thanks to the usual splitting method described in [START_REF]Numerical study of time-splitting spectral discretizations of nonlinear Schrödinger equations in the semiclassical regimes[END_REF], applied to (1.1)-(1.2) with a meshing strategy that ensures that the time step ∆t ε and the space step ∆x ε. We consider the initial datum a 0 (x) = e -25(x-0.5) 2 , v 0 (x) = -1 5 ∂ x ln (e 5(x-0.5) + e -5(x-0.5) ), with x ∈ [-1/2, 3/2]. In all our simulations, the smallest value of ε is 10 -4 . We then discretize the x-interval with J = 2 15 subintervals and the time step is ∆t = ε/100. We refer to the reference solution in a generic way by the notation u ε ref .

With the considered initial datum, new oscillations appear in the reference solution passed a time T * ∼ 0.10, meaning that singularities were created in the limit Euler system (1.10). In Figures 1 and2, we plot particle and current densities ρ ε ref and j ε ref for ε = 10 -4 respectively at time t = 0.05 and t = 0.13, before and past formation of shocks. We present on Figure 3 a zoom close to singularities which shows oscillations on physical observables.

In order to show the asymptotic preserving property of our scheme, we compute relative 1 errors at time t n thanks to the formula We evaluate the error function for different values of J = 2 M , where M is an integer here chosen in [START_REF] Bahouri | Fourier analysis and nonlinear partial differential equations[END_REF][START_REF] Carles | [END_REF], and various scaled Planck parameter ε in [10 -4 , 10 0 ]. We plot on Figures 4,5, 6 and 7 the relative error on physical observables ρ ε and j ε computed with our AP schemes at time t = 0.05 and t = 0.13 respectively before and after the formation of singularities. We clearly see that the error is proportional to (∆x) 2 , independently of ε before formation of singularities. Indeed, the mean slope of the error curves with respect to ∆x is close to the value 2 (subfigures (a)).

err ρ ε (t n ) = ρ ε ref -ρ n,ε 1 / ρ ε ref 1 , ρ n,ε j = |a n,ε j | 2 , err j ε (t n ) = j ε ref -j n,ε 1 / j ε ref 1 , j n,ε j = ε Im(a n,ε j ∇a n,ε j ) + ρ n,ε j v n,ε j , where the 1 norm is u 1 = ∆x J-1 k=1 |u j |.
The independence with respect to ε appears in subfigures (b) where error curves are flat. After the formation of shocks, the behavior of our scheme is very different and the mesh parameters have to be reduced as ε → 0 to get good accuracy.

For a comparison, we present the same study for the classical time splitting scheme applied to (1.1)-(1.2) on Figures 8,[START_REF] Besse | Order estimates in time of splitting methods for the nonlinear Schrödinger equation[END_REF][START_REF] Brenier | Convergence of the Vlasov-Poisson system to the incompressible Euler equations[END_REF][START_REF] Buet | Asymptotic preserving and positive schemes for radiation hydrodynamics[END_REF] Thanks to the spectral accuracy of Fast Fourier Transform, the error levels are smaller than for our AP scheme which is only second order with respect to time and space variables. However, contrary to our AP scheme, the error curves clearly depend on ε, the error being of order O(1) when ε is smaller than 2.10 -3 . In order to have good accuracy, it is required to take ∆x ε.

To make our presentation complete, we present the reconstruction of u ε thanks to equation (1.14) with ε = 0.005. We therefore compute the phase φ ε with amplitude a ε and velocity v ε . Since we have access to values of these quantities at every discrete time t n , we approximate the time dependent integral thank to a simple rectangular quadrature. We see that we can recover a very good approximation of the wave function u ε with very few points before the formation of singularities (see Fig. 12). Obviously, one needs more points to have a good reconstruction of the wave function passed the singularities (see Figs. 13 and14). The previous subsection was devoted to one dimensional simulations. We present here the equivalent error analysis for the two dimensional case. The computational domain is now the square [-0.5, 1.5] 2 discretized with J points in each direction. We still need a reference solution generated with Strang splitting scheme applied to equation (1.1)-(1.2). We take J ref = 2 13 = 8192, so ∆x ∼ 2.5 10 -4 and ∆t = ε/100. In a first test, the initial datum is related to the one chosen for the one dimensional case a 0 (x, y) = e -25 r 2 , v 0 (x) = -1 5 ∇ ln (e 5 r + e -5 r ).

where r = (x -0.5) 2 + (y -0.5) 2 . evaluate the error for different values of J = 2 M , where M is chosen in [START_REF] Bahouri | Fourier analysis and nonlinear partial differential equations[END_REF][START_REF] Buet | Asymptotic preserving and positive schemes for radiation hydrodynamics[END_REF], The error analysis results are equivalent to the ones obtained for one dimensional simulation. We recover an independence of the error with respect to ε before the formation of singularities (see Figs 18 and19). We always need to take finer mesh past the shocks formation (see Figs 20 and21). The last computation concerns a non-symmetric solution. The aim of our simulations is to show that we recover good qualitative properties, even for non radially symmetric data. We consider here as initial datum a Maxwell distribution with two temperatures θ 1 = 0.05 and θ 2 = 0.015 for the initial amplitude given by

and
a 0 = 1 √ 4 1 2π √ θ 1 √ θ 2 exp - (x -0.5) 2 2θ 1 - (y -0.5) 2 2θ 2 .
The initial phase is reduced to φ 0 = 0 (note that from (1.11), ∂ t φ |t=0 = -|a 0 | 2 = 0, so φ becomes instantaneously non-trivial). The scaled Planck factor ε is equal to ε = 0.005 and the simulations are performed with J = 2048 intervals in each direction. To distinguish the evolution for each direction (x, y), we make contour plots of wave function and physical observables but also present slice for plane x = 0.5 and y = 0.5 respectively on bottom right and upper right axes. The ). The oscillatory nature is enhanced at time t = 0.08 mainly in y-direction (see Fig. 23).

We recover the good behavior of physical observables before the formation time of singularities. Actually, we compare the particle densities computed by Strang splitting scheme and our AP scheme at time t = 0.035 on Figures 24 and25. The similar representation for j ε is available on Figures 28 and29. Past the formation time of singularities, the contour plots could be thought as equivalent both for particle and current densities. The only noticeable differences can be seen on x = 0.5 and y = 0.5 plane slices (see Fig. 26, 27, 30 and 31).

iε∂ t u ε + ε 2 2 ∆u ε = V ext u ε ; u ε |t=0 = a 0 e iφ 0 /ε , (4.1)
with V ext = V ext (t, x) real-valued. As noticed in [START_REF] Degond | An asymptotic preserving scheme for the Schrödinger equation in the semiclassical limit[END_REF], using the Madelung transform in this context is interesting only if vacuum can be avoided. Having the idea of Grenier [START_REF] Grenier | Semiclassical limit of the nonlinear Schrödinger equation in small time[END_REF] in mind, the counterpart of (1.8) reads

     ∂ t a ε + ∇φ ε • ∇a ε + 1 2 a ε ∆φ ε = i ε 2 ∆a ε , a ε |t=0 = a 0 , ∂ t φ ε + 1 2 |∇φ ε | 2 + V ext = 0, φ ε |t=0 = φ 0 .
We note that the two equations are decoupled: the second equation is of Hamilton-Jacobi type (known as the eikonal equation in this context), and the solution must be expected to develop singularities in finite time (see e.g. [START_REF] Carles | Semi-classical analysis for nonlinear Schrödinger equations[END_REF]), so this approach is doomed to fail for large time. On the other hand, introducing the artificial diffusion as in (1.11), we obtain

     ∂ t a ε + ∇φ ε • ∇a ε + 1 2 a ε ∆φ ε = i ε 2 ∆a ε -iεa ε ∆φ ε , a ε |t=0 = a 0 , ∂ t φ ε + 1 2 |∇φ ε | 2 + V ext = ε 2 ∆φ ε , φ ε |t=0 = φ 0 . (4.2)
The system is still decoupled, but the good news is that the diffusion introduced into the Hamilton-Jacobi equation smoothes the solution out. Therefore, this system can be considered as a good candidate to obtain an AP scheme for the linear Schrödinger equation, regardless of the presence of vacuum. Instead of giving full details of the analogue of Theorems 1.4, 1.6 and 1.7, as well as of their proofs, we shall simply give a functional framework where the viscous eikonal equation can by solved globally in time. Since we consider the case ε > 0 only, we out the ε 2 term in the viscous eikonal equation, and consider

∂ t φ eik + 1 2 |∇φ eik | 2 + V ext = ∆φ eik ; φ eik (0, x) = φ 0 (x). (4.3)
For k 1, let

SL k = {f ∈ C k (R d ; R), ∂ α f ∈ L ∞ (R d ), ∀1 |α| k}.
For k 2, let

SQ k = {f ∈ C k (R d ; R), ∂ α f ∈ L ∞ (R d ), ∀2 |α| k}.
The key result is the following: The solution of this equation does not depend on x. It is given by w(t, x) = 1/(t-1), hence a contradiction. However, in the periodic case x ∈ T d , this technical issue disappears (see Remark 4.3 below).

Proof. The first step consists in constructing the gradient of φ eik . Differentiating The right hand side is a contraction on

w ∈ C [0, T ]; W k-1,∞ , w L ∞ ([0,T ];W k-1,∞ ) 2 ∇φ 0 W k-1,∞ ,
provided that T > 0 is sufficiently small. This follows from the properties of the heat kernel: for t > 0,

e t∆ f L ∞ f L ∞ ; e t∆ ∇f L ∞ C √ t f L ∞ .
The solution to (4.4) is then global, thanks to the maximum principle (see e.g. [47, Proposition 52.8], which implies that v -v v + , where

∂ t v ± ∓ ∇V ext L ∞ x = 0; v ± (0, x) = ± ∇φ 0 L ∞ . Once such a solution v ∈ C(R + ; W k-1,∞ ) is constructed, set φ eik (t) = φ 0 - t 0 1 2
|v(s)| 2 + V ext (s) -div v(s) ds.

We check that ∂ t (∇φ eik -v) = ∇∂ t φ eik -∂ t v = 0. Therefore, v = ∇φ eik , and φ eik solves (4.3). Rewriting (4.3) as 4.2. Other nonlinearities. One might believe that Theorem 1.7 is bound to the cubic one-dimensional Schrödinger equation, which is completely integrable, and that the two aspects are related. We show that this is not the case, since Theorem 1.7 remains valid for other nonlinearities (with d = 1). Consider now

∂ t φ eik + 1 2 |v| 2 + V ext =
iε∂ t u ε + ε 2 2 ∆u ε = f |u ε | 2 u ε , (t, x) ∈ R + × R d . (4.5)
We suppose that there exists δ > 0 such that f (y) δ, for all y 0. Following [START_REF] Grenier | Semiclassical limit of the nonlinear Schrödinger equation in small time[END_REF], the only modification we have to make to recover the results in Section 2.1 consists in replacing the symmetrizer (2.2) with

S = I 2 0 0 1 4f (a 2 1 +a 2
2 ) I d .

From our assumption on f , S and its inverse S -1 are uniformly bounded, provided that a is bounded. As a matter of fact, the exact assumption in [START_REF] Grenier | Semiclassical limit of the nonlinear Schrödinger equation in small time[END_REF] is f > 0, and all the results in Section 2.1 remain valid under this assumption. The more precise assumption f (y) δ becomes useful to prove that in dimension d = 1, the solution is global. The main point to notice is that the conservation of mass and momentum are the same as for (1.1), like in Proposition 1.1. On the other hand, the conservation of energy becomes

d dt ε∇u ε (t) 2 L 2 + R d F |u ε (t, x)| 2 dx = 0,
where F is an anti-derivative of f , F (y) = y 0 f (r)dr.

By assumption, F (y) δ 2 y 2 + f (0)y, so the potential energy controls the L 4 -norm:

u ε (t) 4 L 4 2 δ R d F |u ε (t, x)| 2 dx -f (0) u ε (t) 2 L 2 2 δ R d F |u ε (t, x)| 2 dx -f (0) u ε (0) 2 L 2 ,
where we have used the conservation of mass in the last inequality. This implies an estimate of the form ε∇u ε (t

) 2 L 2 + u ε (t) 4 L 4
C, with C independent of t 0 and ε ∈ (0, 1]. Therefore, the analysis presented in Section 2.3 can be repeated line by line.

Figure 1 :

 1 Figure 1: Particle and current densities at time T * = 0.05 and for ε = 10 -4 .

Figure 2 :

 2 Figure 2: Particle and current densities at time T * = 0.13 and for ε = 10 -4 .

  Current density j ref

Figure 3 :

 3 Figure 3: Zoom around singularity of density and current at time T * = 0.13 and for ε = 10 -4 .

Figure 4 :

 4 Figure 4: err ρ ε (t = 0.05) for AP scheme

  Error w.r.t ε

Figure 5 :Figure 6 :

 56 Figure 5: err j ε (t = 0.05) for AP scheme

Figure 7 :Figure 8 :

 78 Figure 7: err j ε (t = 0.13) for AP scheme

Figure 9 :Figure 10 :Figure 11 : 12 :

 9101112 Figure 9: err j ε (t = 0.05) for time splitting scheme

Figure 13 :

 13 Figure 13: Re(u ε ) at T = 0.13 for ε = 0.005

Figure 14 : 3 . 3 .

 1433 Figure 14: Re(u ε ) at T = 0.13 for ε = 0.005

  various scaled Planck parameter ε in [5.10 -4 , 10 0 ]. The formation time of singularities is as in one dimension situation located around 0.10. The reference particle and current densities for ε = 5.10 -4 are represented on Figures 15 and 17 before and after formation of singularities.We clearly see that a tiny front appears after shocks both in particle and current densities. The wave function presents strong oscillations with new ones created past the shock creation (see Fig.

17

 17

  )

  (a) ρ ref (b) j ref

Figure 15 :

 15 Figure 15: Contour plot of particle and current densities at time t = 0.05 and for ε = 5.10 -4

Figure 16 :

 16 Figure 16: Contour plot of particle and current densities at time t = 0.13 and for ε = 5.10 -4

13 Figure 17 :Figure 18 :

 131718 Figure 17: Contour plot of real part of u ε at times t = 0.05 and t = 0.13 and for = 5.10 -4

Figure 19 :Figure 20 :

 1920 Figure 19: err j ε (t = 0.05) for AP scheme

  Error w.r.t ε

Figure 21 :

 21 Figure 21: err j ε (t = 0.13) for AP scheme

Figure 22 :

 22 Figure 22: Re(u ε )(t = 0.035)

Figure 23 :Figure 24 :Figure 25 :

 232425 Figure 23: Re(u ε )(t = 0.08)

Figure 26 :Figure 27 :Figure 28 :

 262728 Figure 26: Strang ρ ε (t = 0.08)

Figure 29 :Figure 30 :Figure 31 : 4 . 1 .

 29303141 Figure 29: AP j ε (t = 0.035)

Lemma 4 . 1 .

 41 Let k 2 and φ 0 ∈ SL k , V ext ∈ L ∞ loc (R + ; SL k ). Then (4.3) has a unique solution φ eik ∈ C(R + ; SL k-1 ).

Remark 4 . 2 .

 42 No such global result must be expected in the larger class SQ k . Indeed,suppose d = 1, V ext = 0 and φ 0 (x) = -x 2 /2. If (4.3) had a solution φ eik ∈ C(R + ; QL 2 ), then w = ∂ 2 x φ eik ∈ C(R + ; L ∞ ) would solve ∂ t w + v∂ x w + w 2 = ∂ 2x w; w(0, x) = -1.

( 4 . 3 ) 0 e 0 e

 4300 in space, we have to solve∂ t v + v • ∇v + ∇V ext = ∆v; v(0, x) = ∇φ 0 (x). (4.4)This equation is solved locally in time by a fixed point argument by using Duhamel's formulav(t) = e t∆ ∇φ 0 -t (t-s)∆ (v • ∇v)(s)ds -t (t-s)∆ (∇V ext )(s)ds.

0 e 0 eRemark 4 . 3 .

 0043 ∆φ eik ; φ eik (0, x) = φ 0 (x), Duhamel's formula readsφ eik (t) = e t∆ φ 0 -t (t-s)∆ (|v(s)| 2 )ds -t (t-s)∆ V ext (s)ds, from which we conclude that φ eik ∈ C(R + ; SL k-1 ). In the periodic setting x ∈ T d = (R/(2πZ)) d , it is natural to work in the Wiener algebra W =    f : T d → C, f (x) = j∈Z d b j e ij•x with (b j ) j∈Z d ∈ 1 (Z d ) W k = {f : T d → C, ∂ α f ∈ W, ∀|α| k}.Then Lemma 4.1 is easily adapted by replacing SL k with W k .

Example 4 . 4 .

 44 The assumption f δ > 0 is satisfied in the following cases:• Cubic-quintic nonlinearity: f (y) = y + λy 2 , λ 0.• Cubic plus saturated nonlinearity: f (y) = δy + η y 1+λy , δ, η, λ > 0.

Acknowledgements. C. Besse was partially supported by the French ANR fundings under the project MicroWave NT09_460489. R. Carles was supported by the French ANR project R.A.S. (ANR-08-JCJC-0124-01). C. Besse and F. Méhats were partially supported by the ANR-FWF Project Lodiquas (ANR-11-IS01-0003).


Sultan Demirdis Cornelis Jacominus Van Der Beek

Dorothée Colson

Florence Rullier-Albenque

Yanina Fasano

Takasada Shibauchi

Yuji Matsuda

Shigeru Kasahara

Piotr Gierlowski

Marcin Konczykowski

Electron irradiation of Co, Ni, and P-doped BaFe2As2type iron-based superconductors

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The premise of s ± superconductivity in the multiband iron-based superconductors [1,2], with a sign-changing order parameter between the electron-like and hole-like Fermi-surface sheets [2,3,4], has raised the question of the effect of atomic-scale point-like disorder on superconductivity in these materials [5,6]. In particular, interband scattering of quasiparticles by point-like disorder was suggested to be pair-breaking, and, therefore, to lead to a strong suppression of the critical temperature, the appearance of sub-gap states [7], and a peculiar T 2 dependence of the London penetration depth at low temperature [8,9,10,11,12,13,14].

More recently, the question of the role of point-like disorder has become only more relevant [15]. This is because the anti-ferromagnetic coupling between Cooper pairs on different Fermi surface sheets [2,3,4], and the subtle changes of the Fermi surface with doping [3,4], isovalent chemical substitution [17,18,19,20,21,26,27], and pressure [22,23,24,25], allow for a wide variation of realizations of the order parameter [3,4,27]. In particular, the likelihood of line nodes of the order parameter in isovalently doped BaFe 2 (As 1-x P x ) 2 [32], LaFePO [29], and KFe 2 As 2 [28,30,31] are indicative of a nodal s ± -or a d-wave superconductivity [3,4] in the first two and latter respective cases.

Nodal gap structure not imposed by symmetry, but arising from e.g., the anisotropy of the Fermi surface and/or interband scattering [38], has also been evoked [16]. In the case of an s ± state, the order parameter sign would not only change between bands but also within a given band [16]. One can also have an s-wave gap-structure in which nodes appear for non-zero k z [15,39,40], a situation that might be relevant for Ba(Fe 1-x Co x ) 2 As 2 [35], Ba(Fe 1-x Ni x ) 2 As 2 [12], Na 1-δ (Fe 1-x Co x )As [36], and overdoped Ba 1-x K x Fe 2 As 2 [37].

Disorder effect on exotic superconductivity in iron pnictides

Given that the appearance of gap nodes in an s-wave superconductor does not imply symmetry breaking, indirect experiments yield, as yet, most insight into the order parameter structure. Such experiments include the low temperature behaviour of the penetration depth [8,9,10,11,12,13,14,27,29,32] and the thermal conductivity [30,32,35,37] as probes of the low-energy quasi-particle density of states (QPDOS), but also the sensitivity of the iron-based superconductors to quasi-particle scattering by homogeneous atomic-scale point-like disorder. This can be introduced by, e.g., chemical substitution or energetic particle irradiation. Namely, the occurrence of sign changes of the order parameter implies that even non-magnetic point-like scatterers may act as pair-breakers [5,6,7,15]. For s ± -superconductivity in particular, impurity scattering of quasi-particles between electron-like and hole-like bands (with opposite sign of the order parameter) would lead to a rapid suppression of the critical temperature T c [5,6,7], and open the possibility for non-sign changing s ++ superconductivity due to orbital fluctuations to manifest itself [5,6,41]. The vulnerability of the d-wave symmetry to both interbandand intraband impurity scattering would render this even more fragile than s ± . On the other hand, a less pronounced role of interband compared to interband scattering would reduce the sensitivity of s ± superconductivity to disorder [15,41]. Finally, the presence of "accidental" nodes would result in a rather high initial sensitivity of superconductivity to disorder, with a slower suppression of the critical temperature as function of disorder strength once the interband scattering is sufficiently strong to wash out, or "lift" the nodes on all Fermi surface sheets [16].

Effect of native disorder

Even if the problem at hand is characterized by a multitude of parameters, it is important to at least establish trends as controlled point-like disorder is added to iron-based superconductors. Early experiments have suggested a rather limited sensitivity of superconductivity in the ironbased materials to disorder. For example, the residual electronic term in the specific heat suggests a large low-temperature QPDOS [42,43] away from optimal doping, that may be due to pair-breaking. The ubiquitous T 2 -behaviour of the low-temperature ab-plane penetration depth λ ab (T ) suggest strong pair-breaking [8,9,10,11,12,13,14]. Still, T c 's of the materials in question are high.

In our earlier work [44,45], we have drawn attention to the fact that vortex pinning in the mixed state of the iron-based superconductors is indicative of the importance of native disorder. While the critical current density maximum j c at low fields, and the subsequent j c ∼ B -1/2 drop-off are due to nm-scale heterogeneities of the superconducting properties, leading to a socalled "strong-pinning" contribution j s c to the critical current (see Fig. 1 and Refs. [47,48]), the so-called "weak collective pinning" [49,50] contribution j coll of scattering rates of the doping impurities, such as these are obtained from flux pinning, yield rather large values, which are at odds with a superconducting ground state that would be sensitive to point-like disorder should this be pair-breaking. Moreover, it may be noted that the occurrence of the quasi-particle scattering contribution to flux pinning consistently coincides with that of the T 2 behaviour of the penetration depth. Inversely, it is conspicuously absent in clean materials with a T -linear dependence of λ ab , such as BaFe 2 (As 1-x P x ) 2 [32,45,58].

The relative insensitivity of superconductivity in the iron-based materials may be due to these materials not having a sign-changing order parameter, but may also be attributed to limited interband scattering by the dopant impurities. Given the premise of Coulomb scattering (with a small change in crystal momentum) by the charged dopants [45], this is not unreasonable. Estimates from pinning situate the native point defects rather in the Born limit (see table 1).

Artificial disorder

An obvious manner to test the above ideas is through the artificial introduction of different kinds of point-like defects, either by chemical doping [59,60,61], or by irradiation [62,63]. An example is the substitution of magnetic ions, which suppresses T c rather more effectively than non-magnetic substitutions [59,61]. A major drawback of chemical substitution in the iron-based compounds is, however, that this simultaneously leads to structural changes and/or doping. Energetic particle irradiation therefore seems preferable. Recent work by Tarantini et al. reported on T c -suppression by α-particle irradiation of NdFeAs(O,F), but was prone to criticism, notably in that it induces magnetic impurities and defect clustering. The 3 MeV proton irradiation experiments by Nakajima et al. on Ba(Fe 1-x Co x ) 2 As 2 reported monotonous T c -depression as the residual resistivity increases; the authors provided estimates for the critical pair-breaking parameter as Γ/2πT c ∼ 4-7 and Γ/2πT c ∼ 1.5-2.5, as one goes from underdoped to overdoped [63]. However, apart from point defects, proton and neutron irradiation are known to induce, in situ, point defect cascades and clusters, which may play a different role than that of simple scatterers.

Here, we compare the effect of 2.5 MeV electron irradiation on differently substituted Table 1. Contribution of dopant disorder to elastic scattering parameters of various iron pnictide superconductors, such as estimated from the weak collective pinning contribution to the (flux pinning) critical current density, j coll c . Here, k F is the Fermi wavevector, ξ 0 is the Bardeen-Cooper-Schrieffer coherence length, n d is the atomic point defect density, σ tr is the transport scattering cross-section, δ 0 is the scattering phase angle, Γ = n d [πN (0)] -1 sin 2 δ 0 is the scattering rate, where N (0) = mk F /2πh 2 is the normal-state DOS. BaFe 2 As 2 . Such high-energy electron irradiation is known to produce homogeneously distributed vacancy-interstitial (Frenkel) pairs on all sublattices of a crystal, but, preferentially on that constituted by the lightest atomic species in the material -in the materials studied here, Fe.

Material impurity k F ξ 0 n d σ tr n d ξ 3 0 sin δ 0 Γ Å-1 nm nm -3 Å2
In the high-T c cuprates, electron-irradiation defects are known to be strong unitary scatterers, comparable to Zn substitution. There, they are responsible for the suppression of the critical temperature [64,65,66], and the appearance of a T 2 -temperature dependence of the penetration depth [66]. In semi-metals with extremely low carrier densities, such as Bi, the irradiationinduced defects are charged and lead to effective doping of the material [67]. The mobility of the created point defects at high temperature is a caveat, since it leads to partial annealing and clustering, albeit it to a lesser extent than that produced by proton or neutron irradiation.

-1 10 Single crystalline Ba(Fe 0.925 Co 0.075 ) 2 As 2 : Hysteresis loop of the irreversible magnetization, at 5 K, before (•) and after ( ) irradiation with 5.5×10 19 cm -2 2.5 MeV electrons. Closed arrows depict the direction in which the loop is traversed upon cycling the magnetic field.

Dotted double arrows depict the width of the magnetization loop in the low-field, strong pinning regime, and in the higher field regime in which only the weak collective pinning contribution j coll c is relevant.

It turns out that electron irradiation leads to similar monotonic suppression of T c in Ba(Fe 1-x Co x ) 2 As 2 , Ba(Fe 1-x Ni x ) 2 As 2 , and BaFe 2 (As 1-x P x ) 2 . In Ba(Fe 1-x Co x ) 2 As 2 , the suppression is weakest for the optimally doped material and strongest for underdoped crystals. In the range of intermediate electron fluences studied here, the resistance monotonically increases as function of fluence. In Ba(Fe 1-x Co x ) 2 As 2 and Ba(Fe 1-x Ni x ) 2 As 2 , there is little to no change of the temperature dependence of the superfluid density. We have also measured the critical current density before and after the irradiation. In Ba(Fe 1-x Co x ) 2 As 2 one has a clear increase of the weak collective pinning contribution to j c , while in BaFe 2 (As 1-x P x ) 2 the initially absent weak pinning contribution emerges after irradiation. This finding establishes that the fieldindependent contribution j coll c to the critical current density of iron-based superconductors is, indeed, due to atomic sized point pins. In both materials, this finding allows for an estimate of the density of produced defects per unit irradiation fluence. In Ba(Fe 1-x Ni x ) 2 As 2 , j c is strongly suppressed by the irradiation.

Experimental details 2.1. Single crystalline samples

Single-crystals of Ba(Fe 1-x Co x ) 2 As 2 [68] and Ba(Fe 1-x Ni x ) 2 As 2 [69] were grown using the selfflux method. The high purity starting reagents Ba, FeAs, and CoAs were mixed in the molar ratio 1 : (4 -x) : x, loaded in alumina crucibles and then sealed in evacuated quartz tubes. For each doping level, chemical analysis by an electron probe was performed on several crystals yielding the dopant content within 0.5% absolute accuracy. For this work we studied six Codoping levels and two Ni-doping levels. The optimally doped BaFe 2 (As 1-x P x ) 2 single crystals, with x = 0.33 and x = 0.36, were also grown by the self flux method [70], and characterized using x-ray diffraction and energy dispersive x-ray spectroscopy. No impurity phases were found within the experimental limits of ∼ %1.

Electron irradiation

The 2.5 MeV electron irradiation was performed at the SIRIUS Pelletron facility of the Laboratoire des Solides Irradiés (LSI) at the Ecole Polytechnique in Palaiseau, France [71]. Beam currents varied between 10 and 20 µA, depending on the irradiation run; the beam was swept over a 6 × 6 mm 2 area. Samples of different composition irradiated to the same dose were mounted together and irradiated simultaneously. In order to prevent in-situ defect migration, recombination, and clustering, the irradiations were carried out in a liquid H 2 bath (T = 20 K).

In-situ resistance measurements were carried out on the Ba(Fe 1-x Co x ) 2 As 2 material in order to assess the increase of the residual resistivity as a function of defect density at low temperature, as well as the effect of room-temperature annealing.

Measurements of T c and the surface impedance

T c -values before and after irradiation were determined using differential magneto-optical (DMO) imaging (see Fig. 2 and Ref. [47]). For the lower irradiation doses, in particular, the shift in T c after irradiation is smaller than the transition width. DMO then provides an objective means to quantify the effect of the electron irradiation on different portions of a crystal. The change in T c was further quantified using ex-situ measurements of the temperature dependent microwave surface impedance Z s = R s + iX s . These were performed at the LSI using a superconducting Nb cavity, operating at 26.0 GHz, and cooled to 5 K using a 0.5 W cryocooler cold head. From the temperature dependent resonance frequency f and unloaded quality factor Q 0 , the surface resistance R s and reactance X s were determined as

R s = µ 0 f G 1 Q 0 (T ) - 1 Q0 (1)
and

X s (T ) = µ 0 G [f (T ) -f (0)] + X 0 (2)
respectively. Here f (0) and Q0 are the resonance frequency and the unloaded quality factor in the absence of the sample, and the geometrical factor G = V /4πw √ wd, with V the inner volume of the cavity, d the thickness of the platelet-like single crystal sample, and w its smaller width. The additive constant X 0 was adjusted so that R s = X s in the normal state.

Determination of the critical current density

Measurements of the magnetization M of the Ba(Fe 1-x Co x ) 2 As 2 and Ba(Fe 1-x Ni x ) 2 As 2 single crystals as function of the applied magnetic field H a were performed using a commercial SQUID magnetometer. From the width ∆M (H) of the magnetic hysteresis loops, the fieldand temperature dependence of the sustainable current density j = 3 2 ∆M/w was determined using the Bean model. The numeric prefactor 3 2 is estimated from calculations of E.H. Brandt for rectangular bars of similar aspect ratio as the measured crystals [72]. Data on BaFe 2 (As 1-x P x ) 2 were acquired by M. Konczykowski using the local Hall probe magnetometry technique [73,74]. Briefly, the single crystalline superconducting crystal is centered on top of an array of 11 Hall sensors, of active area 3×3 µm 2 , and spaced by 20 µm, in such a manner that the short er crystal edge is parallel to the array, and the array spans the longer crystal edge. The Hall sensors are used to measure the spatial gradient dB/dx of the local induction component perpendicular to the crystal surface. The gradient dB/dx ( in G/µm ) is directly proportional to the sustainable screening current j, with dB/dx ∼ 1 6 j for a crystal aspect ratio of ∼ 0.1 [72]. Measurements of the flux density gradient were performed both as function of H a after zero field-cooling, and as function of temperature T , after field cooling, and the subsequent reduction of H a by twice the full penetration field.

In what follows, we shall take the screening current at the lowest measurement temperature ( 5 K ) to be representative of the critical current density. At higher temperature, the measured current density can be significantly lower than the pinning critical current density because of flux creep [58,74].

Analysis of the critical current density

The critical current density is further analyzed along the lines of Refs. [44,45]. The hysteresis loops of all charge-doped iron-based superconductors, exemplified by Fig. 1, indicate the presence of two main pinning mechanisms in these materials. The central peak of the magnetization, at low applied fields, is due to the strong-pinning contribution j s c to the critical current, while the constant contribution at higher fields comes from "weak collective pinning" by the dopant atoms. In a critical current density -versus magnetic field plot (see Fig. 6) the strong pinning contribution is responsible for the low-field plateau and the subsequent j c ∝ B -1/2 decrease, while the collective pinning contribution yields the high-field constant j c . It can be written in terms of superconducting parameters, defect density n d , and the elementary pinning force f p of a single defect, as The crystal is then warmed to 300 K, cooled to 21 K, and heated once again to 300 K. After the first anneal at 300 K, the resistance drops back to the value reached after low-temperature irradiation with only 1.4 Ccm -2 (white circle on the curve in the main panel).

j coll c = j 0 n d f 2 p ξ 3 ab ε λ ε 0 2/3 . (3)
Here, j 0 is the depairing current density, ξ ab and λ ab are the ab-plane coherence length and penetration depth respectively, ε 0 ≡ Φ 2 0 /4πµ 0 λ 2 ab is the vortex line energy, ε λ is the penetration depth anisotropy [51], and Φ 0 = h/2e is the flux quantum. The averaging . . . is performed over the vortex core.

The magnitude of j coll c is compatible with scattering of the quasi-particles in the vortex cores as being at the origin of the weak collective pinning contribution in pristine iron-based superconductors, provided that the dopant atoms are the scattering defects. This gives rise to [52,53]

f p = 0.3g(ρ D )ε 0 σ tr πξ 2 ab ξ 0 ξ ab (4)
which depends not only on the Gor'kov impurity parameter ρ D = hv F /2πT c l = ξ 0 /l, but also on the transport cross-section σ tr = (2π/k 2 F ) sin 2 δ 0 (δ 0 is the scattering phase angle, and k F the Fermi wavevector). The analysis of the intermediate field critical current density of several common iron based superconductors allows one to estimate the scattering parameters compiled in Table 1. The weak collective pinning contribution is absent in isovalently doped BaFe 2 (As 1-x P x ) 2 for all doping levels x [see Ref. [58] and Fig. 9(a)]; surprisingly, it is present in Ba(Fe 1-x Ru x ) 2 As 2 [54] and LiFeAs [73].

Results and Discussion

4.1. Ba(Fe 1-x Co x ) 2 As 2 Fig. 4 shows the resistance of an optimally doped Ba(Fe 0.925 Co 0.075 ) 2 As 2 single crystal, measured in-situ during exposure to the 10 µA, 2.5 MeV beam, at T = 21 K. The sample transits to the normal state at a fluence of 1.2 Ccm -2 , after which the resistance increases at a rate of δR/R = 0.14[Ccm -2 ] -1 . After a fluence of 4.6 Ccm -2 was reached, the electron beam was stopped, the sample heated to 300 K, and the resistance measured. The crystal was subsequently cooled down once again to 21 K, and the resistance measured again (Inset to Fig. 4). It appears that after the 300 K anneal, the resistance has dropped to the value first measured after irradiation with a fluence of 1.4 Ccm -2 . Therefore, annealing at 300 K leads to the annihilation and clustering of the point defects produced by the irradiation, with a resulting drop of the irradiation-induced resistance change of 65 %. Since all data presented below concern ex-situ experiments performed after heating the crystals to 300 K, the effect of annealing should be taken into account.

The evolution of the reduced critical temperature T c /T c0 with electron fluence (after 300 K annealing) is shown for various Co-substitution levels x in Fig. 5(a). Here, T c0 is the critical temperature of the crystal before irradiation. It is to be stressed that, given the modest T c changes involved, and the typical dispersion of T c -values between similarly grown samples, even from the same batch, it is essential to compare the values of T c and T c0 on the same crystal before and after irradiation. For the same reason, a plot of the absolute T c versus dose would be meaningless in the regime of small doses. Fig. 5 also shows the effect of 300 K annealing, for one sample. It is seen that the change in T c drops by a factor three after annealing. Here m e is the free electron mass, and m * the effective mass. zδΓ = eδρ/R H m e (with e the electronic charge) is estimated from the change of the resistance due to the irradiation (after annealing), together with the published data for the resistivity ρ and the Hall coefficient R H published in Ref. [68]. In spite of the contributions of multiple bands to the conductivity of Ba(Fe 1-x Co x ) 2 As 2 , this approach can be justified given that the conductivity contribution from the electron sheets dominates the transport properties [68]. Typically, δR/R ∼ 0.05 [Ccm -2 ] -1 for the optimally doped sample (cf Fig. 3). Note, however, that the resistivity change with electron dose depends on x, and that, in order to produce the "simulated" δΓ values of Fig. 5b, we use data measured in-situ on crystals with different x. Estimates of the scattering rate from critical current measurements (see below) yield qualitatively similar results as those shown in the Figure .. The surface impedance of optimally doped Ba(Fe 1-x Co x ) 2 As 2 before and after electron irradiation is depicted in Fig. 3. A monotonous increase of the surface resistance is observed after irradiation; even though in this series of experiments, the increase was not linear as function of fluence. In particular, the crystal irradiated with 2.1 Ccm -2 electrons did not conform to the general trend. The normalized frequency shift of the Nb cavity, which is proportional to the superfluid density n s ∝ λ -2 ab of the superconducting crystals, showed a slight change to a less marked temperature dependence, indicative of disorder having been added to the crystal. The trend is similar to that observed by Hashimoto et al. in Ba 1-x K x Fe 2 As 2 crystals with varying degrees of disorder [75], but is much less marked here.

The low-temperature critical current density of single crystalline Ba(Fe 0.925 Co 0.075 ) 2 As 2 increases by a constant, field-independent amount after electron irradiation (see Fig. 6a). This indicates that it is the weak collective pinning contribution j coll c that is enhanced, presumably by atomic-sized point defects introduced by the irradiation. If one assumes that the relevant defects are Fe vacancy-interstitial pairs, i.e. that σ tr ∼ πD 2 F e ∼ 2.6 Å2 , with D F e = 0.92 Å the ionic radius of Fe 2+ , the increase of the critical current density is consistent, through Eqs. ( 3) and ( 4), with 0.006 dpa Fe / Ccm -2 . The increase of the low-temperature critical current density and the concomitant decrease of T c after irradiation implies a steeper decrease of j c with temperature, and a crossing of j c -values of the pristine and the irradiated crystal at an intermediate temperature.

Ba(Fe

B ( mT ) position ( µm ) µ 0 H a =0.1 T Figure 7.
Temperature dependence of the flux density gradient dB/dx in a BaFe 2 (As 0.64 P 0.36 ) 2 single crystal before and after irradiation with 0.8 × 10 19 electrons cm -2 . Inset: Flux density profile across the crystal after zero-field cooling, application of the applied field µ 0 H a = 100 mT, and subsequent warming. T c versus dose is comparable to that found in underdoped Ba(Fe 1-x Co x ) 2 As 2 (see Fig. 11). As in Ba(Fe 1-x Co x ) 2 As 2 , the temperature dependence of the superfluid density of Ba(Fe 1-x Ni x ) 2 As 2 changes little or not at all, even at the largest irradiation dose of 5.5 × 10 19 electrons cm -2 .

The critical current density of the Ba(Fe 1-x Ni x ) 2 As 2 crystals is strongly suppressed after the irradiation (see Fig. 6b). This is partially because the reduced measurement temperatures are higher in Ba(Fe 1-x Ni x ) 2 As 2 than in its co-substituted counterpart, but mainly because the suppression of the prefactor j 0 ∝ n s e(∆/hk F ) in Eq. ( 3) outweighs the increase in n d (∆ is the superconducting gap amplitude).

BaFe 2 (As

1-x P x ) 2
The Inset to Fig. 7 shows an example of flux density gradients measured on a BaFe 2 (As 0.64 P 0.36 ) 2 single crystal irradiated with 0.9×10 19 electrons cm -2 , after zero field cooling and the application of an external field of 100 mT. The temperature dependence of the flux density gradient obtained after field-cooling in 400 mT, reduction of the applied field to the indicated target value, and subsequent warming is shown in the main panel of the Figure . The flux density gradient is directly proportional to the sustainable current density j. As in Ba(Fe 1-x Co x ) 2 As 2 , the screening current in the irradiated crystal exceeds that of the pristine crystal at low temperature, but drops below it at higher temperature, due to the decrease of T c .

The critical temperatures after electron-irradiation were determined as the extrapolation to zero of the "irreversibility field" H irr (T ) above which the nonlinearity of the current-voltage characteristic -and therefore the critical current density -vanishes. Details on the determination of H irr , can be found in, e.g., Ref. [74]. In all cases, the H irr (T ) values depend very little on the frequency of the ac field used to investigate the screening by the superconducting sample; the corresponding T irr (H) lies very close to the temperatures at which the dc screening current vanishes in Fig. 7. The T c -values resulting from the extrapolation of H irr (T ) to zero are gathered in Fig. 11. The T c -depression as function of electron fluence is comparable to that in Ba(Fe 1-x Ni x ) 2 As 2 . Note that the resistivity increase of the irradiated BaFe 2 (As 1-x P x ) 2 crystals corresponds to δR/R ∼ 0.16 [Ccm -2 ] -1 , three times higher than for Ba(Fe 1-x Co x ) 2 As 2 .

Figure 9 shows the hysteretic loops of the local flux density gradient on the surface of a pristine (a) and electron-irradiated BaFe 2 (As 1-x P x ) 2 single crystal (b), for different temperatures. The pristine single crystal shows astroid-shaped magnetic hysteresis loops that are indicative of strong flux pinning by extended point-like defects only. After electron irradiation, the hysteresis loops open up due to the appearance of the field-independent collective pinning contribution from atomic-sized point defects. We interpret the appearance of a non-zero j coll c as being due to vortex pinning by vacancy-interstitial pairs introduced by the irradiation. Moreover, j coll c monotonically increases (at low temperature) as function of irradiation fluence (Fig. 10), which allows for a direct comparison with the theory for quasi-particle scattering mediated collective vortex pinning. Under the assumption that the only unknown parameter, σ tr , is again approximated by the ionic cross-section of a Fe-vacancy, a fit to Eq. (3) yields a defect density of 0.0035 dpa / Ccm -2 . This number is taken as more precise than that obtained for Ba(Fe 1-x Co x ) 2 As 2 in subsection 4.1, in which only two points were available. Note that the hypothesis of simple voids (non-scattering point defects) would necessitate an unphysical 1 dpa / Ccm -2 to explain the magnitude of the critical current density change.

Discussion; Comparison between differently substituted materials

The modification of the critical current density of Ba(Fe 1-x Co x ) 2 As 2 and BaFe 2 (As 1-x P x ) 2 by electron irradiation demonstrates that the field-independent contribution to j c apparent in sub-T to T-fields is indeed due to collective pinning by atomic-sized point defects. The very reasonable number obtained for defect generation by the electron beam, 0.35% dpa / Ccm -2 , indicates that it is, in a first approximation, safe to assume that the most relevant produced defects are Fe vacancies, and that these vacancies are responsible for quasi-particle scattering. This interpretation is strengthened by the fact that scattering rates in Ba(Fe 1-x Co x ) 2 As 2 estimated from the change of the resistivity are very much comparable to Γ = n d /[πN (0)]/ sin 2 δ 0 estimated from the defect density. A comparison with Table 1, shows that the scattering cross-section of Fe vacancies should be comparable to that of Co and Ru impurities. Indeed, electronirradiation of Ba(Fe 1-x Co x ) 2 As 2 and Ba(Fe 1-x Ni x ) 2 As 2 does not lead to appreciable changes of the temperature dependence of the superfluid density of those materials, which suggests that the effect of the dopant impurities overwhelms that of the defects added by the irradiation.

Irrespective of the precise changes in resistivity incurred by the irradiation, a striking feature of our results is the similar suppression of T c by the electron irradiation of underdoped Ba(Fe 1-x Co x ) 2 As 2 , Ba(Fe 1-x Ni x ) 2 As 2 and BaFe 2 (As 1-x P x ) 2 for a given dose. The T csuppression in optimally doped Ba(Fe 1-x Co x ) 2 As 2 is a factor 2 less. Given the clear evidence for line nodes in BaFe 2 (As 1-x P x ) 2 , this poses the question possible nodal structures of the order parameter in Ba(Fe 1-x Ni x ) 2 As 2 , and, possibly, Ba(Fe 1-x Co x ) 2 As 2 as well. In the latter two materials, it was proposed that a nodal line might exist on the α-(hole-like) sheet, at finite k z [35,12].

It is tempting to compare our results to that obtained by other types of irradiation, as well as by doping. It turns out that the suppression of T c for optimally doped Ba(Fe 1-x Co x ) 2 As 2 is remarkably comparable to that reported for Co-doping of Ba 1-x K x Fe 2 As 2 in Ref. [76], and much less than that obtained by Co-doping of KFe 2 As 2 studied by the same authors. The decrease of T c observed after electron irradiation is also comparable, as far as optimally doped Ba(Fe 1-x Co x ) 2 As 2 is concerned, to the Co-doping experiments on Ba 1-x K x Fe 2 As 2 by Li et al., while the T c suppression in underdoped Ba(Fe 1-x Co x ) 2 As 2 , Ba(Fe 1-x Ni x ) 2 As 2 , and BaFe 2 (As 1-x P x ) 2 is similar to the Cu and Zn doping by those authors. Even if the comparison with chemical doping is tenuous, the observed trend clearly indicates a much weaker sensitivity of the materials studied here than what is expected for the scenario of s ± superconductivity with strong interband scattering, and sets the three materials even further from the d-wave scenario.

A comparison with the previous irradiation studies reveals that T c -suppression after electron irradiation is somewhat weaker than the results obtained by Nakajima by 3 MeV proton irradiation of optimally doped Ba(Fe 1-x Co x ) 2 As 2 [63], and much weaker than in the α-particle irradiation of NdFeAs(O,F) by Tarantini et al. [62].

Summary and Conclusions

Ba(Fe 1-x Co x ) 2 As 2 crystals of different dopng levels x, as well as Ba(Fe 1-x Ni x ) 2 As 2 and BaFe 2 (As 1-x P x ) 2 crystals, have been subjected to high energy (2.5 MeV) electron irradiation at low temperature (21 K). In spite of substantial annealing and clustering of the defects, the enhancement of the collective pinning contribution to the critical current density due to quasi-particle scattering in the vortex cores in Ba(Fe 1-x Co x ) 2 As 2 , and its appearance in BaFe 2 (As 1-x P x ) 2 in which it is absent in the pristine material, demonstrate that the irradiation produces atomic-sized point defects. The increase in magnitude of the low-temperature critical current density is consistent with the pinning by Fe vacancies, created by the irradiation. In Ba(Fe 1-x Co x ) 2 As 2 and Ba(Fe 1-x Ni x ) 2 As 2 , the introduction of these supplementary point defects does not lead to a significant modification of the temperature dependence of the superfluid density, which suggests that the defects are overwhelmed by pre-existing disorder (notably, the presence of the dopant atoms). Surprisingly, the critical temperature is similarly suppressed in all three materials. This, in spite of the fact that the order parameter in BaFe 2 (As 1-x P x ) 2 is thought to have line nodes, while this possibility is much less certain in the other two materials. The results therefore lend credence to evidence for line nodes obtained from c-axis penetration depth [12] and thermal conductivity measurements [35].

  Figure 1.Single crystalline Ba(Fe 0.925 Co 0.075 ) 2 As 2 : Hysteresis loop of the irreversible magnetization, at 5 K, before (•) and after ( ) irradiation with 5.5×10 19 cm -2 2.5 MeV electrons. Closed arrows depict the direction in which the loop is traversed upon cycling the magnetic field.Dotted double arrows depict the width of the magnetization loop in the low-field, strong pinning regime, and in the higher field regime in which only the weak collective pinning contribution j coll

Figure 2 .

 2 Figure 2. Transition from the superconducting to the normal state of a Ba(Fe 0.925 Co 0.075 ) 2 As 2 single crystal before (a,c) and after irradiation with 2.1 Ccm -2 2.5 MeV electrons (b,d), as imaged by the differential magneto-optical (DMO) technique [47]. The images (a,b) show the progressive admission of an ac magnetic field of magnitude 1 Oe, applied perpendicularly to the sample surface, as the temperature is raised. In these flux density maps, areas of high luminous intensity I correspond to the value of the applied field, while dark areas correspond to zero field, (i.e. full screening). The dark rectangles in the upper left hand panels of (a) and (b) correspond to the sample outline, i.e., full Meissner expulsion of the magnetic field at the lowest temperature. The white squares indicate the positions at which the "transmittivity" (or "local ac susceptibility") defined as [I(T ) -I(T T c )]/[I(T T c ) -I(T T c )] depicted in the lower panels (c,d) was determined.

Figure 3 .

 3 Figure 3. Surface impedance Z s of Ba(Fe 0.925 Co 0.075 ) 2 As 2 crystals before and after various low-temperature irradiation runs with 2.5 MeV electrons, and subsequent annealing at 300 K. Data points and thin lines show the surface resistance R s and reactance X s respectively, for various electron doses. (b) Shift δf (T ) of the resonant frequency of the superconducting Nb cavity (operating at 26.0 GHz) as function of temperature, normalized with respect to the low-temperature extrapolated δf (0), for Ba(Fe 0.925 Co 0.075 ) 2 As 2 crystals irradiated with various fluences of 2.5 MeV electrons.

Figure 4 .

 4 Figure 4. Ba(Fe 0.925 Co 0.075 ) 2 As 2 : Resistance as function of dose of 2.5 MeV electrons, measured at 21 K. The sample transits to the normal state after a dose of 1.2 Ccm -2 , after which the resistance increases at a rate of δR/R = 0.14[Ccm -2 ] -1 . The Inset shows the effect of annealing. The resistance of the crystal is represented as function of time. The protocol comprises initial measurements at 21 K following irradiation with 4.6 Ccm -2 2.5 MeV electrons.The crystal is then warmed to 300 K, cooled to 21 K, and heated once again to 300 K. After the first anneal at 300 K, the resistance drops back to the value reached after low-temperature irradiation with only 1.4 Ccm -2 (white circle on the curve in the main panel).

Figure 5 .

 5 Figure 5. Ba(Fe 1-x Co x ) 2 As 2 : evolution of T c (normalized by the initial critical temperature T c0 ) as function of (a) the fluence of 2.5 MeV electrons (b) the estimated normal state scattering rate Γ (normalized by T c0 ), for different doping levels x. The small data points show the T c -values after annealing at 300 K. The large green triangle shows the drop of T c following low-temperature irradiation of optimally doped Ba(Fe 0.925 Co 0.075 ) 2 As 2 , the arrow shows the effect of annealing at 300 K.

Fig. 5 (

 5 b) shows the same data, as function of the product of the induced change in the scattering rate δΓ/2πT c and the electronic effective mass enhancement z = m * /m e .

Figure 6 .

 6 Figure 6. Sustainable current density j of (a) optimally doped Ba(Fe 1-x Co x ) 2 As 2 (x = 0.075) and (b) optimally doped Ba(Fe 1-x Ni x ) 2 As 2 ( x = 0.035) before (open symbols) and after irradiation at 21 K with 5.5 × 10 19 cm -2 2.5 MeV electrons (closed symbols). The measurements were performed at 5 K ( • , • ), 7.5 K ( , ), 10 K ( , ), 11 K ( , ), 13 K ( ), 16 K ( ), and 17.5 K ( , ).

  1-x Ni x ) 2 As 2 Surface impedance measurements and DMO imaging on the irradiated Ba(Fe 1-x Ni x ) 2 As 2 crystals show a more rapid depression of T c with dose than in Ba(Fe 1-x Co x ) 2 As 2 . If anything,

Figure 8 . 2 (

 82 Figure 8. Irreversibility field BH irr (T ) as function of temperature, for BaFe 2 (As 0.64 P 0.36 ) 2 single crystals irradiated with 0.8 × 10 19 ( ), 1.6 × 10 19 (•), and 4.5 × 10 19 electrons cm -2 (•, , ). B irr (T ) is determined as the locus of temperatures T above which the nonlinearity of the crystal current-voltage characteristic vanishes at given induction H a . T c is determined from the extrapolation of H irr (T ) to zero field.

Figure 9 .

 9 Figure 9. (a) Hysteretic loops of the local flux density gradient versus local induction B, measured on the surface of a pristine BaFe 2 (As 0.67 P 0.33 ) 2 single crystal, at various temperatures (indicated). The astroid-shaped hysteresis loops are determined by the sole strong-pinning contribution to the critical current. (b) ibid, measured on the surface of a BaFe 2 (As 0.7 P 0.3 ) 2 single crystal irradiated at 23 K with 1.8 × 10 19 electrons cm -2 , at the indicated temperatures. The opening of the loops at higher flux densities reveal the emergence of a weak collective pinning contribution by the atomic-sized point defects introduced by the irradiation.

Figure 10 .

 10 Figure 10. BaFe 2 (As 1-x P x ) 2 (with x = 0.33, x = 0.36): Electronfluence dependence of the weak collective pinning contribution j coll c to the critical current density. The drawn line denotes a fit to Eqs.(3,4). Assuming that the point defects most relevant for quasi-particle scattering are Fe vacancies induced by the irradiation, the parameter value 5.1 × 10 8 would correspond to 0.0035 dpa / Ccm -2 .

2 x 36 Figure 11 .

 23611 Figure 11. (a) Dose-dependence of the T c -change of Ba(Fe 1-x Co x ) 2 As 2 , Ba(Fe 1-x Ni x ) 2 As 2 ,and BaFe 2 (As 1-x P x ) 2 , after irradiation with 2.5 MeV electrons and annealing at 300 K. The upper scale shows the density of point defects (presumably Fe vacancies) added by the irradiation, such as determined from the dose-dependence of the critical current density of BaFe 2 (As 1-x P x ) 2 (see Fig.10). (b) Dose dependence of T c /T c0 , where T c0 is the critical temperature of the pristine crystal.

Acknowledgements

We acknowledge support from the the grant MagCorPnic of the Réseau Thématique de Recherche Avancée "Triangle de la Physique du Plateau de Saclay", and by the Agence Nationale de la Recherche grant "PNICTIDES". Part of the work was also possible thanks to the support of the ECOS-Sud-MINCyt France-Argentina bilateral program, Grant No.A09E03. Work in Japan was supported by KAKENHI from JSPS, and by Grant-in-Aid for the Global COE program The Next Generation of Physics, Spun from Universality and Emergence from MEXT, Japan.