Laurent Gosse
email: l.gosse@ba.iac.cnr.it

A well-balanced and asymptotic-preserving scheme for the one-dimensional linear Dirac equation

Keywords: Dirac equation, One-dimensional relativistic quantum mechanics, Asymptotic-preserving and well-balanced numerical methods Mathematics Subject Classification (2000) MSC 35Q41, 65M06, 81Q05

The numerical approximation of one-dimensional relativistic Dirac wave equations is considered within the recent framework consisting in deriving local scattering matrices at each interface of the uniform Cartesian computational grid. For a Courant number equal to unity, it is rigorously shown that such a discretization preserves exactly the L 2 norm despite being explicit in time. This construction is wellsuited for particles for which the reference velocity is of the order of c, the speed of light. Moreover, when c diverges, that is to say, for slow particles (the characteristic scale of the motion is non-relativistic), Dirac equations are naturally written so as to let a "diffusive limit" emerge numerically, like for discrete 2-velocity kinetic models. It is shown that an asymptotic-preserving scheme can be deduced from the aforementioned well-balanced one, with the following properties: it yields unconditionally a classical Schrödinger equation for free particles, but it handles the more intricate case with an external potential only conditionally (the grid should be such that c∆ x → 0). Such a stringent restriction on the computational grid can be circumvented easily in order to derive a seemingly original Schrödinger scheme still containing tiny relativistic features. Numerical tests (on both linear and nonlinear equations) are displayed.

Introduction

The fundamental equation of relativistic quantum mechanics is the Dirac equation which combines the special theory of relativity with quantum effects; see also [START_REF] S M Morsink | Black hole radiation of Dirac particles in 1+1 dimensions[END_REF][START_REF] Succi | Lattice Boltzmann equation for quantum mechanics[END_REF] for an application within the Jackiw-Tetelboim "R=T" model of 1 + 1 general relativity. In one space dimension, the dynamics of such spin- 1 2 particles is described by a 2-component spinor wave-function but nevertheless display interesting features.

Scalar and vector potentials in linear 1+1 Dirac systems

In 1D, the Dirac equation for a relativistic particle of mass m > 0 moving in an external potential V (x) ∈ R 2 reads, where H 0 stands for the "free Hamiltonian",

i∂ t Ψ = [H 0 + V (x)]Ψ , H 0 = -icσ 1 ∂ x • + mc 2 h σ 3 • .
The vector Ψ (t, x) ∈ C 2 is usually called a Dirac spinor. By neglecting magnetic fields, the remaining external potential V is split into:

V = V t + σ 1 V e + σ 2 V p + σ 3 V s ,
where V e , the so-called space-term can always be absorbed by means of a gauge transform. The other components V t , V p and V s are referred to as the time-potential, scalar and pseudo-scalar terms, respectively. Hereafter, the focus will be drawn onto the case V p ≡ 0, V s ,V t are real scalar functions of x ∈ R. To fix ideas, the term V s may correspond to a confining potential, like a space-dependent effective mass of an electron [START_REF] Hiller | Solution of the one-dimensional Dirac equation with a linear scalar potential[END_REF], whereas V t , the so-called Lorentz time-potential cannot have a confining effect because of the Klein paradox [START_REF] Bosanac | Solution of Dirac equation for a step potential and the Klein paradox[END_REF] (see also [START_REF] Domingez-Adame | Solvable linear potentials in the Dirac equation[END_REF][START_REF] Domingez-Adame | A one-dimensional relativistic screened Coulomb potential[END_REF][START_REF] Solomon | An exact solution of the Dirac equation for a time-dependent Hamiltonian in 1-1 dimension space-time[END_REF][START_REF] Succi | Numerical solution of the Schrödinger equation using discrete kinetic theory[END_REF][START_REF] Succi | Lattice Boltzmann equation for quantum mechanics[END_REF][START_REF] Weinberger | All you need to know about the Dirac equation[END_REF]). Parameters c > 0 and σ i , i = 1, 2, 3, denote the speed of light and classical Pauli matrices, respectively:

σ 1 = 0 1 1 0 , σ 2 = 0 -i i 0 , σ 3 = 1 0 0 -1 .
1.2 Conservation of the L 2 -norm and non-relativistic limit

In the context of quantum mechanics applications, the preservation of the L 2 norm is important because it is related to the conservation of the probability of presence for particles. For the linear Dirac equations, the Crank-Nicolson scheme is a popular choice [START_REF] Alvarez | Linearized Crank-Nicolson scheme for nonlinear Dirac equations[END_REF][START_REF] Alvarez | The numerical study of a nonlinear one-dimensional Dirac equation[END_REF] because, by Cayley's theorem [START_REF] Hairer | Geometric numerical integration. Structure-preserving algorithms for ordinary differential equations[END_REF][START_REF] Guo | On the Relation between Stability of Continuous-and Discrete-Time Evolution Equations via the Cayley Transform[END_REF], it preserves exactly this L 2 norm: here we show in §2.2 that the well-balanced approach [START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws[END_REF] allows to improve on what is explained in [START_REF] Wessels | Discretizing the one-dimensional Dirac equation[END_REF] by furnishing an original, less costly, explicit scheme endowed with the same preservation property. An original MUSCL-type reconstruction [START_REF] Gosse | MUSCL reconstruction and Haar wavelets[END_REF] is presented, too, in §2. [START_REF] Askar | Explicit integration method for the timedependent Schrödinger equation for collision problems[END_REF]. This is numerically illustrated on both linear and nonlinear benchmarks: see §4.1 and §4.2, respectively. Usually, such well-balanced discretizations deliver automatically asymptotic-preserving (AP) schemes when the limiting process only retains flux terms, like in [START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws[END_REF][START_REF] Gosse | An asymptotic-preserving well-balanced scheme for the hyperbolic heat equations[END_REF]: here, by studying the nonrelativistic limit of Dirac equations with several potentials [START_REF] Bechouche | Semi-(non)relativistic limits of the Dirac equation with external time-dependent electromagnetic fields[END_REF], we see that this is not true anymore when these zero-order potential terms survive in the limit. Indeed, it is shown in §3.3 and §3.4 that despite the second order finite difference is unconditionally recovered, obtaining the correct contribution of the potential asks for a stringent restriction on the computational grid. The somewhat simpler case of non-relativistic limit of free particles (including explicit schemes) is studied in §3.3. However, it is possible to circumvent it and the correct asymptotic behavior is achieved. Numerical validations are shown in §4.3. Existing asymptotic-preserving strategies were mostly developed in order to handle the classical limit of both linear and non-linear Schrödinger equations: see [START_REF] Carles | Numerical aspects of the nonlinear Schrödinger equation in the semiclassical limit in a supercritical regime[END_REF][START_REF] Degond | An asymptotic preserving scheme for the Schrodinger equation in the semiclassical limit[END_REF]. a strong convergence result for free particles is furnished in Appendix A.

2 Derivation of L 2 -conservative well-balanced approximations

It appears that the presentation advocated in [START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws[END_REF], and which relies on the use of local scattering matrices, is closely related to ideas published in other contexts, see [START_REF] Bahlouli | Solution of one-dimensional Dirac equation via Poincaré map[END_REF][START_REF] Cotaescu | Applying the Dirac equation to derive the transfer matrix for piecewise constant potentials[END_REF].

Computation of local scattering matrices

Accordingly, we start with a model slightly more involved than what appears in [START_REF] Wessels | Discretizing the one-dimensional Dirac equation[END_REF]:

∂ t Ψ + c 0 1 1 0 ∂ x Ψ = - i h m(x)c 2 1 0 0 -1 + V t (x) Ψ . (2.1)
Denoting Ψ = (ψ + , ψ -), one diagonalizes the convective part of (2.1) by means of Riemann invariants, which reduce to diagonal variables by linearity. The resulting problem reads, for h = 1 and c 2 being lumped into the space-dependent mass m(x),

ζ ± = 1 2 (ψ + ± ψ -), ∂ t ζ ± ± c∂ x ζ ± = -i (m(x)ζ ∓ + V t (x)ζ ± ) , (2.2)
In order to derive a set of piecewise-constant numerical approximations ζ ∆ x ± , one deduces from [START_REF] Wessels | Discretizing the one-dimensional Dirac equation[END_REF] that the right-hand side must be upwinded if the explicit Euler time-integrator is used. Following [START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws[END_REF], Chapter 8, the first step for building a wellbalanced scheme for (2.2) is to consider the wave equation for which low-order terms are "localized" on a uniform grid with a space-step ∆ x:

∂ t ζ ± ± c∂ x ζ ± = -i∆ x ∑ j∈Z (m(x)ζ ∓ + V t (x)ζ ± ) δ (x -x j-1 2 ), (2.3)
where shorthand notation x j = j∆ x is used for j ∈ Z, together with t n = n∆t, n ∈ N and ∆t is a time-step constrained by a Courant number ν = c∆t ∆ x . Piecewise constant approximations ζ ∆ x ± are defined by values generated by a time-marching process:

ζ ∆ x ± (t n , x j ) = ζ n j,± ∈ C.
The presence of Dirac masses at each interface x j-1 2 yields the following scheme:

ζ n+1 j,+ = ζ n j,+ - c∆t ∆ x (ζ n j,+ -ζ n j-1 2 ,+ ), ζ n+1 j,-= ζ n j,-+ c∆t ∆ x (ζ n j+ 1 2 ,--ζ n j,-). (2.4)
Coupling zero-order terms are treated by means of jump relations induced by the Dirac masses inserted in (2.3). Formal consistency comes from a weak convergence,

x j x j-1 2 x j+ 1 2 ζ n j,± ζ n j+1,± ζ n j-1,± ζ n j-1,+ ζ n j,- ζ n j,+ ζ n j+1,- ζ n j-1 2 ,+ ζ n j-1 2 ,- ζ n j+ 1 2 ,- ζ n j+ 1 2 ,+ 1
∆ x ∑ j∈Z δ (x -x j-1 2 ) = ∆ x ∑ j∈Z δ x -( j - 1 2 )∆ x ⇀ 1 as ∆ x → 0.
This way, interface values ζ n j∓ 1 2 ,± are deduced from ζ n j∓1,± by means of local scattering matrices (see Fig. 2.1) which expression is given below. Lemma 2.1 For any j ∈ Z, ∆ x > 0, the smooth solution ζ± (x) of the boundary-value problem for the stationary equations in x ∈ (0, ∆ x), with ζ+ (0

) = ζ L + , ζ-(∆ x) = ζ R -, ±c∂ x ζ± = -i m(x j-1 2 ) ζ∓ + V t (x j-1 2 ) ζ± , satisfies the following relation, with constant values V t = V t (x j-1 2 ) c , ω = m(x j-1 2 ) c : ζ+ (∆ x) ζ-(0) = exp(i∆ xV t ) cosh(∆ xω) -i tanh(∆ xω) -i tanh(∆ xω) exp(-i∆ xV t ) cosh(∆ xω) ζ+ (0) ζ-(∆ x) . (2.5)
The 2 × 2 matrix in (2.5) is the scattering matrix for (2.4): it relates the outgoing states ζ∓ (0/∆ x) to the incoming ones ζ± (0/∆ x). These states are said to be incoming because they are prescribed on incoming characteristics pointing inside the domain.

One may factorize this matrix with sech(ω∆ x) = 1 cosh(ω∆ x) , the hyperbolic secant:

S(x) = sech(ω∆ x) exp(i∆ xV t ) -i sinh(∆ xω) -i sinh(∆ xω) exp(-i∆ xV t ) , S j± 1 2 def = S(x j± 1 2
).

Proof: The computation is standard because the stationary equations are linear, only the forward/backward boundary data make it differ from a usual ODE system:

c d dx ζ+ ζ- = -i V t ω -ω -V t ζ+ ζ- .
By means of a diagonalization process and a logarithmic integration, one finds:

ζ+ (∆ x) = exp(i∆ xV t ) ζ+ (0) cosh(ω∆ x) -i ζ-(0) sinh(ω∆ x) , ζ-(∆ x) = exp(i∆ xV t ) i ζ+ (0) sinh(ω∆ x) + ζ-(0) cosh(ω∆ x) . (2.6)
In the second equation, the unknown is ζ-(0), hence the system becomes:

ζ-(0) = exp(-V t ∆ x) cosh(ω∆ x) ζ-(∆ x) -i tanh(ω∆ x) ζ+ (0), ζ+ (∆ x) = exp(i∆ xV t ) ζ+ (0) cosh(ω∆ x) -sinh 2 (ω∆ x) cosh(ω∆ x) -i ζ-(∆ x) tanh(ω∆ x).
By rearranging terms, one finds easily (2.5).

Hyperbolic rotation, L 2 -estimate and convergence

In the special case for which V t ≡ 0, one observes that the relations (2.6) imply that ζ± (∆ x) (right states) result from the action of a hyperbolic rotation onto ζ± (0) (left states). Planar hyperbolic rotations are sometimes called "Lorentz transformations". Going back to the general case V t ≡ 0, the scheme (2.4) rewrites:

ζ n+1 j,+ = (1 -ν)ζ n j,+ + ν exp(i∆ xV t (x j-1 2
))

cosh(∆ xω j-1 2 ) ζ n j-1,+ -i tanh(∆ xω j-1 2 )ζ n j,- ζ n+1 j,-= (1 -ν)ζ n j,-+ ν exp(-i∆ xV t (x j+ 1 2 )) cosh(∆ xω j+ 1 2 ) ζ n j+1,--i tanh(∆ xω j+ 1 2 )ζ n j,+ (2.7)
Let's hereafter denote S j-1 2 the 2 × 2 scattering matrix appearing in (2.4): in the particular case where the Courant number ν = 1, which means that the computational grid is chosen so that c∆t = ∆ x (like in [START_REF] Wessels | Discretizing the one-dimensional Dirac equation[END_REF]), ζ n+1 j,± are the scattering values from the data of ζ n j∓1,± . Thus a simple criterion for the preservation of the L 2 norm reads: Lemma 2.2 Assume V t (x) ∈ L ∞ ∩L 2 (R) is a real function and c∆t = ∆ x, there holds:

∀n ∈ N, ζ ∆ x (t n , •) L 2 (R) = ζ ∆ x (0, •) L 2 (R) .
Remark 2.1 The preservation of the L 2 norm with an explicit scheme like (2.7) is not straightforward: indeed, stable upwind discretizations with Courant numbers ν < 1 always contain a certain amount of numerical viscosity which, in the nonlinear case, dissipate enough entropy to prevent oscillation onset in the numerical solution. Mathematical entropies being just convex functions, in the context of quantum mechanics applications which ask for conservation of the probability of presence, the scheme must moreover preserve (and not dissipate) a particular entropy, the L 2 -norm.

Proof: The property is equivalent to stating that for any j ∈ Z, the scattering matrix S j-1 2 is unitary. By multiplying it by the complex conjugate of its transpose, it is straightforward to check that one recovers the identity matrix of R 2 . Hence:

ζ n+1 j,+ ζ n+1 j-1,- * ζ n+1 j,+ ζ n+1 j-1,- = ζ n j,+ ζ n j-1,- * S * j-1 2 T S j-1 2 ζ n j,+ ζ n j-1,- , which directly yields |ζ n+1 j,+ | 2 + |ζ n+1 j-1,-| 2 = |ζ n j,+ | 2 + |ζ n j-1,-| 2 .
To derive the uniform estimate, it remains to multiply by ∆ x and to sum on j ∈ Z.

⊓ ⊔ Thanks to the linearity of the system (2.1), the aforementioned conservation of the L 2 norm, which yields weak compactness in L 2 , is enough to ensure the weak convergence of a subsequence of approximate solutions ζ ∆ x ± toward the exact ones [START_REF] Thaller | Advanced Visual Quantum Mechanics[END_REF] for smooth V t . Moreover, by uniqueness of the limit, all the sequence converges. More accurate bounds are established in §A where strong L 2 loc convergence is shown.

A MUSCL-type reconstruction compatible with the S-matrix

The L 2 -preservation property rigorously holds only for a Courant number ν = 1, thus one may want to make the scheme (2.7) less vulnerable to the damage of numerical dissipation. One standard way to proceed is to set up the so-called MUSCLreconstruction, relying on a conveniently limited piecewise-linear extrapolation of the cell-average values ζ n j,± . One issue is that such piecewise-linear data ask now for the resolution of "Generalized Riemann Problems" (GRP) at each interface x j± 1 2 , see [START_REF] Berthon | Space-time Generalized Riemann Problem Solvers of Order k for Linear Advection with Unrestricted Time Step[END_REF][START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction[END_REF]: even in cases where it is doable, it doesn't fit with treating the right-hand side by means of S-matrices. However, a new type of MUSCL scheme was devised in [START_REF] Gosse | MUSCL reconstruction and Haar wavelets[END_REF] and, along with heavy notation, it proceeds as follows:

1. Reconstruction: given the sequence of spinors (ζ n j,± ) j∈Z ∈ C 2 , one first performs a linear interpolation. At both borders of each cell

C j = (x j-1 2 , x j+ 1 2
), the following values are derived by means of local slopes σ n j and a slope-limiter φ : R → [0, 2],

ζ n,R j-1 2 ,± = ζ n j,± - σ n j ∆ x 2 , ζ n,L j+ 1 2 ,± = ζ n j,± + σ n j ∆ x 2 .
More precisely, mass conservation is ensured with the choice:

ζ n,R j-1 2 ,± = ℜ(ζ n j,± ) - φ (r n j ) 2 ℜ ζ n j+1,± -ζ n j,± +iℑ(ζ n j,± ) -i φ (r n j ) 2 ℑ ζ n j+1,± -ζ n j,± , φ (r) = r + |r| 1 + |r| ,
and the ratios to which slope-limiters are applied are,

r n j = ℜ ζ n j+1,± -ζ n j,± ℜ ζ n j,± -ζ n j-1,± , rn j = ℑ ζ n j+1,± -ζ n j,± ℑ ζ n j,± -ζ n j-1,±
.

Finally, the reconstructed interface states read, for each cell C j :

ζ n,R j-1 2 ,± = ζ n j,± - φ (r n j ) 2 ℜ ζ n j+1,± -ζ n j,± -i φ (r n j ) 2 ℑ ζ n j+1,± -ζ n j,± , ζ n,L j+ 1 2 ,± = ζ n j,± + φ (r n j ) 2 ℜ ζ n j+1,± -ζ n j,± + i φ (r n j ) 2 ℑ ζ n j+1,± -ζ n j,± .
2. Fluxes computation: Yet we have discontinuities at both the center x j and interfaces x j± 1 2 : assuming the right-hand side of (2.2) is again localized only at interfaces, the formulation (2.3) still holds. Thus the discontinuities at the center x j are simply advected until they reach the interfaces x j± 1 2 , and there they are scattered by the S-matrix: see Fig. 2.2. Observe that (2.7) can be put in a matrix form as:

ζ n+1 j,+ ζ n+1 j-1,- = (1 -ν) ζ n j,+ ζ n j-1,- + νS j-1 2 ζ n j-1,+ ζ n j,- , ν = c∆t ∆ x . (2.8)
So, by taking into account the (new) discontinuity at the cell's center,

ζ n+1 j,+ ζ n+1 j-1,- = ζ n j,+ ζ n j-1,- -min(ν, 1 2 )      ζ n,L j+ 1 2 ,+ ζ n,R j-1 2 ,-   -S j-1 2   ζ n,L j-1 2 ,+ ζ n,R j+ 1 2 ,-      + max(0, ν -1 2 )      ζ n,R j+ 1 2 ,+ ζ n,L j-1 2 ,-   -S j-1 2   ζ n,R j-1 2 ,+ ζ n,L j+ 1 2 ,-      .
3. Overall algorithm: Hence, by linearity, one can form the following 2 states,

ζ n j+ 1 2 ,+ = min(ν, 1 2 )ζ n,L j+ 1 2 ,+ -max(0, ν -1 2 )ζ n,R j+ 1 2 ,+ , ζ n j+ 1 2 ,-= min(ν, 1 2 )ζ n,R j+ 1 2 ,- -max(0, ν -1 2 )ζ n,L j+ 1 2 ,- , (2.9)
and the MUSCL scheme rewrites in a slightly less involved form:

ζ n+1 j,+ ζ n+1 j-1,- = ζ n j,+ ζ n j-1,- -      ζ n j+ 1 2 ,+ ζ n j-1 2 ,-   -S j-1 2   ζ n j-1 2 ,+ ζ n j+ 1 2 ,-      .
(2.10)

The scheme (2.9)-(2.10) will be tested in §4.2.2.

Asymptotic-Preserving property for non-relativistic limit

The numerical time-marching scheme (2.7) is well adapted to the case where c is of the order of the characteristic scale of the computational grid, that is to say, for relativistic particles. In case the particles (usually represented with wavepackets endowed with a particular group velocity) move much slower, the aforementioned computational grids, for which c∆t ≃ ∆ x, don't really fit thus one is interested to explore the behavior of the numerical process (2.7) in the limit c → +∞.

Formal non-relativistic limit of 1D Dirac equation

Hereafter we consider the particular case of (2.1) with a mass m > 0 independent of x. In order to study the non-relativistic limit of a Dirac particle with (positive) energy E ≃ mc 2 , one retains only the slower time-variation like exp(-i(Emc 2 )t/h) by considering the rescaled quantities, ψ± (t, x) = exp(imc 2 t/h)ψ ± (t, x), which satisfy:

∂ t ψ+ + c∂ x ψ-= i hV t (x) ψ+ , ∂ t ψ-+ c∂ x ψ+ = 2i mc 2 h ψ-+ i hV t (x) ψ-.
When c diverges, the second equation is dominated by the rule, ψ-= h 2imc ∂ x ψ+ , which, once inserted in the first one, yields the classical Schrödinger equation:

i∂ t ψ+ + h2 2m ∂ xx ψ+ = V t (x) ψ+ . (3.1)
The rigorous theory of non-relativistic limits of linear multidimensional Dirac equation with external potential is given in [START_REF] Bechouche | Semi-(non)relativistic limits of the Dirac equation with external time-dependent electromagnetic fields[END_REF][START_REF] Hiller | Solution of the one-dimensional Dirac equation with a linear scalar potential[END_REF].

Rescaled diagonal variables and local scattering matrix

With a slight abuse of notation, we still denote the Riemann invariants associated with the "rescaled spinor" ζ ± = 1 2 ( ψ+ ± ψ-). A new (small) parameter ε = 1 c is introduced for ease of reading, and a rescaled system is derived:

∂ t ζ ± ± 1 ε ∂ x ζ ± = 2i ± m ε 2 (ζ + -ζ -) + V t (x) 2 ζ ± . (3.2)
By substituting V t 2 → V t and following the same roadmap as in Lemma 2.1, one gets:

Lemma 3.1 Let V t > 0, ε > 0
λ ε = (εV t + m/ε) 2 -(m/ε) 2 → 2V t m, ε → 0,
1. the local scattering matrix reads:

S = 1 D ε -im ελ ε sin(2λ ε ∆ x) 1 1 -im ελ ε sin(2λ ε ∆ x) . (3.3) where D ε = cos(2λ ε ∆ x) -i εV t +m/ε λ ε sin(2λ ε ∆ x). 2.
Moreover, S is unitary: (S * ) T S is the 2 × 2 identity matrix.

Proof: The proof of (i) proceeds by integrating the following ODE system:

1 ε d dx ζ+ ζ- = 2i m ε 2 + V t -m ε 2 m ε 2 -( m ε 2 + V t ) ζ+ ζ- , x ∈ (0, ∆ x).
By means of a standard diagonalization procedure, one finds that ±λ ε are corresponding eigenvalues and a tedious logarithmic integration delivers:

ζ+ (∆ x) ζ-(∆ x) = cos(2λ ε ∆ x) + i εV t +m/ε λ ε sin(2λ ε ∆ x) -im ελ ε sin(2λ ε ∆ x) im ελ ε sin(2λ ε ∆ x) cos(2λ ε ∆ x) -i εV t +m/ε λ ε sin(2λ ε ∆ x) ζ+ (0) ζ-(0) . Since ζ-(∆ x)
belongs to the incoming boundary data, the second line must be inverted in order to derive the expression of the outgoing state ζ-(0):

ζ-(0) = ζ-(∆ x) -im ελ ε sin(2λ ε ∆ x) ζ+ (0) cos(2λ ε ∆ x) -i εV t +m/ε λ ε sin(2λ ε ∆ x) .
When plugging ζ-(0) into the expression of ζ+ (∆ x), a simplification occurs:

ζ+ (∆ x) = =1 cos 2 (2λ ε ∆ x) + =1 (εV t + m/ε) 2 λ 2 ε - (m/ε) 2 λ 2 ε sin 2 (2λ ε ∆ x) cos(2λ ε ∆ x) -i εV t +m/ε λ ε sin(2λ ε ∆ x) ζ+ (0) - im ελ ε sin(2λ ε ∆ x) cos(2λ ε ∆ x) -i εV t +m/ε λ ε sin(2λ ε ∆ x) ζ-(∆ x).
It is now straightforward to derive the expression of the scattering matrix. Now, for proving (ii), one just computes (S * ) T S and simplifies the expression of D ε • D * ε :

(S * ) T S = 1 cos 2 (2λ ε ∆ x) + (εV t + m/ε) 2 (εV t + m/ε) 2 -(m/ε) 2 =1+ (m/ε) 2 λ 2 ε sin 2 (2λ ε ∆ x) ×   1 + (m/ε) 2 λ 2 ε sin 2 (2λ ε ∆ x) 0 0 1 + (m/ε) 2 λ 2 ε sin 2 (2λ ε ∆ x)   .
With this expression, one deduces easily that (S * ) T S is the 2 × 2 identity matrix.

The special case of free particles

In the former calculation, it appears that a source of complication is the interaction of both source terms, namely the particle's mass m and the time-potential V t . However, one may restrict attention to the simpler case of non-relativistic limit of free particles, for which V t ≡ 0. Clearly, this assumption implies that λ ε = 0, hence cos(2λ ε ∆ x) = 1 and εV t +m/ε

λ ε sin(2λ ε ∆ x) = 2∆ xm ε .
The (constant) scattering matrix S simplifies into:

S = 1 1 -2im∆ x/ε 1 -2im∆ x/ε -2im∆ x/ε 1 .

Asymptotic-Preserving implicit schemes

One can follow the early computations in [START_REF] Gosse | An asymptotic-preserving well-balanced scheme for the hyperbolic heat equations[END_REF] dealing with the diffusive regime: it consists in observing that the outgoing states rewrite in the equivalent form,

   ζ+ (∆ x) = ζ-(∆ x) + 1 1-2i∆ xm/ε ζ+ (0) -ζ-(∆ x) , ζ-(0) = ζ+ (0) + 1 1-2i∆ xm/ε ζ-(∆ x) -ζ+ (0) .
Inserting these scattering states into the expression of the scheme (2.4) yields:

   ζ n+1 j,+ = ζ n j,+ -∆t ε∆ x (ζ n+1 j,+ -ζ n+1 j,-) + ∆t ε∆ x-2i∆ x 2 m ζ n+1 j-1,+ -ζ n+1 j,-, ζ n+1 j,-= ζ n j,-+ ∆t ε∆ x (ζ n+1 j,+ -ζ n+1 j,-) + ∆t ε∆ x-2i∆ x 2 m ζ n+1 j+1,--ζ n+1 j,+ , (3.4)
In order to shed light onto the numerical dynamics of ψn j,+ , one sums up both lines of this scheme and unconditionally obtains that, in the limit ε → 0,

i ψn+1 j,+ = i ψn j,+ - ∆t 2∆ x 2 m ψn+1 j+1,+ -2 ψn+1 j,+ + ψn+1 j-1,+ ,
a consistent finite-difference approximation of the classical Schrödinger equation, treated implicitly because the explicit Euler scheme is unconditionally unstable [START_REF] Chan | Stable Explicit Schemes for Equations of the Schrdinger Type[END_REF].

Remark 3.1 Based on the results of [START_REF] Alvarez | Linearized Crank-Nicolson scheme for nonlinear Dirac equations[END_REF][START_REF] Alvarez | The numerical study of a nonlinear one-dimensional Dirac equation[END_REF][START_REF] Lubich | Integrators for Quantum Dynamics: A Numerical Analyst's Brief Review, in Quantum Simulations of Complex Many-Body Systems: From Theory to Algorithms[END_REF][START_REF] Markowich | Numerical approximation of quadratic observables of Schrödinger-type equations in the semiclassical limit[END_REF], it may be a better choice to use a second-order Crank-Nicolson time-integrator (instead of an implicit Euler method) in (3.4) in order to preserve the L 2 -norm of the (piecewise constant) numerical approximation ψ∆x + . The computations become more involved, though, as one can see:

               ζ n+1 j,+ + ∆t 2ε∆ x (ζ n+1 j,+ -ζ n+1 j,-) - ∆t/2 ε∆ x-2i∆ x 2 m ζ n+1 j-1,+ -ζ n+1 j,- = ζ n j,+ -∆t ε2∆ x (ζ n j,+ -ζ n j,-) + ∆t/2 ε∆ x-2i∆ x 2 m ζ n j-1,+ -ζ n j,-, ζ n+1 j,--∆t 2ε∆ x (ζ n+1 j,+ -ζ n+1 j,-) - ∆t/2 ε∆ x-2i∆ x 2 m ζ n+1 j+1,--ζ n+1 j,+ = ζ n j,-+ ∆t 2ε∆ x (ζ n j,+ -ζ n j,-) + ∆t/2 ε∆ x-2i∆ x 2 m ζ n j+1,--ζ n j,+ .
(3.5)

An explicit and conditionally stable AP scheme ?

The derivation of Asymptotic-Preserving schemes is usually all about the behavior, as a small parameter goes to zero, of discrete derivatives in the space variable, x. The discretization in time being kept as the Euler forward method, for instance. In the present context, this raises an issue because applying straightforwardly such an elementary time-integrator leads to unconditional instability: in the former subsection, we thus chose to set up time-implicit schemes for the sake of stability, but a drawback is their CPU cost. Hence one may wonder whether it may be possible to proceed, like in [START_REF] Gosse | An asymptotic-preserving well-balanced scheme for the hyperbolic heat equations[END_REF], in such a way that a stable time-explicit scheme for the Schrödinger equation can be recovered in the limit ε → 0. Basically, explicit, or barely implicit (i.e. not asking for inversions of linear systems) methods essentially reduce to:

-Leapfrog or Dufort-Frankel methods, which are not self-starting, [START_REF] Askar | Explicit integration method for the timedependent Schrödinger equation for collision problems[END_REF][START_REF] Markowich | Numerical approximation of quadratic observables of Schrödinger-type equations in the semiclassical limit[END_REF][START_REF] Wu | Dufort-Frankel-Type Methods for Linear and Nonlinear Schrödinger Equations[END_REF]] -Euler forward method, but perturbed by supplementary dissipation, [START_REF] Chan | Stable Explicit Schemes for Equations of the Schrdinger Type[END_REF][START_REF] Mickens | Stable explicit schemes for equations of Schrödinger type[END_REF] Setting up the AP process described in (3.4) within either Leapfrog or Dufort-Frankel time-integration methods appears to produce only numerical instability. However, the second strategy, based on supplementary diffusion terms can lead to stable algorithms, so we present it now, following both [START_REF] Chan | Stable Explicit Schemes for Equations of the Schrdinger Type[END_REF][START_REF] Mickens | Stable explicit schemes for equations of Schrödinger type[END_REF].

1. the following forward-Euler scheme (compare with (3.4)) is unstable,

     ζ n+1 j,+ -ζ n j,+ ∆t + 1 ε∆ x (ζ n+1 j,+ -ζ n+1 j,-) = 1 ε∆ x-2i∆ x 2 m ζ n j-1,+ -ζ n j,-, ζ n+1 j,--ζ n j,- ∆t -1 ε∆ x (ζ n+1 j,+ -ζ n+1 j,-) = 1 ε∆ x-2i∆ x 2 m ζ n j+1,--ζ n j,+ , (3.6)
at least because, by formally imposing ε = 0 and summing, an unstable discretization of the free Schrödinger equation arises. However, it has the advantage of being easy to implement since the implicit part can be inverted analytically. 2. one may think that this instability comes exclusively from the change of equation's type in the limit ε → 0, and so that modifying the time-differentiation like in [START_REF] Chan | Stable Explicit Schemes for Equations of the Schrdinger Type[END_REF][START_REF] Mickens | Stable explicit schemes for equations of Schrödinger type[END_REF] can fix this issue. Keeping finite-differences in space, one substitutes:

ζ n+1 j,± -ζ n j,± ∆t ζ n+1 j,± -ζ n j,± ∆t 1 + ∆t(α -iβ ) or ζ n+1 j,± -ζ n j,± ∆t + ∆ x 2 (α -iβ ) .
According to [START_REF] Chan | Stable Explicit Schemes for Equations of the Schrdinger Type[END_REF], the second choice leads to CFL restrictions in ∆t ≃ O(∆ x 4 ), which is prohibitively costly, so we restrict ourselves to the first one, less demanding on ∆t, even if it can become unstable for ∆t, ∆ x becoming very small. 3. within our notation, there is supplementary diffusion if β > 0, but it is needed only for small values of ε, so we propose the following modified time-differentiation,

ζ n+1 j,± -ζ n j,± ∆t 1 + ∆t(α -iβ )(1 -ε γ ) , 0 < γ ≤ 1, (3.7)
meaning that diffusion is activated only as ε ≪ 1. It remains to plug it inside (3.6).

Position densities were computed with the resulting scheme (3.6)-(3.7), β = 1 and ε = 1, 0.1, 0.01, see Fig. 3.1. Initial data was sampled on a uniform grid of 2 8 points,

ζ ± (t = 0, x) = 1 2 exp i tanh(x) -(x/4) 2 ,
x ∈ (-16, 16).

The data for Schrödinger, treated with the second-order Leapfrog method, is thus 2ζ ± (t = 0, •). Numerical position densities become closer as ε is decreased, up to ε ≃ 0.001 where convergence stalls because of diffusion. Accordingly, one cannot really speak about explicit AP schemes because the scheme obtained in the limit ε → 0 doesn't match exactly the numerical approximation one gets, for instance, out of a Leapfrog discretization: it still contains the supplementary diffusion which is necessary for stability. Hence we shall not push the study of such conditionally stable explicit Asymptotic-Preserving discretizations any further.

Difficulties with the complete case

In the general case V t ≡ 0, the full scattering matrix (3.3) must be considered, and the eigenvalue λ ε strongly couples both the source terms m/ε 2 and εV t : in order to ease the reading, we shall limit ourselves to the case of a constant time-potential V t > 0 because it already highlights most of the corresponding issues. Thanks to this simplifying assumption, the scattering matrix S is still independent of the index j ∈ Z (the modifications necessary for handling a space-dependent time-potential are straightforward) and the outgoing states read now:

   ζ+ (∆ x) = ζ-(∆ x) + 1 D ε ζ+ (0) -cos(2λ ε ∆ x) -iεV t λ ε sin(2λ ε ∆ x) ζ-(∆ x) , ζ-(0) = ζ+ (0) + 1 D ε ζ-(∆ x) -cos(2λ ε ∆ x) -iεV t λ ε sin(2λ ε ∆ x) ζ+ (0) .
Clearly, these expressions split naturally between a first part which generates the second-order divided difference (like in the preceding section), and a second one which contains the time-potential term. Plugging into the upwind scheme (2.4),

                     ζ n+1 j,+ = ζ n j,+ -∆t ε∆ x (ζ n+1 j,+ -ζ n+1 j,-) + ∆t εD ε ∆ x ζ n+1 j-1,+ -cos(2λ ε ∆ x)ζ n+1 j,- + i∆t εD ε ∆ x εV t λ ε sin(2λ ε ∆ x)ζ n+1 j,- ζ n+1 j,-= ζ n j,-+ ∆t ε∆ x (ζ n+1 j,+ -ζ n+1 j,-) + ∆t εD ε ∆ x ζ n+1 j+1,--cos(2λ ε ∆ x)ζ n+1 j,+ + i∆t εD ε ∆ x εV t λ ε sin(2λ ε ∆ x)ζ n+1 j,+ (3.8)
Taking into account that λ ε = O(1) uniformly in ε, there are 2 singular limits:

1. the term cos(2λ ε ∆ x) εD ε ∆ x behaves as follows as ε → 0,

cos(2λ ε ∆ x) ∆ x ε cos(2λ ε ∆ x) -i ε 2 V t +m λ ε sin(2λ ε ∆ x) → 1 ∆ x ε -2i∆ xm tan(2λ ε ∆ x) 2λ ε ∆ x .
Assuming that λ ε ∆ x → 0 as ∆ x → 0, which is immediate thanks to the uniform bound on λ ε , one sees that this term is endowed with the correct behavior because

tan(2λ ε ∆ x) 2λ ε ∆ x → 1, ∆ x → 0.
2. the term εV t sin(2λ ε ∆ x)

ελ ε D ε ∆ x
behaves as follows as ε → 0,

V t sin(2λ ε ∆ x) ∆ xλ ε cos(2λ ε ∆ x) -i∆ x εV t + m ε sin(2λ ε ∆ x) → 2V t 2∆ xλ ε cos(2λ ε ∆ x) sin(2λ ε ∆ x) -2i ∆ xm ε .
Now, still assuming that λ ε ∆ x → 0 as ∆ x → 0, one sees that 2∆ xλ ε tan(2λ ε ∆ x) → 1, as expected. But the imaginary term 2i m∆ x ε doesn't vanish unless m∆ x ≤ ε 1+α , α > 0, which is a stringent restriction on the computational grid. Surprisingly, in order to correctly recover the potential term in the non-relativistic limit, one must impose a restriction on ∆ x/ε. This phenomenon seems to be new: it doesn't manifest itself in the diffusive limits previously investigated: see [START_REF] Gosse | An asymptotic-preserving well-balanced scheme for the hyperbolic heat equations[END_REF][START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws[END_REF].

Fortunately, the issue raised by the imaginary term displaying an incorrect behavior in the limit ε → 0 can be easily circumvented by substituting D ε by ℜ(D ε ) in the denominator, its real part (which is endowed with the correct behavior as ε → 0 in absence of restrictions on ∆ x). Another fix may be the following one:

D ε Dε = cos(2λ ε ∆ x) -i εV t + max(1, 1 ε )m λ ε sin(2λ ε ∆ x).
A simple asymptotic-preserving scheme in the non-relativistic limit of (3.2) is deduced from (3.8) by,

                     ζ n+1 j,+ = ζ n j,+ -∆t ε∆ x (ζ n+1 j,+ -ζ n+1 j,-) + ∆t εD ε ∆ x ζ n+1 j-1,+ -cos(2λ ε ∆ x)ζ n+1 j,- + i∆t εℜ(D ε )∆ x εV t λ ε sin(2λ ε ∆ x)ζ n+1 j,- ζ n+1 j,-= ζ n j,-+ ∆t ε∆ x (ζ n+1 j,+ -ζ n+1 j,-) + ∆t εD ε ∆ x ζ n+1 j+1,--cos(2λ ε ∆ x)ζ n+1 j,+ + i∆t εℜ(D ε )∆ x εV t λ ε sin(2λ ε ∆ x)ζ n+1 j,+ (3.9)
thanks to the fact that ℜ(D ε ) retains only the correct (real) part of the denominator.

Remark 3.2 For the L 2 -norm, a mid-point time-integrator for (3.9) is an appealing choice thanks to the Cayley theorem. If the time-potential is not constant in x, the scheme (3.9) involves both values V t (x j± 1 2

) acting on each ζ n j,∓ , respectively. Numerical computations involving a Crank-Nicolson time-integrator like (3.5) are in §4.3.

The formal Schrödinger limit is obtained by adding both equations of (3.9),

ψn+1 j = ψn j - ∆t εD ε ∆ x ψn+1 j-1 -2 cos(2λ ε ∆ x) ψn+1 j + ψn+1 j-1 + i∆t sin(2λ ε ∆ x) ℜ(D ε )λ ε ∆ x V t 2 ψn+1 j , D ε = cos(2λ ε ∆ x) -i εV t + m/ε λ ε sin(2λ ε ∆ x).
under the innocuous conditions ε → 0, λ ε ∆ x → 0 (and after re-substituting V t → V t 2 ). Such a first-order "implicit Euler" discretization, which appears to be an original Schrödinger scheme containing tiny relativistic features, dissipates the L 2 norm.

Proposition 1 Assume the time-potential V t ∈ L ∞ (R) is a real smooth function, then the scheme (3.10) is asymptotic-preserving in the non-relativistic limit ε → 0.

Proof: With the former computations at hand, there are 2 cases to look at:

-Assume first that the time-potential V t (x) ≥ 0, the main novelty is that both λ ε and D ε depend on space and one must consider:

     λ ε (x j-1 2 ) = (εV t (x j-1 2 ) + m/ε) 2 -(m/ε) 2 ∈ R, D ε (x j-1 2 ) = cos 2λ ε (x j-1 2 )∆ x -i εV t (x j-1 2 )+m/ε λ ε (x j-1 2 ) sin 2λ ε (x j-1 2 )∆ x .
The scheme (3.9) is amended as follows:

                               ζ n+1 j,+ = ζ n j,+ -∆t ε∆ x (ζ n+1 j,+ -ζ n+1 j,-) + ∆t εD ε (x j-1 2 )∆ x ζ n+1 j-1,+ -cos(2λ ε (x j-1 2 )∆ x)ζ n+1 j,- + i∆t εℜ(D ε (x j-1 2 ))∆ x εV t (x j-1 2 ) λ ε (x j-1 2 ) sin(2λ ε (x j-1 2 )∆ x)ζ n+1 j,- ζ n+1 j,-= ζ n j,-+ ∆t ε∆ x (ζ n+1 j,+ -ζ n+1 j,-) + ∆t εD ε (x j+ 1 2 )∆ x ζ n+1 j+1,--cos(2λ ε (x j+ 1 2 )∆ x)ζ n+1 j,+ + i∆t εℜ(D ε (x j+ 1 2 ))∆ x εV t (x j+ 1 2 ) λ ε (x j+ 1 2 ) sin(2λ ε (x j+ 1 2 )∆ x)ζ n+1 j,+ (3.10)
Thanks to the smoothness of x → V t (x), the interface coefficients are still consistent because V t (x j∓ 1 2 ) -V t (x j ) → 0 strongly as ∆ x → 0.

-In case V t (x) has no definite sign, it can happen that λ ε (x j-1 2

) becomes purely imaginary for some j ∈ Z. But this doesn't create any substantial issue because all the (circular) trigonometric relations which were used formerly still hold for (hyperbolic) trigonometric functions: for all x ∈ R,

cos(ix) = exp(i 2 x) + exp(-i 2 x) 2 = cosh(x), sin(ix) = exp(i 2 x) -exp(-i 2 x) 2i = i sinh(x),
thus one recovers cos 2 (ix) + sin 2 (ix) = 1, together with all the limits as x → 0.

Numerical results in various settings

A test-case for linear problems

Here we simply aim at displaying the behavior of the explicit scheme (2.7) on a standard test-case for (2.1) involving both a pseudo-scalar and a time potential. More precisely, the pseudo-scalar term V s (x) = αx 2 2 has a confinement effect, free from the Klein paradox. The time potential corresponds to a wall modeled by means of an indicator function, V t (x) = 3 2 χ 1.25<x<1.75 . Initial data for ζ ± (t = 0, x) read

1 σ √ 2π exp - (x + 1.5) 2 2σ 2 , σ = 0.125.
The computational domain is x ∈ (-3, 3) and N = 2 8 points are used to grid it. On Fig. 4.1, results are shown at different times with a Courant number ν = 1, allowing for an exact (up to machine accuracy) conservation of the L 2 norm. We show the moduli added to the pseudo-scalar potential, V s (x) + |ζ + | 2 + |ζ -| 2 in order to emphasize its effects on the wave packet's dynamics. Zitterbewegung is observed especially when the movement is interrupted by the wall, and a small tunneling effect takes place.

Formal extension to nonlinear problems

Well-balanced upwind scheme

Following some examples displayed in [START_REF] Bournaveas | Theory and numerical approximations for a nonlinear Dirac system[END_REF], we intend to test the scheme (2.7) against exact solutions of a nonlinear Dirac equation (compare with (2.1), no time-potential):

∂ t Ψ + 0 1 1 0 ∂ x Ψ + i 1 + |ψ -| 2 -|ψ + | 2 -1 + |ψ + | 2 -|ψ -| 2 •Ψ = 0. ( 4.1)
According to [START_REF] Alvarez | Linearized Crank-Nicolson scheme for nonlinear Dirac equations[END_REF][START_REF] Alvarez | The numerical study of a nonlinear one-dimensional Dirac equation[END_REF][START_REF] Bournaveas | Theory and numerical approximations for a nonlinear Dirac system[END_REF], this equation admits a family of exact solutions given by:

ψ + (t, x) = exp(-iΛt)A(x), ψ -(t, x) = i exp(-iΛt)B(x),
where

Λ ∈ (-1, 1), µ = √ 1 -Λ 2 , µ A = µ 2(1 + Λ ) µ B = µ 2(1 -Λ ), and
A(x) = µ A cosh(µx) 1 + Λ cosh(2µx) , B(x) = µ B sinh(µx) 1 + Λ cosh(2µx)
.

In order to tackle (4.1) by means of (2.7), it suffices therefore to cancel the timepotential V t and to include the nonlinearity |ψ -| 2 -|ψ + | 2 , "frozen at each time t n ", in the pseudo-scalar potential V s . This yields a "nonlinearly varying mass" instead of the linear space-varying mass m(x) showing up in (2.1). Numerical results for the modulus of Ψ with N = 2 9 points griding the interval x ∈ (-32, 32) and a Courant number equal to 1, ∆t = ∆ x, are displayed on Fig. 4.2, left. The time T = 2π Λ is chosen so that the exact solutions are equal to the initial data. Hence, despite the aforementioned derivations don't rigorously extend to nonlinear cases, it is found numerically that the scheme (2.7) behaves rather correctly on these exact solutions, the L 2 norm is conserved too, up to machine accuracy (see Fig. 4.2, right).

Testing the well-balanced (WB) MUSCL extension

As said before, the L 2 -preservation property holds essentially for ν = 1, that is, in the absence of numerical diffusion (Lemma 2.2). One may wonder whether the MUSCL reconstruction (2.9)-(2.10) is able to improve this situation. Of course, as explained in Remark 2.1, the L 2 norm corresponds to an entropy for the semi-linear hyperbolic system (see [START_REF] Zhang | Global strong solution to a nonlinear Dirac-type equation in one dimension[END_REF][START_REF] Zhang | Global solution to a cubic nonlinear Dirac equation in 1 + 1 dimensions[END_REF] for applications of hyperbolic techniques to Dirac systems) thus it is dissipated by finite-difference methods for the sake of stability in the presence of steep gradients. On the left of Fig. 4.3, we display the sensitivity of the time-evolution of the L 2 norm with respect to the Courant number for both the upwind (2.7) and the MUSCL (2.9)-(2.10) schemes on the semilinear problem (4.1). The improvement in terms of L 2 -norm is of the order of 15% roughly, so this type of numerical scheme is too much dissipative for Courant numbers below 0.8. On the right of Fig. 4.3, we show that, for ν = 0.8 the L 2 -error grows linearly in time while displaying some oscillations. Moreover, the decay of these errors when ∆ x passes from 2 -7 to 2 -10 appears to be linear, too. The parameters were kept identical to those of Fig. 4.2.

4.3 Non-relativistic limit and Asymptotic-Preserving (AP)

We confront numerically the scheme (3.9) against a direct Schrödinger computation: more precisely, a Crank-Nicolson time-integrator is applied both to (3.2), with the same numerical fluxes leading to (3.9), and to the usual non-relativistic equation,

i∂ t ψ + 1 2m ∂ xx ψ = V t (x)ψ, V t (x) = x 2 2 , (4.2)
in the domain x ∈ (-5, 5) with N = 2 9 grid points. Initial data correspond to a Gaussian pulse centered in x = - 5 4 and with wave number k = 1. Numerical densities of presence are nearly identical as soon as ε < 0.001 and they roughly match the ones generated by the direct Schrödinger computation (see Fig. 4.4, bottom line). Oppositely, when the Dirac equation is in relativistic regime ε ≃ 1, noticeable differences show up as, for small time, the initial pulse splits into 2 components (Fig. 4.4, top, left), and despite the potential's effect, the numerical solutions are quite different at later time (Fig. 4.4, top, right). This was to be expected as the solutions for Dirac and Schrödinger equations have no reason to look like each other in relativistic regime. Following Remark 3.2, a Crank-Nicolson approach of the type (3.5) was implemented for this benchmark in order to show its numerical stability even in the presence of a singular limit ε → 0; practically, it seems that simulating the "relaxing system" (3.2) is less demanding in terms of numerical boundary conditions. We chose to implement Neumann boundary conditions, even if it is of little importance because the quantum particle is unable to reach the borders in presence of the steep parabolic potential (this situation is similar to the one of Fig. 4.1).

Conclusion

The classical defect of finite-difference methods in the context of quantum mechanics is their numerical viscosity which generally perturbs the L 2 -norm of the approximate solution (one remedy is the setup of second-order implicit time-integrators [START_REF] Hairer | Geometric numerical integration. Structure-preserving algorithms for ordinary differential equations[END_REF][START_REF] Lubich | Integrators for Quantum Dynamics: A Numerical Analyst's Brief Review, in Quantum Simulations of Complex Many-Body Systems: From Theory to Algorithms[END_REF][START_REF] Wessels | Discretizing the one-dimensional Dirac equation[END_REF]). Here we showed that, at least in the relativistic scaling, the well-balanced approach yields a numerical scheme (2.7) on an inhomogeneous linear problem which preserves exactly the L 2 -norm while being completely explicit in time. When it comes to the (singular) non-relativistic limit of this scheme, there are 2 different situations:

1. for free particles (no external time-potential), the well-balanced scheme (2.7) can be reformulated into an asymptotic-preserving one (3.8) by mimicking former computations designed in the context of discrete kinetic models [START_REF] Gosse | An asymptotic-preserving well-balanced scheme for the hyperbolic heat equations[END_REF], 2. for particles submitted to an external time-potential, the coupling between different source terms leads to a "strange imaginary term" in D ε which produces an incorrect behavior on the potential only (the Laplacian term is always correct). However, it is easy to fix this issue by taking the real part of this D ε coefficient.

A very different strategy for quantum mechanics equations, better suited for L 2 -norm preservation, is the split-step methods [START_REF] De Frutos | Split-Step Spectral Schemes for Nonlinear Dirac Systems[END_REF][START_REF] Huang | A time-splitting spectral scheme for the Maxwell-Dirac system[END_REF][START_REF] Markowich | Numerical approximation of quadratic observables of Schrödinger-type equations in the semiclassical limit[END_REF]; these algorithms too can furnish asymptotic-preserving strategies in the non-relativistic limit, see Remark 4.3 in [26]. This scheme has all the strengths and weaknesses which characterize the Fourierbased time-splitting schemes:

-Spectral accuracy and L 2 -norm conservation on periodic smooth solutions; -Possible issues to reproduce Jost "distorted plane wave" solutions (the ones which constitute the continuous spectrum of the stationary operator) in a scattering environment [START_REF] Duchêne | Wave operator bounds for one-dimensional Schrdinger operators with singular potentials and applications[END_REF]. Reflection/transmission coefficients may result flawed unless a fine grid is used with potentials V (x) being either oscillating (see [START_REF] Jin | Mathematical and computational methods for semiclassical Schrodinger equations[END_REF], page 183) or sharply-varying (thus inducing high frequencies).

Multi-dimensional problems are usually handled by dimensional splitting, but it is delicate matters in case one wants to reproduce accurately wave interactions in a context of weak solutions. A first step may be to consider the simpler case of socalled "massless Dirac equations" for which the 4-component spinor reduces to a 2-component one, see for instance [START_REF] Noelle | Hyperbolic systems of conservation laws, the Weyl equation, and multidimensional upwinding[END_REF].

A Discrete H 1 (R)-bounds and strong convergence of free particles Lemma A.1 Assume c∆t ≤ ∆ x, x → V t (x),V s (x) being Lipschitz and ζ ± (t = 0,•) ∈ H 1 (R), then numerical approximations generated by (2.7 satisfy:

∀n ∈ N, ζ ∆x (t n ,•) L 2 (R) ≤ ζ ∆x (0,•) L 2 (R) . If moreover, ∂ x V t = ∂ x V s ≡ 0, then: ∀n ∈ N, ∂ x ζ ∆x (t n ,•) L 2 (R) ≤ ∂ x ζ ∆x (0,•) L 2 (R) .
Proof: Hereafter we denote

|ζ ± | = |ζ + | 2 + |ζ -| 2 the Euclidian vectorial norm and recall that |S j+ 1 2 ζ ± |.
-As a consequence of the compact writing (2.8), one has

ζ n+1 j,+ ζ n+1 j-1,- = (1 -ν) ζ n j,+ ζ n j-1,- + νS j-1 2 ζ n j-1,+ ζ n j,-
, and deduces easily an L 2 estimate for ν ≤ 1 because each scattering matrix S j-1 2 , j ∈ Z, is unitary.

  ∑ j∈Z ∆ x ζ n+1 j,+ ζ n+1 j-1,- 2   1 2 ≤ (1 -ν) ∑ j∈Z ∆ x ζ n j,+ ζ n j-1,- 2 1 2 + ν ∑ j∈Z ∆ x S j-1 2 ζ n j-1,+ ζ n j,- 2 1 2 ≤ (1 -ν) ∑ j∈Z ∆ x ζ n j,+ ζ n j-1,- 2 1 2 + ν ∑ j∈Z ∆ x ζ n j-1,+ ζ n j,- 2 1 2 ≤ ∑ j∈Z ∆ x ζ n j,+ ζ n j-1,- 2 1 2 ≤ ∑ j∈Z ∆ x ζ n=0 j,+ ζ n=0 j,- 2 1 2
.

The first step is just Minkowski's inequality and linearity of the L 2 norm for 0 ≤ ν ≤ 1. The second one is the S-matrix being unitary. The third is shifting sum's indexes in order to produce

∑ j∈Z ζ n j-1,+ ζ n j,- 2 = ∑ j∈Z ζ n j-1,+ | 2 + |ζ n j,- 2 = ∑ j∈Z ζ n j,+ | 2 + |ζ n j-1,- 2 = ∑ j∈Z ζ n j,+ ζ n j-1,- 2 .
This reordering of the sums is specific to the L 2 norm; therefore one needs both the S-matrix being unitary and the ability to rearrange the summation in order to derive the final estimate.

-When the scattering matrix is space-dependent, this doesn't automatically yield an estimate on the divided-differences: indeed, one sees on Fig. 4.1 that the presence of external potentials can increase the Lipschitz constant of the numerical approximations. One begins by "triangulating" as follows: .

We denoted by Lip x (S) the Lipschitz constant in x of the S-matrix, which depends on the Lipschitz constant of both V t (x) (the time potential) and V s (x) (the pseudo-scalar term). One cannot deduce a discrete H 1 estimate because of that last term. However, if V t and V s are constant, Lip x (S) ≡ 0 and: Lemma A.1 states that the L 2 -norm (a particular entropy) is dissipated in time as soon as the Courant number ν < 1 by the upwind finite-differencing process. This is generally seen as a drawback for such a method; however, in certain circumstances and for short-time computations, dissipation can be interesting when an initial data oscillating at a frequency close to the grid's cutoff frequency π/∆ x is prescribed. These data usually induce wave-packets endowed with flawed group velocities (see [START_REF] Micu | Uniform boundary controllability of a semi-discrete 1-D wave equation with vanishing viscosity[END_REF][START_REF] Zuazua | Propagation, observation, and control of waves approximated by finite difference methods[END_REF] and references therein) of the order of ∆ x and energy conservation implies that they never disappear. Instead, running a numerical scheme containing a little bit of numerical viscosity allows to progressively smoothen high frequencies and restore wave-packets with group velocities compatible with the computational grid.

Fig. 2 . 1

 21 Fig. 2.1 Schematic view of the scheme (2.4) with Riemann problems for (2.3) and standing waves.

2 Fig. 2 . 2

 222 Fig. 2.2 Schematic view of the MUSCL reconstruction and scattering matrices

Fig. 3 . 1

 31 Fig. 3.1 Position densities generated by (3.6)-(3.7) at time t = 0.5 for ε = 1,0.1,0.01 (left) and timegrowth of the L 2 difference between them and Leapfrog approximation for the free Schödinger equation.

Fig. 4 . 1

 41 Fig. 4.1 Snapshots of (2.7) at t = 0,1.54,2.08,3.54 of a pulse on a parabola hitting a wall with c = 1.

Fig. 4 . 2

 42 Fig. 4.2 Solution of (4.1) for T = 2π Λ , Λ = 3 4 (left), L 2 -norm change (right).

2 L2Fig. 4 . 3

 243 Fig. 4.3Relative dissipation of L 2 -norm for ν = 0.9,0.65,0.4 and ∆ x = 2 -8 (left), growth of L 2 errors for 4 decreasing values of ∆ x = 2 -7 ,2 -8 ,2 -9 ,2 -10 but ν = 0.8 kept constant (right), both functions of time.

Fig. 4 . 4

 44 Fig. 4.4 Comparison of Crank-Nicolson results on both (3.2) and Schr ödinger equation (4.2) at T = 1 (left), T = 3 (right) and for ε = 1 (top), ε = 0.01 (bottom).

  + 1)∆t • Lip x (S) • ∑ j∈Z

Lemma A. 2 2 . 2 = 2 |

 2222 Under the hypotheses of Lemma A.1, assume both V t ,V s are constant, then:∀n ∈ N, ζ ∆x (t n+1 ,•) -ζ ∆x (t n ,•) L 2 (R) ≤ c∆t ∂ x ζ ∆x (0,•) L 2 (R) + Lip ∆x (S) ζ ∆x (0,•) L 2 (R) .Proof: It consists in evaluating the size of the time-variation and taking advantage of previous H 1 estimates:Now, from the expression (2.5), one sees that for any j ∈ Z, if ∆ x = 0, S j-1 Id, the identity matrix.Hence |Id -S j-1 ≤ ∆ x • Lip ∆x (S) and this is enough to ensure time-equicontinuity thanks to (A.1). ⊓ ⊔ Lemmas A.1 and A.2 allow to conclude that strong convergence holds in L 2 loc (R + * × R) (and almost everywhere) as ∆ x → 0, c∆t ≤ ∆ x for free particles of constant mass m without external potential.

Acknowledgements

The author thanks Dr. S. Malaguti (Ferrara) for having suggested this problem.


Quentin Ehlinger

Laurent Joly
email: laurent.joly@univ-lyon1.fr

Olivier Pierre-Louis

Giant slip at liquid-liquid interfaces using hydrophobic ball bearings

Keywords: numbers: 68.05.-n, 83.50.Lh, 68.03.Cd, 47.61.-k

Liquid-gas-liquid interfaces stabilized by hydrophobic beads behave as ball bearings under shear and exhibit giant slip. Using a scaling analysis and Molecular Dynamics simulations we predict that, when the contact angle θ between the beads and the liquid is large, the slip length diverges as Rρ -1 (πθ) -3 where R is the bead radius, and ρ is the bead density.

Starting with the seminal work of Navier in 1823 [1], the study of interfacial slip at liquid-solid interfaces has a long history. Slip is usually not observed at macroscopic scales, but the recent downsizing of hydrodynamic flows in micro and nanofluidic devices [2,3] has paved the way for measurements at small length scales. These measurements revealed that slip occurs at the nanometer scale on flat substrates [4][5][6][7][8]. However, nature has produced surfaces such as plant leaves, with complex topographic structures (bumps, hairs, etc.) on the top of which water drops can be deposited without collapse in the so-called Cassie-Baxter (or fakir) state [9][10][START_REF] Shirtcliffe | Newton[END_REF]. Man-made devices based on this principle led to an increase of slip by orders of magnitude, using e.g. substrates with nanopillar arrays [12][13][14][15][16][17]. One may therefore wonder whether a similar strategy could be used to increase slip in liquid-liquid interfaces. Indeed, a few studies considering the possibility of slip at the bare interface between two simple liquids have shown that slip, if present, was negligible above the molecular scale [18,19]. Surprisingly, virtually no work has been devoted to optimizing slip at liquid/liquid interfaces. A large liquid-liquid slip would open new possibilities in micro and nanofluidics, permitting different liquids in contact to flow independently with reduced interfacial friction. In the following, we show that it is possible to achieve giant slip between two liquids using a bed of hydrophobic beads, as depicted in Fig. 1(a).

Following the principle of usual ball bearings, we suggest to use an interface where hydrophobic beads maintain a gas layer between two bulk liquids. The stability of this interface is ensured by wetting forces, instead of being controlled by the contact between the beads and the solid substrate in usual ball bearings. Such a metastable configuration reminds of pillar-based superhydrophobic surfaces, used to amplify liquid-solid slippage [12][13][14][15][16][17]. In both cases the liquid rests mainly on a gas layer, leading to a low level of friction. However there are two important differences between these systems: to the advantage of hydrophobic ball bearings, beads are able to roll, thereby reducing friction at the liquid/bead interface. However beads will always penetrate inside the liquid, thereby inducing viscous dissipation, and consequently decreasing slippage. As seen in Fig. 1(a,b), the penetration depth is directly controlled by the wetting angle θ of the liquid at the bead surface. As a consequence, θ is expected to have a strong influence on the efficiency of liquid/liquid bearings. In the following, we start by quantifying analytically the influence of the wetting angle on liquid/liquid slip in this system. We then confirm the obtained scaling law by means of Molecular Dynamics (MD) simulations. Finally, we discuss the expected amplitude of slip in experimental systems.

i i y x -v v y x -V V L L θ R h V -V (b) (c)
In order to discuss the behavior of liquid-liquid bearings under shear, we shall consider a simplified geometry with a periodic array of beads. The unit cell contains a single bead, as sketched in Fig. 1(b). Shear is forced by a velocity difference 2V between two parallel walls separated by a distance 2L y . The x and y axes are indicated in the schematic, and the z axis is orthogonal to the plane of the schematic. The system is assumed to be periodic along the x and z axes. At liquid-solid interfaces, the Navier partial slip boundary condition (BC) [1,6], ∆v = b γ, relates the velocity jump at the interface ∆v to the far-field shear rate γ inside the liquid, where b is the slip length. The slip length between two liquids with the same viscosity may be defined using the same equation without ambiguity (see Supplementary material). Therefore, one can use the velocity v i extrapolated at the median plane passing though the center of the bead as in Fig. 1

(c), leading to b = 2v i γ , (1)
where γ = (V -v i )/L y is the macroscopic shear rate in the liquid. Focusing on small scales hydrodynamics, we investigate the limit of small Reynolds numbers and negligible gravity effects. We consider a steady state. In the referential of the bead center of mass, the solid-liquid and liquid-gas interfaces have null normal velocity (however, they can exhibit a non-vanishing tangential velocity). In such a fixed geometry, a minimum dissipation principle is available for Stokes flow (i.e. the limit of vanishing Reynolds numbers) with no-slip boundary condition as discussed historically by Helmholtz [20], Rayleigh [21], and Korteweg [22]. This principle is easily extended to a variable slip length along the boundary [23], and therefore applies to our system, where the interface is a sum of portions of liquid-gas interface with infinite slip and liquid-solid interface with zero or partial slip. Since the walls are only an artificial setup to enforce shear, we consider a no-slip BC at the wall-liquid interface in the following. However, at the physically relevant liquid-bead interface, we shall discuss both cases of perfect slip and vanishingly small slip.

Let us start with the case of a perfect slip. The steady state corresponds to the minimization of the total viscous dissipation:

Q = η 2 d 3 r [∂ µ v ν + ∂ ν v µ ] 2 , (2)
where µ, ν represent x, y, or z, and the Einstein summation convention is assumed. Our strategy is to evaluate Q as a function of the interface velocity v i , and to select the velocity v i * for which Q is minimized. We therefore need to analyze of the flow field v = (v x , v y , v z ) present in the system. The first contribution to v is due to the direct shear flow v D imposed by the boundary conditions. In the upper part of the system, with y > 0, we have

v D = v i + y γ x. ( 3
)
The second contribution is the backflow v B caused by the response of the system to the direct flow. We assume that the liquid-gas interface remains flat, so that the backflow is mainly caused by the presence of the bead. Its main component is along y:

v B = ∂ x y s v Dx ŷ, (4)
where y s (x, z) is the position of the liquid boundary (liquid-solid plus liquid-gas interfaces). Since θ → π, we have ∂ xx y s ≈ -1/R and y s ≈ R above the bead, leading to

∂ x v By ≈ - 1 R v i + R γ . (5)
Using Eq. ( 2), the total dissipation now reads:

Q = η 2 d 3 r γ - 1 R v i + R γ 2 . ( 6
)
The backflow penetrates the liquid in a domain of volume ℓ 3 , where ℓ ≈ 2R(π -θ) is the lateral extent of the solidliquid contact region. Hence, the backflow terms in the integrand of Eq. ( 6) should be integrated in this volume only. Using the separation of length scales ℓ ≪ R ≪ L y , we then find

Q = η L x L z L y (V -v i ) 2 + 8ηαR(π -θ) 3 (v i ) 2 , (7)
where α is a dimensionless prefactor accounting for the system geometry. We may now find the slip velocity v i * which minimizes the dissipation from the equation

∂ v i Q 2D | v i =v i * = 0. Finally, using v i * in Eq. (1), we ob- tain b = 1 4ρRα(π -θ) 3 , (8)
where ρ = 1/(L x L z ) is the bead density at the interface. The no-slip limit is more subtle, due to the expected divergence of viscous dissipation at the triple line [24,25], as the intrinsic liquid/bead slip length b 0 vanishes. This contribution enters into play as a divergence of ∂ xx y s at the triple line, replacing 1/R by (π -θ)/b 0 in Eq. ( 5). Assuming that the backflow extends to a distance b 0 along the triple-line perimeter πℓ then suggests a dissipation 3 as in Eq. (7). A more precise analysis based on the expansion for θ → π of the exact two-dimensional solution of Stokes flow in a wedge [25] confirms this scaling, but also indicates the presence of a logarithmic correction. Indeed, as θ → π the dissipation per unit length of triple line reads ∼ ηv i2 (π -θ) 2 ln(ℓ 0 /b 0 ), where ℓ 0 is a macroscopic cutoff. Assuming ℓ 0 ∼ ℓ, we obtain a dissipation

∼ ηv i2 πℓb 2 0 [(π -θ)/b 0 ] 2 ∼ ηv i2 R(π -θ)
∼ ηv i2 R(π -θ) 3 ln[(π -θ)R/b 0 ], from which we compute b ∼ 1/{R(π -θ) 3 ln[(π -θ)R/b 0 ]}.
We therefore expect the scaling behavior in Eq. ( 8) to be weakly affected (i.e. within logarithmic corrections) by the boundary condition between the liquid and the bead. Note that our analysis relies on the assumption that θ does not depend on the triple-line velocity. This should be achieved for low capillary numbers [24]. In this limit, θ does not depend on η, and thus Eq.( 8) indicates that b should be independent of η.

In order to test the predictions of Eq. ( 8), we performed molecular dynamics (MD) simulations, using the LAMMPS package [26]. We consider once again a periodic array of beads, the unit cell of which is presented in Fig. 2. All atoms interacted via a Lennard-Jones (LJ) pair potential: 6 ], where i and j are the types of interacting atoms, and r and ε ij are respectively their distance and interaction energy. The zero-potential distance σ is the same for all atom types. In addition, we chose the fluid-fluid interaction energy as the reference energy ε f f = ε. In the following, we will use LJ units of energy ε, distance σ, and time τ = σ m/ε, where m is the atomic mass. In order to increase the liquid/vapor surface tension and to decrease the vapor density as compared to a simple LJ liquid, we used a liquid of LJ dimers [28], i.e. pairs of LJ atoms bonded by a harmonic potential. We used a dissipative particle dynamics (DPD) thermostat to maintain the liquid atoms at a temperature k B T = 0.73ε while preserving hydrodynamics [29]. The bead was modeled as a rigid spherical shell of radius R = 7.07σ (with a thickness larger than the LJ cutoff 2.5σ), cut inside a fcc crystal with density 1/σ 3 . The confining walls were made of three (010) layers of frozen atoms on a fcc lattice, with the same density 1/σ 3 . The contact angle is controlled by the fluid-solid interaction energy in the LJ potential. We used ε f w = 0.7ε between fluid and wall atoms, so as to ensure a small contact angle, and consequently a no-slip BC [30]. In contrast, the fluid-bead interaction energy ε f b was tuned in order to change the contact angle θ. The latter was then computed from in situ measurements of the vapor gap width. Finally, the distance between the top of the bead and the walls was ≈ 2R in all simulations. This was large enough to remove in practice the influence of the confining walls, and to allow for a linear average velocity profile in the liquid which enables the measurement of the shear rate γ. After a period of equilibration, the system was sheared by imposing opposite wall velocities ±V along x. All along the simulation, both walls were used as pistons to impose a vanishing pressure inside the system.

V ij = 4ε ij [(σ/r) 12 -(σ/r)
To measure b, we first used the kinematic definition, Eq. ( 1). In order to extract ∆v and γ, we fitted the numerical velocity profiles in the linear regions between the bead and the walls, as shown in Fig. 2. However, for large contact angles, the shear rate γ becomes increasingly small, and its measurement in the fluid is inaccurate. We therefore resort to a different definition of the slip length based on a dynamic interpretation [31]. On the one hand, the bulk viscous shear stress is η γ, where η is the liquid viscosity. On the other hand, the interface shear stress equals the interfacial friction force f per unit area, which is proportional to the velocity jump: f = λ∆v. This latter relation defines the interfacial friction coefficient λ. Equating the bulk and interface shear stresses, one recovers the partial slip BC with b = η/λ. Therefore, the slip length can also be computed as the ratio between the bulk liquid viscosity and the interfacial friction coefficient. The liquid viscosity η ≈ 5.52 LJ units, was determined under the same thermodynamic conditions in independent shear flow simulations. The friction coefficient λ was computed as the ratio of the measured shear force per unit area f and the measured velocity jump ∆v. For the analysis of the non-trivial θdependence of b, we used a fixed box size L x = L z = 4R, corresponding to ρ = 1/(16R 2 ). In Fig. 3, the resulting slip length is plotted as a function of θ. Using the dynamic method, all data points collapse on a single master curve b ∝ (π -θ) -3 , as predicted by Eq. ( 8). The fit to Eq. ( 8) provides a value for the dimensionless prefactor α = 0.31±0.02. The results of the kinematic method are similar to those of the dynamic method. However, as anticipated above, the kinematic method is less precise at small shear rates, leading to larger scatter of the data. The comparison between the results of the two methods is provided in Supplementary Material. We have also checked that b ∼ 1/ρ, as predicted by Eq. ( 8). Deviations from this latter scaling are observed when hydrodynamic interactions come into play at high bead-density and when π -θ is not small. Details are discussed in Supplementary Material.

A recurrent issue in MD simulations is the need to use a very large forcing in order to extract the system response out of thermal noise. Such large driving forces could lead to a nonlinear response, questioning the description of the interface behavior by means of a slip length. Our shear rates were ∼ 10 7 to 10 9 Hz, which is very large as compared to experimentally relevant shear rates, typically smaller than 10 5 Hz. As shown in Fig. 3, we varied V over almost one order of magnitude and checked that the measured slip lengths remained unaffected. We may therefore conclude that the system response is linear, and can be safely extrapolated to the much lower experimental shear rates. In addition, it is clear from the adequacy between Eq.( 8) and MD simulations that thermal fluctuations play a negligible role in the emergent macroscopic behavior, at least in the explored parameter range.

Finally, note that we do not control the slip at the liquid-bead interface in our MD simulations. It is known that both the contact angle and the curvature of the solid surface may have a strong impact on the liquid-bead slip [30,32]. However, as discussed above, the variation of the liquid-bead slip length b 0 is not crucial for our main result Eq. ( 8), and should lead to minor logarithmic corrections, the discussion of which are beyond the scope of the present Letter.

Interfaces similar to that of Fig. 1(a) have actually already been observed experimentally using liquid marbles floating on liquid substrates [33,34]. However, to our knowledge no analysis of slip has been performed for these interfaces. An important way to reduce slip in these systems is to decrease the bead density, which directly enters into the expression of b in Eq.( 8). In addition, lowering the density may help to prevent contact and friction between the beads, which is expected have negative consequences on slip. One possible way to control the density is to use electrostatic repulsion between the beads along the same lines as the in formation of colloidal Wigner crystals at water-oil interfaces [35]. Let us assume that one could build such a system, with typical densities ρ varying from 0.01R -2 to 0.1R -2 . Assuming 0.03R -2 (corresponding to a distance of about two diameters between the beads), Eq.( 8) predicts b ≈ 20R with θ = 120 • , and b ≈ 200R with θ = 150 • (using e.g. superhydrophobic beads). With a bead radius R = 10 µm, b would therefore be on the scale of millimeters for large contact angles (b ≈ 0.2 and 2 mm for θ = 120 and 150 • respectively).

In our opinion, a strong limitation of this giant slip is the stability of the interface under pressure variations. Two processes are involved. First, the Laplace-pressureinduced curvature of the liquid-gas interface imposes an additional backflow, increasing dissipation. This effect is similar to the effect of pressure in pillar-based superhydrophobic surfaces [36,37]. Defining the curvature κ of the liquid-gas interface, and assuming that the related backflow extends in the liquid up to a distance similar to the inter-bead distance, we obtain an additional dissipation ∼ ηρ -3/2 κ 2 (v i ) 2 in Eq. (7). Recalling the Laplace law κ = ∆P/2γ, where ∆P is the pressure variation and γ is the surface tension, we obtain a correction of the slip length b -1 = b -1 | ∆P =0 + α ′ ρ -1/2 (∆P ) 2 /4γ 2 with α ′ a constant of order one, showing that b decreases when increasing |∆P |. The correction is negligible when |∆P | < 2γρ 1/4 b -1/2 , where the zero-pressure slip length b is given by Eq. ( 8). This condition becomes stringent when b is large, and for θ = 150 • with R = 10 µm, one obtains |∆P | < 8 mbar. The second limitation related to the pressure is the possibility of contact between the two opposite liquid-gas interfaces, leading to the collapse of the system. Balancing the typical interface displacement κρ -1 with the typical average separation ∼ R, the collapse will be avoided when ∆P < 2γRρ. Using R = 10 µm and ρ = 0.03R -2 , one finds again P < 4 mbar. The stability can be enhanced, at the cost of a smaller slip length, by using higher densities and smaller beads. For instance, with R = 0.1 µm and ρ = 0.1R -2 , one finds P < 1 bar, but b ≈ 6 µm for θ = 150 • .

In conclusion, giant slip can be obtained between two liquids using hydrophobic ball bearings. We hope that this work will motivate experimental investigations of the dynamic properties of such interfaces. Liquid-liquid bearings open new pathways for micro-and nano-fluidics research. One major direction could be to build fluidic devices without walls, where different liquids could flow with regard to each other while maintaining an extremely low level of friction, and preventing mixing by diffusion between the different channels.

OPL wishes to thank J. Yeomans and H. Kusumaatmaja for discussions.

FIG. 1 .

 1 FIG. 1. (color online) Schematics. (a) Using hydrophobic beads as a ball bearing to enhance liquid/liquid slip. (b) Unit cell of the periodic system analyzed in the text. (c) Expected macroscopic flow in the interface region.

FIG. 2 .

 2 FIG.2. (color online) Snapshot of a typical molecular dynamics system (created using VMD[27]), with the associated steady-state velocity profile on the right.

FIG. 3 .

 3 FIG. 3. (color online) Slip length b as a function of the measured wetting angle θ, using the dynamic definition of the slip length. Data set keys indicate the shear velocity V , in LJ units (σ/τ ). All data sets collapse on a master curve b ∝ (πθ) -3 (dashed line).