F Dhaouadi
email: dhaouadi@insa-toulouse.fr

N Favrie

S Gavrilyuk
email: sergey.gavrilyuk@univ-amu.fr

Extended Lagrangian approach for the defocusing Non-Linear Schrödinger equation

Keywords: Non-Linear Schrödinger equation, augmented Lagrangian, hyperbolic equations, Godunov type methods

We study the defocusing Non-Linear Schrödinger (NLS) equation written in hydrodynamic form through the Madelung transform. From the mathematical point of view, the hydrodynamic form can be seen as the Euler-Lagrange equations for a Lagrangian submitted to a differential constraint corresponding to the mass conservation law. The dispersive nature of NLS equation poses some major numerical challenges. The idea is to introduce a two-parameter family of extended Lagrangians, depending on a greater number of variables, whose Euler-Lagrange equations are hyperbolic and accurately approximate NLS equation in a certain limit. The corresponding hyperbolic equations are studied and solved numerically using Godunov type methods. Comparison of exact and asymptotic solutions to the NLS equation ('grey' solitons and dispersive shocks) and the corresponding numerical solutions to the extended system was performed. A very good accuracy of such a hyperbolic approximation was observed.

Introduction

In continuum mechanics, the most of dissipationless models can be derived from a constrained Lagrangian that depends only on the instantaneous values of macroscopic variables and not on their derivatives. The imposed constraints correspond to the conservation of main physical quantities (density, entropy,...). However, it is often necessary to construct a Lagrangian that also depends on gradients of macroscopic variables and (or) their material time derivatives. The mathematical nature of the corresponding Euler-Lagrange equations completely changes : usually, we shift from a 'hyperbolic' system of equations to a 'dispersive' one. The examples of dispersive models are numerous : shallow water equations (Serre-Green-Naghdi (SGN) equations [START_REF] Green | A derivation of equations for wave propagation in water of variable depth[END_REF] or their generalization accounting for capillary effects [START_REF] Dias | On the fully-nonlinear shallow-water generalized Serre equations[END_REF]), bubbly fluid flows [START_REF] Gavrilyuk | Generalized vorticity for bubbly liquid and dispersive shallow water equations[END_REF], solid mechanics for gradient elasticity [START_REF] Forest | Micromorphic approach for gradient elasticity, viscoplasticity, and damage[END_REF], non-linear optics [START_REF] Gurevich | Dissipationless shock waves in media with positive dispersion[END_REF]. For further examples, see a review article [START_REF] El | Dispersive shock waves and modulation theory[END_REF]). Only few of them are integrable, so the problem of construction of analytical or numerical solutions to the mathematical models is particularly challenging.

Currently, numerical methods used for dispersive systems of equations are extremely time-consuming. For example, when treating the SGN equations, it is necessary at each time step to invert an elliptic operator in the computational domain [START_REF] Métayer | A numerical scheme for the Green-Naghdi model[END_REF]. Also, it is not completely unclear how to impose the transparent boundary conditions for general dispersive systems allowing us to avoid 'parasitic reflection' from artificial numerical boundaries [START_REF] Besse | Discrete transparent boundary conditions for the mixed KdV-BBM equation[END_REF]. Recently, a new approach was developed for the SGN equations [START_REF] Favrie | A rapid numerical method for solving Serre-Green-Naghdi equations describing long free surface gravity waves[END_REF] allowing us to avoid these problems. This approach is based on the formulation of an extended Lagrangian that allows us to transform the dispersive equations into a set of hyperbolic equations with relaxation terms. The obtained equations become local and there is no more need to invert non-local operators. In particular, this makes the computation time sharply reduced. Also, for hyperbolic systems of equations the question about transparent boundary conditions can be answered as fully as possible at least in linear approximation.

We extend here the approach of [START_REF] Favrie | A rapid numerical method for solving Serre-Green-Naghdi equations describing long free surface gravity waves[END_REF] to the case where the dispersion is described by the gradients of unknowns. Such models appear, in particular, in the study of capillary fluids [START_REF]La capillarité interne[END_REF], [START_REF]Capillarité interne en mécanique des millieux continus[END_REF], [START_REF] Casal | A representation of liquid-vapor interfaces by using fluids of second grade[END_REF], [START_REF] Eglit | A generalization of the model of an ideal compressible fluid[END_REF], [START_REF] Gavrilyuk | Media with equations of state that depend on derivatives[END_REF], [START_REF] Benzoni-Gavage | On the well-posedness for the Euler-Korteweg model in several space dimensions[END_REF] and quantum mechanics [START_REF] Aftalion | Équation de Schrödinger non linéaire: des condensats de Bose-Einstein aux supersolides[END_REF], [START_REF] Hoefer | Dispersive and classical shock waves in Bose-Einstein condensates and gas dynamics[END_REF], [START_REF] Carles | Gross-Pitaevskii and Korteweg[END_REF], [START_REF] Jenkins | Regularization of a sharp shock by the defocusing nonlinear Schrödinger equation[END_REF], [START_REF] El | Dispersive shock waves and modulation theory[END_REF]. Here we will study the NLS equation (more exactly, its hydrodynamic form through the Madelung transform [START_REF] Madelung | Quantentheorie in Hydrodynamischer Form[END_REF]) which is mathematically equivalent to capillary fluids. An additional advantage compared to the general dispersive systems whose Lagrangian depends on the gradients of unknowns is that it is integrable [START_REF] Zakharov | On the complete integrability of a nonlinear Schroödinger equation[END_REF]. It allows us, in particular, to have a large family of exact solutions in order to quantify the accuracy of the approximation related to the introduction of the extended Lagrangian.

In section 2, a short reminder of the defocusing NLS equation and the Madelung transform is presented. In section 3, the extended Lagrangian formulation of this equation is given. Section 4 covers the main numerical resolution steps of the obtained equations. In section 5, we present two test cases. First, we compare the accuracy of our approach in the case of a 1-D 'grey' soliton which is an exact solution to NLS equation. Then, we compare our results for discontinuous initial data with solutions of modulation equations to the NLS equation (Whitham's equations) [START_REF] Gurevich | Nonstationary structure of a collisionless shock wave[END_REF], [START_REF] Gurevich | Expanding self-similar discontinuities and shock waves in dispersive hydrodynamics[END_REF], [START_REF] Gurevich | Dissipationless shock waves in media with positive dispersion[END_REF]. In both cases, a very good agreement of numerical and analytical results is obtained. Technical details are given in the Appendices.

In this paper, we shall use bold lowercase letters for vectors and bold uppercase letters for any second order tensor. For any vectors u and v, we shall denote u • v their scalar product and u ⊗ v their tensor product. For any scalar λ and vector u we shall denote ∇λ the gradient vector (column) of λ, ∂u ∂x the jacobian matrix of u. The divergence of a second order tensor A is the vector div(A) such that for any constant vector a: div(A) • a = div(Aa).

In particular, for any vector fields u, v:

div(Au) = div(A) • u + tr A ∂u ∂x , div (u ⊗ v) = vdiv (u) + ∂v ∂x u.
Also, for any function f , the notation ḟ means the material derivative of f along the velocity field u:

ḟ = ∂f ∂t + u • ∇f.

Defocusing NLS equation

Non-linear Schrödinger's equation (NLS) appears in a wide range of contexts, for example in non-linear optics [START_REF] Chiao | Self-trapping of optical beams[END_REF], surface gravity waves [START_REF] Grimshaw | Envelope solitary waves[END_REF], quantum fluids [START_REF] Roberts | The nonlinear Schrödinger equation as a model of superfluidity[END_REF], [START_REF] Hoefer | Dispersive and classical shock waves in Bose-Einstein condensates and gas dynamics[END_REF]. In general it serves as a fundamental model for description of non-linear and dispersive phenomena. Integrability of the NLS equation in both focusing (f (ρ) < 0) and defocusing (f (ρ) > 0) cases was shown in [START_REF] Zakharov | On the complete integrability of a nonlinear Schroödinger equation[END_REF]. A general form of the NLS equation can be written as:

i ψ t + 2 2m ∆ψ -mf |ψ| 2 ψ = 0. (1)
Here, f is a function of |ψ| 2 that characterizes the non-linearity of the equation. If we consider the parameter = m , the NLS equation can be written as :

i ψ t + 2 2 ∆ψ -f (|ψ| 2 )ψ = 0. (2)
In this equation, acts only as a scaling factor. In fact, if we consider the transformed variables x = x and t = t , we can rewrite (2) in the new variables as :

iψ t + 1 2 ∆ ψ -f (|ψ| 2 )ψ = 0, ( 3
)
where ∆ is the Laplacian with respect to the space coordinates x . In the following we take the scaling factor = 1.

Madelung's transform

Equation ( 2) can be rewritten in equivalent 'hydrodynamic' form. For this, let us consider Madelung's transform [START_REF] Madelung | Quantentheorie in Hydrodynamischer Form[END_REF], [START_REF] Kamchatnov | Nonlinear periodic waves and their modulations[END_REF] defined by : ψ(x, t) = ρ(x, t)e iθ(x,t) . (

Here, ρ ≥ 0 and θ are real valued functions of time and space. Plugging this expression into equation (1) permits us to separate imaginary and real parts of the equation. This yields the following system:

         ρ t + div(ρu) = 0 u t + (u • ∇) u + ∇ f (ρ) - ∆( √ ρ) 2 √ ρ = 0 (5)
Here, the considered 'velocity' field u is defined as :

u = ∇θ (6)
It is important to note that under such a definition curl(u) = 0. This form is said to be 'hydrodynamic' due to its similarity with the Euler equations. We will only consider the defocusing NLS equation for which the function f is given in a special form :

f (ρ) = ρ. ( 7
)
In this case, system (5) can be rewritten as :

       ρ t + div(ρu) = 0, (ρu) t + div (ρu ⊗ u + Π) = 0, (8)
with

Π = ρ 2 2 - 1 4 ∆ρ Id + 1 4ρ ∇ρ ⊗ ∇ρ. (9)
For such a set of equations, one can construct a Lagrangian, for which [START_REF] Besse | Discrete transparent boundary conditions for the mixed KdV-BBM equation[END_REF] are the Euler-Lagrange equations [START_REF]La capillarité interne[END_REF], [START_REF]Capillarité interne en mécanique des millieux continus[END_REF], [START_REF] Eglit | A generalization of the model of an ideal compressible fluid[END_REF], [START_REF] Casal | A representation of liquid-vapor interfaces by using fluids of second grade[END_REF], [START_REF] Gavrilyuk | Media with equations of state that depend on derivatives[END_REF], [START_REF] Benzoni-Gavage | On the well-posedness for the Euler-Korteweg model in several space dimensions[END_REF] :

L = Ωt ρ |u| 2 2 - ρ 2 2 - 1 4ρ |∇ρ| 2 2 dΩ. ( 10
)
Equations ( 8) admit the energy conservation law of the form :

∂E ∂t + div(Eu + Πu - 1 4 ρ∇ρ) = 0 (11)
where

E = ρ |u| 2 2 + ρ 2 2 + 1 4ρ |∇ρ| 2 2 (12)
Compared to classical models, for which the potential energy does not depend explicitly on ∇ρ, we have an additional term in the energy flux (-1 4 ρ∇ρ ) known as interstitial working [START_REF] Dunn | On the thermomechanics of interstitial working[END_REF].

Exact solutions

We are looking for traveling wave solutions :

ρ(x, t) = ρ(ξ) ; u(x, t) = u(ξ) ; ξ = x -U t. ( 13
)
Under this assumption, successive integrations of (8) yield the following equations for the hydrodynamic quantities ρ and u [START_REF] Gurevich | Dissipationless shock waves in media with positive dispersion[END_REF], [START_REF] El | Decay of an initial discontinuity in the defocusing NLS hydrodynamics[END_REF]:

         dρ dξ 2 = 4(ρ -b 1 )(ρ -b 2 )(ρ -b 3 ) ; u = U + q ρ U 2 = b 1 + b 2 + b 3 ; q 2 = b 1 b 2 b 3 (14)
In these expressions, b i are arbitrary constants. Without loss of generality, we can assume that b 1 > b 2 > b 3 > 0. One can obtain a solution of (8) in terms of the Jacobi elliptic function dn :

ρ = b 1 -(b 1 -b 3 )dn 2 b 1 -b 3 (x -U t) , s , u = U + q ρ . ( 15
)
Here s is the parameter of the elliptic function satisfying the relation :

s 2 = b 2 -b 3 b 1 -b 3 , 0 < s 2 < 1, (16)
and dn(v, s) = 1 -s 2 sn 2 (v, s), sn(v, s) = sin(ϕ(v, s)),

where ϕ(v, s) is obtained implicitly from the equation :

v = ϕ 0 dθ 1 -s 2 sin 2 (θ)
.

For each fixed value of 0 < s < 1, solution ( 15) is a periodic wave of amplitude a and wavenumber k given by :

a = b 2 -b 3 2 , k = π K(s) 2a s 2 . ( 17
)
Here K(s) is the complete elliptic integral of the first kind defined as :

K(s) = π 2 0 dθ 1 -s 2 sin 2 (θ) . ( 18
)
For the limiting values s → 0 and s → 1, the oscillations vanish and two different behaviors are observed for each limit. In the case s → 0, the amplitude of the oscillations vanishes. The wave number remains at a finite value given by k → 2 √ b 1 -b 3 . In the case s → 1, we get K(s) → ∞ and therefore the oscillations vanish as their wavenumber k → 0. Since dn(v, s)| s→1 = 1/cosh(v) the solution behaves like a soliton of amplitude (b 1 -b 3 ).

Below, we will present briefly some exact and asymptotic solutions to the NLS equation based on the representation [START_REF] Benzoni-Gavage | On the well-posedness for the Euler-Korteweg model in several space dimensions[END_REF].

One-dimensional 'grey' soliton

The first case corresponds to the family of solutions obtained from (15) in the limit s → 1 :

ρ = b 1 - b 1 -b 3 cosh 2 √ b 1 -b 3 (x -U t) u = U - b 1 √ b 3 ρ (19)
This family of solutions is called 'grey' solitons [START_REF] Aftalion | Équation de Schrödinger non linéaire: des condensats de Bose-Einstein aux supersolides[END_REF]. It consists of a localized wave pulse of amplitude (b 1 -b 3 ) that propagates with a constant velocity U while maintaining its shape. The parameters b 1 and b 3 define the limit values of the soliton such that :

lim |x|→∞ ρ(x, t) = b 1 ; ρ min = b 3 ; lim |x|→∞ u(x, t) = U -b 3 ; u min = U - b 1 √ b 3 . ( 20
)
The shape of the soliton is shown on Figure 1.

Dispersive shock waves

In classical hydrodynamics described by dispersionless hyperbolic Euler equations, shock waves (strong discontinuities) can appear. In the case of dispersive hydrodynamics, these singularities are resolved by the apparition of an oscillatory wave train in a region of space which expands over time. This is referred to as a dispersive shock wave (DSW) [START_REF] Gurevich | Nonstationary structure of a collisionless shock wave[END_REF], [START_REF] Gurevich | Expanding self-similar discontinuities and shock waves in dispersive hydrodynamics[END_REF], [START_REF] El | Decay of an initial discontinuity in the defocusing NLS hydrodynamics[END_REF], [START_REF] El | Unsteady undular bores in fully nonlinear shallowwater theory[END_REF], [START_REF] El | Dispersive shock waves and modulation theory[END_REF]. The amplitude of the DSW is described by Whitham's modulation equations [START_REF] Whitham | Linear and nonlinear waves[END_REF].

Let us consider the Riemann problem for NLS equation :

   ρ(x) = ρ L x < 0 ρ(x) = ρ R x > 0    u(x) = u L x < 0 u(x) = u R x > 0 (21)
Such an initial discontinuity gives rise to either a DSW or a rarefaction wave on each side. The evolution of amplitude of the DSW is described in terms of the Riemann invariants r i (r 1 > r 2 > r 3 > r 4 ) of the corresponding Whitham's averaged equations [START_REF] Pavlov | Nonlinear Schroödinger equation and the Bogolyubov-Whitham method of averaging[END_REF], [START_REF] Gurevich | Dissipationless shock waves in media with positive dispersion[END_REF] :

∂r i ∂t + V i ∂r i ∂x = 0, i = 1, 2, 3, 4. (22)
The characteristic velocities are given by:

V 1 = U (r) + 1 2 (r 1 -r 2 ) 1 - r 2 -r 4 r 1 -r 4 E(s) K(s) -1 , V 2 = U (r) - 1 2 (r 1 -r 2 ) 1 - r 1 -r 3 r 2 -r 3 E(s) K(s) -1 , V 3 = U (r) + 1 2 (r 3 -r 4 ) 1 - r 2 -r 4 r 2 -r 3 E(s) K(s) -1 , V 4 = U (r) - 1 2 (r 3 -r 4 ) 1 - r 1 -r 3 r 1 -r 4 E(s) K(s) -1
, where the complete elliptic integral of the second kind E(s) is defined as

E(s) = π 2 0 1 -s 2 sin 2 (θ)dθ.
The variables b i and U are linked to the r i via the relations [START_REF] Pavlov | Nonlinear Schroödinger equation and the Bogolyubov-Whitham method of averaging[END_REF], [START_REF] Gurevich | Dissipationless shock waves in media with positive dispersion[END_REF] :

b 1 = 1 16 (r 1 + r 2 -r 3 -r 4 ) 2 , ( 23
) b 2 = 1 16 (r 1 + r 3 -r 2 -r 4 ) 2 , ( 24
) b 3 = 1 16 (r 1 + r 4 -r 2 -r 3 ) 2 , ( 25
)
U = 1 4 (r 1 + r 2 + r 3 + r 4 ) , ( 26
)
s 2 = (r 1 -r 2 )(r 3 -r 4 ) (r 1 -r 3 )(r 2 -r 4 ) . ( 27
)
If we consider a self-similar evolution of the Riemann invariants, r i depend only on τ = x/t and equations [START_REF] Gurevich | Nonstationary structure of a collisionless shock wave[END_REF] reduce to :

r i (V i -τ ) = 0 ( 28
)
It means that one of the Riemann invariants, r j changes in space and time and its characteristic speed is V j = τ while the three other invariants are constants determined by the the initial conditions on both sides of the initial discontinuity. Outside of the DSW region, it was shown in [START_REF] Gurevich | Dissipationless shock waves in media with positive dispersion[END_REF] that the Whitham equations for the defocusing NLS degenerate into Euler equations for shallow water flows. This transition occurs when two of the Riemann invariants r i merge together, leading to either s 2 = 0 or s 2 = 1. The two remaining invariants behave like Riemann invariants for the shallow water equations :

r ± = u ± 2 √ ρ ( 29
)
The solution is then determined via matching of the Riemann invariants at the DSW fronts [START_REF] Gurevich | Nonstationary structure of a collisionless shock wave[END_REF], [START_REF] Gurevich | Dissipationless shock waves in media with positive dispersion[END_REF].

Since this procedure depends strongly on the structure of the flow, we will consider the case in which a rarefaction wave on the left and a DSW on the right are created (see Figure 2). At the oscillatory front (corresponding to the leading edge of the DSW, s 2 = 0) we get [START_REF] El | Decay of an initial discontinuity in the defocusing NLS hydrodynamics[END_REF], [START_REF] El | Dispersive shock waves and modulation theory[END_REF]:

r 1 = r 2 , V 1 = V 2 , r 3 = r + (R), r 4 = r -(R). ( 30
)
At the soliton front (corresponding to the trailing edge of the DSW, s 2 = 1) we get :

r 2 = r 3 , V 2 = V 3 , r 1 = r + (0), r 4 = r -(0). ( 31
)
Here, r ± (0) and r ± (R) are the values of the invariants r ± at the states '0' (the constant state found solving the Riemann problem for the non-dispersive shallow water equations) and 'R' (the state on the right of the initial discontinuity). The profile of the solution is shown in Figure 2. The oscillatory part of the solution is plotted according to the following algorithm :

1. Set the values of ρ L , ρ R , u L , u R . 2. Calculate the values of ρ 0 u 0 , r 1 , r 3 , r 4 . 3. Calculate r 2 (s) = (r 1 (r 3 -r 4 ) + s 2 (r 1 -r 3 )r 4 )/(r 3 -r 4 + (r 1 -r 3 )s 2 ). 4. Calculate the functions b 1 (s), b 2 (s), b 3 (s), U (s). 5. Calculate τ (s) = V 2 (s).
6. Choose a time instant t. The DSW is shown as a parametric plot of (15

) : ρ(s, t) = ρ(τ (s), t), 0 ≤ s ≤ 1.
7. The low and upper boundaries of the oscillatory profile are described by :

ρ inf (s) = b 3 (s), ρ sup (s) = s 2 b 1 (s) + (1 -s 2 )b 3 (s). τ 2 τ 1 τ 3 τ 4 τ=x/t ρ 0 ρ R ρ L ρ
τ 1 = u R + 8ρ 0 -8 √ ρ 0 ρ R + ρ R 2 √ ρ 0 - √ ρ R , τ 2 = u R + √ ρ 0 , τ 3 = u 0 - √ ρ 0 , τ 4 = u L - √ ρ L . τ 2 τ 1 τ = x/t r + (0) r + (R) r 2
Figure 3: Asymptotic profile of invariant r 2 for Whitham's system. In the DSW region τ 2 < τ < τ 1 , r 2 (dashed line) varies while the other r i are constants.

Extended Lagrangian formulation

In this part, the extended Lagrangian formulation is given for the NLS equation. The Euler-Lagrange equations and hyperbolicity study are given in Appendices A, B in a more general framwork, from which the application to NLS equation can be easily retrieved. This general approach will allow us to apply our method to other physical problems with the density gradient dependent energies.

Let us consider a new dual to ρ variable η The idea is to substitute ∇ρ by ∇η and to guarantee the convergence of η to ρ in a certain limit. In order to do that, let us consider the extended Lagrangian :

L e = Ωt L e (u, ρ, η, η, ∇η) dΩ, ( 32
)
where a two-parameter family of local Lagrangians L e is given by :

L e (u, ρ, η, η, ∇η) = ρ |u| 2 2 + β 2 ρ η2 - ρ 2 2 - 1 4ρ |∇η| 2 2 - λ 2 ρ η ρ -1 2 . ( 33
)
This extended Lagrangian is constructed from the original Lagrangian (10) by replacing ∇ρ by ∇η and adding two terms :

• The first term λ 2 ρ η ρ -1 2
is a penalization term. When λ → ∞, the difference (η/ρ -1) vanishes.

• The second term β 2 ρ η2 is necessary in order to regularize the time evolution of η and to ensure the hyperbolicity of the new governing equations. It can be perceived as a small micro-inertia term.

We can write now the extended Hamilton's action:

a = t1 t0 L e dt = t1 t0 Ωt L e (u, ρ, η, w, p) dΩ dt, ( 34
)
where w = η, p = ∇η.

Applying a technique for variation of the functionals with differential constraints (in our case, this is the mass conservation law) [START_REF] Dell'isola | Variational Models and Methods in Solid and Fluid Mechanics[END_REF], one can obtain the following system of equations, written in conservative form (for details, see Appendix A):

∂ρ ∂t + div(ρu) = 0, (35)
∂ρu ∂t + div ρu ⊗ u + ρ 2 2 - 1 4ρ |p| 2 + ηλ(1 - η ρ ) Id + 1 4ρ p ⊗ p = 0, (36)
∂ρη ∂t + div(ρηu) = ρw, ( 37
)
∂ρw ∂t + div ρwu - 1 4ρβ p = λ β 1 - η ρ , ( 38
)
∂p ∂t + div ((p • u -w) Id) = 0; curl(p) = 0. ( 39
)
This system admits the energy conservation law :

∂E e ∂t + div E e u + Pu - 1 4ρ wp = 0,
where

E e = ρ |u| 2 2 + β 2 ρ η2 + ρ 2 2 + 1 4ρ |p| 2 2 + λ 2 ρ η ρ -1 2 , P = ρ 2 2 - 1 4ρ |p| 2 + ηλ(1 - η ρ ) Id + 1 4ρ p ⊗ p.
It is interesting to note that, even for the extended system, the energy conservation equation always contains, as in the case of original system (see [START_REF]Capillarité interne en mécanique des millieux continus[END_REF]) the interstitial working term written here in the form -1 4ρ wp.

The main point in using the approach of an extended Lagrangian is a possibility to replace the original dispersive system by an approximate hyperbolic system of equations [START_REF] Favrie | A rapid numerical method for solving Serre-Green-Naghdi equations describing long free surface gravity waves[END_REF]. This will allow us to use standard and rapid numerical methods for hyperbolic systems of equations. The system (35)-(39) is hyperbolic (see proof in Appendix B). The characteristic speeds c in x-direction are given by:

c = u, (c -u) 2 ± = 1 4βρ 2 + ρ + λη 2 ρ 2 ± | -1 4βρ 2 + ρ + λη 2 ρ 2 | 2 . ( 40
)
Linearizing the governing equations on the constant solution ρ = ρ 0 , u = 0, w = 0, p = 0, η = ρ 0 and looking for the solutions which are proportional to e i(kx-ωt) , where k is the wave number and ω is the frequency, one can obtain the dispersion relation expressed here in the form c p = c p (k), where c p = ω k is the phase velocity :

(c ± p ) 2 = 1 4βρ 2 0 + ρ 0 + λ + λ βρ 2 0 k 2 ± 1 4βρ 2 0 + ρ 0 + λ + λ βρ 2 0 k 2 2 -4 λ βρ 0 k 2 + ρ 0 + λ 4βρ 2 0 2 . ( 41
)
One can easily see that in the linear case the phase velocity approaches to the characteristic velocity in the limit k → ∞.

Estimation of β and λ

The extended Lagrangian approach (32) requires a choice of β and λ. This choice can be based, for example, on the fact that the dispersion relations for both extended and original systems must remain close in a specific range of wave numbers. For the equilibrium state defined by ρ = ρ 0 , u = u 0 = 0, the dispersion relation for the original Schrödinger equation is :

c 2 p = ρ 0 + k 2 /4. ( 42
)
When (β, λ) → (0, +∞), the convergence of these dispersion relations is not uniform, that is, the curves almost coincide for low wave numbers but start to stray away from each other starting from a certain threshold wave number k max (β, λ). This threshold wave number must be chosen such that the wave numbers that may be present in the solution are contained in [0, k max ].

Numerical resolution

We consider the numerical resolution of equations (39) in the 1-D case. The system can be written as :

∂U ∂t + ∂F ∂x = S, (43)
where U, F and S are respectively the vector of conservative variables, flux vector and source term given by :

U =       ρ ρu ρη ρw p       , F =       ρu ρu 2 + ρ 2 2 + ηλ(1 -η ρ ) ρηu ρwu -1 4ρβ p pu -w       , S =        0 0 ρw λ β 1 -η ρ 0        .
Due to a non-zero term source, a splitting strategy is applied [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF], [START_REF] Favrie | A rapid numerical method for solving Serre-Green-Naghdi equations describing long free surface gravity waves[END_REF]. Hence, at each time step, the numerical resolution is split into a hyperbolic part (no source terms) :

∂U ∂t + ∂F ∂x = 0, (44)
and an ordinary differential equation (ODE) part :

dU dt = S. (45)

Hyperbolic step

Equation ( 44) is solved by using a classical Godunov type scheme for hyperbolic equations :

U n+1 i = U n i + ∆t ∆x F * i+ 1 2 -F * i-1 2 , 1 ≤ i ≤ N. ( 46
)
Here N is the number of cells of uniform size ∆x, and ∆t is the time step satisfying the Courant -Friedrichs-Levy condition. The intercell fluxes F * i+ 1 2 are computed using the HLL-Rusanov approximate Riemann solver [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF] :

F * i+ 1 2 U n i , U n i+1 = 1 2 F(U n i+1 ) -F(U n i ) - 1 2 κ n i+ 1 2 U n i+1 -U n i , ( 47
)
where

κ n i+ 1 2
is by using the Davis approximation [START_REF] Davis | Simplified second-order Godunov-type methods[END_REF] :

κ n i+1/2 = max j (|c j (U n i )|, |c j (U n i+1 )|), (48)
where c j are the eigenvalues of the extended system.

ODE step

The source terms treatment is reduced to a second order ordinary differential equation with constant coefficients which can be solved exactly in our case. The relaxation equations are given by :

                               dρ dt = 0 dρu dt = 0 dρη dt = ρw dρw dt = λ β 1 - η ρ dp dt = 0
Therefore, the exact solution is given by :

           ρ n+1 = ρn u n+1 = ūn p n+1 = pn η n+1 = ρn + (η n -ρn ) cos(Ω∆t) + wn Ω sin(Ω∆t) w n+1 = Ω(ρ n -ηn ) sin(Ω∆t) + wn cos (Ω∆t)
where Ω = λ βρ 2 .

Results

All numerical results of this section are obtained by the MUSCL-Hancock extension to the Godunov scheme with MIN-MOD limiter [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF].

One-dimensional grey soliton

We consider the initial conditions corresponding to the solitary wave solutions of NLS equation :

ρ(x, 0) = b 1 - b 1 -b 3 cosh 2 √ b 1 -b 3 x u(x, 0) = U - b 1 √ b 3 ρ(x, 0) , ( 49
) η(x, 0) = ρ(x, 0), w(x, 0) = -ρ(x, 0)u x (x, 0), p(x, 0) = ρ x (x, 0).
We use periodic boundary conditions in the computational domain. Thus we impose the following relation between numerical fluxes in the first and last numerical cells:

F * 1 2 = F * N + 1 2 = F * (U n N , U n 1 ) . ( 50
)
It means that when the soliton passes through one of the boundaries, it continuously reappears on the other side. This permits us to run simulations for a longer time without having to use large domains.

The period of such a configuration is the required time for the soliton to reach back its initial position.

The results are shown on Figure 5: One can see that the shape and the position of the simulated soliton are in perfect agreement with the exact solution. Furthermore, Figure 6 shows the solution for different mesh sizes ∆x = 0.004, ∆x = 0.002, ∆x = 0.0008 and ∆x = 0.0002, with a focus on the soliton peak. We can see that for more refined meshes, we get more accurate values for the soliton position and amplitude, which highlights the convergence of the scheme.

Dispersive Riemann problem

We present here the results for dispersive Riemann problem. In order to get the case where a DSW propagates to the right and a rarefaction wave to the left, the initial values must satisfy the condition [START_REF] El | Decay of an initial discontinuity in the defocusing NLS hydrodynamics[END_REF] :

r + (L) > r + (R) > r -(R) > r -(L) (51)
From the numerical point of view, the initial step function induces an infinite gradient which causes spurious effects to appear near the initial discontinuity. Hence, we will use a rather regular step-like function defined by :

ρ 0 (x) = ρ L + ρ R 2 + ρ L -ρ R 2 tanh x δ , ( 52
)
u 0 (x) = u L + u R 2 + u L -u R 2 tanh x δ . ( 53
)
This spreads the abrupt discontinuity over the region [-2δ, 2δ] as shown in Figure ( 7) : We plot the quantities ρ and u, obtained numerically, as functions of the auto-similar variable x/t. The results are shown in figure [START_REF] Besse | Discrete transparent boundary conditions for the mixed KdV-BBM equation[END_REF] :

x ρ -2δ 2δ 0 ρ L ρ R x smoothed step theoretical step
x/t ρ τ 4 τ 3 τ 2 τ 1 ρ L ρ R ρ 0 x/t numerical simulation Whitham envelope x/t u τ 4 τ 3 τ 2 τ 1 u L u R u 0 numerical simulation Whitham envelope
ρ 0 = 1 4 ( √ ρ R + √ ρ L + 1 2 (u L -u R )) 2 and u 0 = 1 2 (u L +u R )+ √ ρ L - √ ρ R . The initial values used for this simulation are ρ L = 2, ρ R = 1, u L = u R = 0.
The parameters are: β = 2.10 -5 , λ = 300, ∆x = 0.000667.

Clearly, the overall structure of the solution complies with the asymptotic one shown on Figure 2. The amplitude of the oscillations shows a very good agreement with the asymptotic DSW profile from Whitham's theory of modulations. The oscillations develop in both regions: τ > τ 1 and τ < τ 1 . These oscillations are part of the solution but they do not appear in the asymptotic (in time) limit because they vanish when t → ∞ as their amplitude decreases with time as a ∝ t -1/2 [START_REF] Gurevich | Expanding self-similar discontinuities and shock waves in dispersive hydrodynamics[END_REF]. This behavior is also displayed with the same power law in our results. By measuring the amplitude a of the first oscillation at the vicinity of τ 4 , we could plot the function f (t) = at 3/2 (see Figure [START_REF] Favrie | A rapid numerical method for solving Serre-Green-Naghdi equations describing long free surface gravity waves[END_REF] . The plot shows that f (t) is a linear function of time which implies that a √ t is constant. Another important detail is to check the position of the soliton which arises at the vicinity of τ = τ 2 . Let τ s (t) be the soliton position at time t. We will compare τ s (t) with the theoretical position τ 2 . The relative error is given by :

err pos (t) = τ s (t) -τ 2 τ 2 . ( 54
)
The time evolution of this error is shown on Figure

Conclusions and perspectives

In order to solve the NLS equation, a new numerical method based on the formulation of the 'extended Lagrangian' approach was developed. The Euler-Lagrange equations corresponding to this extended Lagrangian are hyperbolic and approximate the NLS equation with a good accuracy, when it was compared with exact and approximate solutions to the NLS equations. The equations being hyperbolic grants the possibility of developing classical Godunov type schemes which makes computational time reasonable and opens a possibility to extend this work to the multidimensional case. Indeed, the extended model is hyperbolic, the equations are local and the scheme could be parallelized very efficiently using domain decomposition. This property may also simplify the development of 'transparent' boundary conditions which prevent the reflection of oscillations and small perturbations that reach the boundaries.

A Extended Lagrangian

Let us consider Hamilton's action :

a = t1 t0 L e dt (55)
where L e is a generic extended Lagrangian given by :

L e = Ωt ρ |u| 2 2 -W (ρ, η, η, ∇η) dΩdt. ( 56
)
Here [t 0 , t 1 ] is the time interval, Ω t is the material volume, and W is a generic potential. In particular, for the approximation of the Schrödinger equation it will be chosen as :

W (ρ, η, η, ∇η) = ρ 2 2 + 1 4ρ |∇η| 2 2 + λρ 2 η ρ -1 2 - βρ 2 η2 . ( 57
)
However, in the following we will consider a generic function W because the calculations are less cumbersome in a general case. The variation of Hamilton's action is submitted to the constraint which is the mass conservation law :

∂ρ ∂t + div (ρu) = 0. ( 58
)
Two types of variations will be considered : one corresponds to the variation of η as independent variable (and will be denoted by δη), and the second one corresponds to the virtual displacement of the continuum, respecting the mass conservation law (this variation will be denoted by δx). In particular, in the last case the Eulerian variations of the density, velocity and η are related to δx as (see for details [START_REF] Dell'isola | Variational Models and Methods in Solid and Fluid Mechanics[END_REF]):

δρ = -div (ρδx) , δu = (δx) - ∂u ∂x δx, δ η = δu • ∇η. (59)
Here 'hat' means to the variation at fixed Eulerian variables, and 'dot' denotes the material derivative. Both variations are vanishing at the boundary of [t 0 , t 1 ]×Ω t . We consider first the variation of Hamilton's action with respect to η. We need to express the variation of η :

δ η = δ ∂η ∂t + u • δ∇η. ( 60
)
Schwarz's theorem allows us to interchange the order of derivatives :

δ ∂η ∂t = ∂δη ∂t , δ∇η = ∇(δη). (61)
Hence :

δ η a = t1 t0 Ωt δ η ρ |u| 2 2 -W (ρ, η, η, ∇η) dΩdt = t1 t0 Ωt - ∂W ∂η δη + ∂W ∂ η δ η + ∂W ∂∇η δ∇η dΩdt = t1 t0 Ωt - ∂W ∂η δη + ∂W ∂ η δ ∂η ∂t + u • δ∇η + ∂W ∂∇η δ∇η dΩdt.
Using the Gauss-Ostrogradski's theorem and taking into account the fact that the variation δη is vanishing on the boundary of [t 0 , t 1 ] × Ω t , we finally obtain : This yields classical Euler-Lagrange equation for η:

δ η a = t1 t0 Ωt - ∂W ∂η + ∂ ∂t ∂W ∂ η + div ∂W ∂ η u + div
- ∂W ∂η + ∂ ∂t ∂W ∂ η + div ∂W ∂ η u + div ∂W ∂∇η = 0. ( 62

B Hyperbolicity of the extended Lagrangian model

We consider the 1D case where all the variables are functions of only (x, t) instead of (x = (x, y, z) T , t).

The velocity field u and the gradient p of w will be denoted by u and p, respectively. The governing equations for the extended Lagrangian model can be rewritten as : In order to study the hyperbolicity of this system, we write it in quasi-linear form :

∂ρ
∂U ∂t + A(U) ∂U ∂x = q (75)
where:

U =       ρ u w p η       A(U) =        u ρ 0 0 0 1 + λη 2 ρ 3 u 0 0 λ ρ 1 -2η ρ p 4βρ 3 0 u -1 4βρ 2 0 0 p -1 u 0 0 0 0 0 u        q =        0 0 1λ βρ 1 -η ρ 0 w       
The eigenvalues c of the matrix A are :

c = u, (c -u) 2 ± = 1 4βρ 2 + ρ + λη 2 ρ 2 ± -1 4βρ 2 + ρ + λη 2 ρ 2 2 2 . ( 76
)
The right-hand side of (76) is always positive. However, the roots can be multiple if

1 4βρ 2 = ρ + λη 2 ρ 2 .
One can easily prove that even in the case of multiple roots one always has five linear independent eigenvectors. So, the system is always hyperbolic.

Linearizing the governing equations on the constant solution ρ = ρ 0 , u = 0, w = 0, p = 0, η = ρ 0 and looking for the solutions which are proportional to e i(kx-ωt) , where k is the wave number and ω is the frequency, one can obtain the dispersion relation expressed here in the form c p = c p (k), where c p = ω k is the phase velocity :

(c ± p ) 2 = 1 4βρ 2 0 + ρ 0 + λ + λ βρ 2 0 k 2 ± 1 4βρ 2 0 + ρ 0 + λ + λ βρ 2 0 k 2 2 -4 λ βρ 0 k 2 + ρ 0 + λ 4βρ 2 0 2 . ( 77
)

Figure 1 :

 1 Figure 1: Different shapes of the variables ρ and u for the 'grey' soliton solution, for different values of the parameters b 1 and b 3 at t = 0

Figure 2 :

 2 Figure 2: Asymptotic profile of the solution to NLS equation (continuous line) for the Riemann problem ρ L = 2, ρ R = 1 , u L = u R = 0. The boundaries τ i , i = 1, 2, 3, 4 delimit the DSW and the rarefaction wave regions. The modulation of the DSW profile between τ 2 and τ 1 is described by the rarefaction wave solution to the Whitham system (bold dashed line). The oscillatory profile is shown at t = 70. The values of τ i , i = 1, 2, 3, 4, are given by τ 1= u R + 8ρ 0 -8 √ ρ 0 ρ R + ρ R

Figure 4 :

 4 Figure 4: The dispersion relation (42) (continous line) and (41) for the extended Lagrangian (dashed lines) for for β = 10 -4 and different values of λ.

Figure 5 :

 5 Figure 5: Numerical profiles of ρ (left) and u (right) for the grey soliton at t = 0 (dot-dashed line) and at t = 2T (continuous line). The used domain is L = [-20, 20] with ∆x = 0.0002, the period is T = D/U = 20. Parameters used for the simulation are b 1 = 1.5, b 3 = 1, U = 2, = 1, β = 10 -4 , λ = 500.

Figure 6 :

 6 Figure6: Magnified view on the numerical 'grey' soliton peak, for different mesh sizes, at t = 2T . The theoretical amplitude and position are ρ = 1 and x = 0. Parameters used for the simulation are the same as in figure[START_REF] Gurevich | Dissipationless shock waves in media with positive dispersion[END_REF].

Figure 7 :

 7 Figure 7: Magnified view over the smoothed step (continuous line) for δ = 0.1, ∆x = 0.000667.

Figure 8 :

 8 Figure 8: Comparison of the numerical results (thin line) with the Whitham modulational profile of the DSW (thick line) at t = 70. The upper figure shows ρ = f (x/t) and the lower figure shows u = f (x/t). The displayed τ i are the theoretical boundaries of the DSW and the rarefaction wave. The values of ρ 0 and u 0 are the theoretical values in the central plateau given by ρ0 = 1 4 ( √ ρ R + √ ρ L + 1 2 (u L -u R )) 2 and u 0 = 1 2 (u L +u R )+ √ ρ L -√ ρ R .The initial values used for this simulation are ρ L = 2, ρ R = 1, u L = u R = 0.

Figure 9 :

 9 Figure 9: Vanishing oscillations at the vicinity of the singular point τ = τ 4 . The figure on the left shows that these oscillations decrease in time. The figure on the right shows that the amplitude of the first oscillation a is such that at 3/2 is linear. This implies the power law a ∝ t -1/2 .

Figure 10 :

 10 Figure 10: Relative error on the position of the first soliton plotted as function of time.

  ∂W ∂∇ηδη dΩdt = 0, for any δη.

)

  Now, we consider the variation of Hamilton's action with respect to δx. Using (59) and Gauss-Ostrogradski's formula, we will arrive to :

	δa = -	t1 t0	Ωt	∂ρu ∂t	+ div ρu ⊗ u + ρ	∂W ∂ρ	-W Id +	∂W ∂∇η	⊗ ∇η	• δxdΩdt = 0, for any δx.
	This implies :
					∂ρu ∂t	+ div ρu ⊗ u + ρ	∂W ∂ρ	-W Id +	∂W ∂∇η	⊗ ∇η = 0	(63)
	Thus we have the following system of equations :
					∂ρ ∂t	+ div(ρu) = 0
					∂ ∂t		∂W ∂ η	+ div	∂W ∂ η u + div	∂W ∂∇η	-	∂W ∂η	= 0
					∂ρu ∂t	+ div ρu ⊗ u + ρ	∂W ∂ρ	-W Id +	∂W ∂∇η	⊗ ∇η = 0
	Adding the new variables :
											w = η =	∂η ∂t	+ u • ∇η, p = ∇η,	(64)
	we can rewrite the governing equations in the form :
						∂ρ ∂t	+ div(ρu) = 0	(65)
						∂ρu ∂t	+ div ρu ⊗ u + ρ	∂W ∂ρ	-W Id +	∂W ∂p	⊗ p = 0	(66)
						∂η ∂t	+ u • ∇η = w	(67)
						∂ ∂t		∂W ∂w	+ div	∂W ∂w	u +	∂W ∂p	-	∂W ∂η	= 0	(68)
						∂p ∂t	+	∂p ∂x	u +	∂u ∂x	T	p -∇w = 0, curl(p) = 0.	(69)
	In the particular case (57), we obtain :
			∂ρ ∂t		+ div(ρu) = 0	(70)
			∂ρu ∂t	+ div ρu ⊗ u +	ρ 2 2	-	1 4ρ	|p| 2 + ηλ(1 -	η ρ	) Id +	1 4ρ	p ⊗ p = 0	(71)
			∂ρη ∂t	+ div(ρηu) = ρw	(72)
			∂ρw ∂t	+ div ρwu -	1 4ρβ	p =	λ β	1 -	η ρ	(73)
			∂p ∂t		+ div ((p • u -w) Id) = 0; curl(p) = 0.	(74)

Acknowledgments This work was partially supported by l'Agence Nationale de la Recherche, France (grant numbers ANR-11-LABEX-0092, and ANR-11-IDEX-0001-02).