Alexandre Poulain
email: alexandre.poulain@sorbonne-universite.fr

Katharina Schratz
email: katharina.schratz@sorbonne-universite.fr

Convergence, error analysis and longtime behavior of the Scalar Auxiliary Variable method for the nonlinear Schrödinger equation

Keywords: 2010 Mathematics Subject Classification. 35Q55, 65M70, 65M12, 65M15 Keywords and phrases. Nonlinear Schrödinger equation, Scalar Auxiliary Variable, Pseudospectral method, Error bounds

come

Convergence, error analysis and longtime behavior of the Scalar Auxiliary Variable method for the nonlinear

Schrödinger equation

Introduction

We consider the Gross-Pitaevskii [START_REF] Gross | Structure of a quantized vortex in boson systems[END_REF] equation (NLS) set on the d-dimensional torus Ω = T d = (R/2πZ) d (where d ≤ 3) i∂ t u(t, x) = -∆u(t, x) + V (x)u(t, x) + f |u(t, x)| 2 u(t, x), t ∈ (0, T ] (

with initial conditions u(0, x) = u 0 (x), a real-valued interaction potential V (x) and nonlinearity f (|u| 2 ).

The Hamiltonian energy associated to equation (1.1) takes the form

H(u, u) = 1 2 Ω |∇u| 2 + V (x) |u| 2 + F |u| 2 dx,
where F |u| 2 is defined by F |u| 2 = f |u| 2 . Note that the Hamiltonian H(u(t), u(t)) as well as the probability density u(t, •) 2 L 2 (Ω) is preserved by the system (1.1). In the following we will denote by E 1 the sum of the nonlinear and potential part of the Hamiltonian

E 1 = 1 2 Ω V (x) |u| 2 + F |u| 2 dx.
Using the decomposition u(t, x) = p(t, x) + iq(t, x), equation (1.1) can be furthermore rewritten as the Hamiltonian system ∂ t p = -∆q + δE 1 [t] δq , ∂ t q = ∆p -δE 1 [t] δp ,

with the associated Hamiltonian

H(p, q) = 1 2 Ω |∇p| 2 + |∇q| 2 + V (x) |p| 2 + |q| 2 + F |p| 2 , |q| 2 dx. (1.3)
In this notation, E 1 takes the form

E 1 = 1 2 Ω V (x) |p| 2 + |q| 2 + F |p| 2 , |q| 2 dx.
Due to their importance in numerous applications, reaching from Bose-Einstein condensation over nonlinear optics up to plasma physics, nonlinear Schrödinger equations are nowadays very well studied numerically. In the last decades a large variety of different numerical schemes has been proposed [START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF][START_REF] Antoine | Computational methods for the dynamics of the nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF][START_REF] Bao | Numerical solution of the Gross-Pitaevskii equation for Bose-Einstein condensation[END_REF][START_REF] Gonzales | Integral equation method for the continuous spectrum radial Schrödinger equation[END_REF][START_REF] Gonzales | Integral equation method for coupled Schrödinger equations[END_REF]. Thanks to their simplicity and accuracy, a popular choice thereby lies in so-called splitting methods, where the right hand side of (1.1) is split into the linear and nonlinear part, respectively, see, e.g., [START_REF] Bao | On time-splitting spectral approximations for the Schrödinger equation in the semiclassical regime[END_REF][START_REF] Besse | Order estimates in time of splitting methods for the nonlinear Schrödinger equation[END_REF][START_REF]The nonlinear Schrödinger equation and applications in Bose-Einstein condensation and plasma physics[END_REF] and the references therein. The popularity of splitting methods also stems from their structure preservation. They conserve exactly the L 2 norm of the solution and allow for near energy conservation over long times, see, e.g., [START_REF] Faou | Geometric numerical integration and Schrödinger equations[END_REF]. However, in [START_REF] Ostermann | Low regularity exponential-type integrators for semilinear Schrödinger equations[END_REF] the authors show that in certain applications splitting methods suffer from severe order reduction such as in case of non-linearities with non-integer exponents. The latter arises for instance in context of optical dark and power law solitons with surface plasmonic interactions [START_REF] Crutcher | Derivation of the Effective Nonlinear Schrödinger Equations for Dark and Power Law Spatial Plasmon-Polariton Solitons Using Nano Self-Focusing[END_REF]. As a solution to that issue, the authors proposed in [START_REF] Ostermann | Low regularity exponential-type integrators for semilinear Schrödinger equations[END_REF] a new class of low regularity exponential-type integrators for NLS. In this article we use a different approach based on the so-called Scalar Auxiliary Variable (SAV) method which was originally proposed to design structure-preserving numerical schemes for gradient flows [START_REF] Shen | The scalar auxiliary variable (SAV) approach for gradient flows[END_REF][START_REF]A new class of efficient and robust energy stable schemes for gradient flows[END_REF]. Very recently it also became popular in context of Hamiltonian systems [START_REF] Antoine | Scalar Auxiliary Variable/Lagrange multiplier based pseudospectral schemes for the dynamics of nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF][START_REF] Fu | A structure-preserving algorithm for the fractional nonlinear Schrödinger equation based on the sav approach[END_REF][START_REF] Cai | Structure-preserving algorithms for the twodimensional sine-Gordon equation with Neumann boundary conditions[END_REF][START_REF] Feng | High-order mass-and energy-conserving sav-gauss collocation finite element methods for the nonlinear schrödinger equation[END_REF]. The main advantage of the SAV method lies in the fact that it preserves a modified Hamiltonian on the discrete level. Due to its generality, it can be applied to a large class of equations involving any kind of nonlinearity. The resulting numerical schemes are linearly implicit and allow for efficient calculations.

The main idea behind the SAV method is to introduce a scalar variable r(t) = √ E 1 + E c that will become an unknown at the discrete level and where the arbitrary constant E c > 0 is used to obtain E 1 + E c > 0. We must stress that one as to be very careful with the choice of the constant E c . Indeed, it is well known that even for the cubic non-linearity, i.e. f (|p| 2 , |q| 2 ) = β p 2 + q 2 2 with β < 0 (focussing NLSE), the hamiltonian energy (1.3) is not bounded from below a priori. In the following analysis, we implicit assume that it exists a constant E c such that E 1 + E c > 0, which is often the case in the study of Bose-Einstein condensate as pointed out by Antoine et al. [START_REF] Antoine | Scalar Auxiliary Variable/Lagrange multiplier based pseudospectral schemes for the dynamics of nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF]. In practice, we compute the term E 1 explicitly and therefore one can adapt the constant E c during the simulation. The system is supplemented by an equation describing the time evolution of r(t). In case of the nonlinear Schrödinger equation (1.1) the continuous SAV model takes the form

     ∂ t p
= -∆q + r(t)g 1 (p, q), ∂ t q = ∆p -r(t)g 2 (p, q), ∂ t r(t) = 1 2 [(g 1 (p, q), ∂ t q) + (g 2 (p, q), ∂ t p)] ,

(

where (•, •) denotes the standard L 2 scalar product and g 1 (p, q) = 1

E 1 [t] + E c δE 1 [t] δq , g 2 (p, q) = 1 E 1 [t] + E c δE 1 [t] δp .
Associated to this SAV model we find the Hamiltonian

H(p, q) = 1 2 Ω |∇p| 2 + |∇q| 2 dx + |r| 2 ,
which is conserved by the SAV model (1.4). In the following, we assume that for i = 1, 2

g i (p, q) ≤ C (|p| + |q|) β + 1 , g i (p, q) ≤ C (|p| + |q|) β + 1 , (1.5)
for some β, β > 0.

Remark 1.1 In this paper, we focus on the Gross-pitaevskii equation under the form (1.1). Even though the choice of the nonlinearity and, therefore, the precise form of E 1 , depends on the structure of the considered Schrödinger equation, we highlight that the SAV scheme is, in its design, general enough to work for a large number of applications. Indeed, as long as there exists a constant E c such that for all times t ≥ 0 it holds that E 1 + E c > 0 we can apply the SAV method. Therefore, modifications such as the effect of dipole-dipole interactions, rotating GPE (see Antoine et al. [START_REF] Antoine | Scalar Auxiliary Variable/Lagrange multiplier based pseudospectral schemes for the dynamics of nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF]), or even time dependent potentials V = V (t, x) can be taken into account. For an extensive overview on applications and generalisations of the nonlinear Schrödinger equation, we refer the interested reader to the review article of Bao and Cai [START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF] and the references therein.

Following the works of Antoine et al. [START_REF] Antoine | Scalar Auxiliary Variable/Lagrange multiplier based pseudospectral schemes for the dynamics of nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF] and Fu et al. [START_REF] Fu | A structure-preserving algorithm for the fractional nonlinear Schrödinger equation based on the sav approach[END_REF], we analyze a fully discrete SAV scheme for the nonlinear Schrödinger equation (1.1) based on a Crank-Nicholson time discretization of the NLS SAV model (1.4) coupled with a pseudo-spectral discretization for the spatial discretization. Energy conservation properties of the SAV method for nonlinear Schrödinger equations were recently derived in [START_REF] Antoine | Scalar Auxiliary Variable/Lagrange multiplier based pseudospectral schemes for the dynamics of nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF][START_REF] Fu | A structure-preserving algorithm for the fractional nonlinear Schrödinger equation based on the sav approach[END_REF] and their convergence was extensively tested numerically. Very recently, Feng et al. [START_REF] Feng | High-order mass-and energy-conserving sav-gauss collocation finite element methods for the nonlinear schrödinger equation[END_REF] use the SAV method to design arbitrary high order space-time finite element scheme for the nonlinear Schrödinger equation. While their method uses a finite element discretization in space, we propose in this work to use a Fourier pseudospectral discretization. The main contribution of this article lies in establishing global error estimates on the fully discrete Fourier-PseudoSpectral Crank-Nicholson NLS SAV scheme (CN-SAV-SP in short). More precisely, we derive weak and strong convergence and prove second order error estimates for the fully discrete scheme. Our theoretical convergence analysis is inspired by the analysis of the SAV method in the context of gradient flows [START_REF] Shen | Convergence and error analysis for the scalar auxiliary variable (SAV) schemes to gradient flows[END_REF]. We underline our convergence results with numerical experiments and compare the SAV scheme with classical splitting methods. Our numerical findings suggest that in certain cases, such as in case of non-linearities involving a non-integer exponent, the SAV scheme preserves its second order energy conservation property while classical splitting methods suffer from sever order reduction. We also conduct numerical experiments showing that the SAV scheme is able to compute correctly ground states of Bose-Einstein condensates.

Outline of the paper. In the first part of the paper, we carry out a fully discrete error analysis of the SAV scheme and establish second order convergence estimates, see Theorem 5.4. Our theoretical convergence results are then numerically underlined in the second part of the paper, see Section 6. We also denote by H m p (Ω) the subset of H m (Ω) that consists of 2π-periodic functions that are in H m (Ω). We denote by L p (0, T ; V ) the Bochner spaces i.e. the spaces with values in Sobolev spaces [START_REF] Adams | Sobolev spaces[END_REF]. The norm in these spaces is defined for all Bochner measurable functions η by

Notations. Let L p (Ω), W m,p (Ω) with H m (Ω) = W m,2 (Ω),
η L p (0,T ;V ) = T 0 η p V dt 1/p , η L ∞ (0,T ;V ) = ess sup t∈(0,T ) η V .
The standard L 2 inner product is denoted by (•, •) Ω and the duality pairing between

H 1 (Ω) = H -1 (Ω) and H 1 (Ω) by < •, • > Ω . The dual space H -1 (Ω) is endowed with the norm φ H -1 (Ω) = sup η∈H 1 (Ω) {< φ, η > Ω , η 2 1 ≤ 1}.
Remark 1.2 Even though our model problem (1.1) is equipped with periodic boundary conditions, our analysis holds for homogeneous Dirichlet or Neumann boundary conditions.

2 Numerical scheme

Time and space discretisation of the SAV model

We use a standard Fourier pseudospectral method [START_REF] Canuto | Spectral methods, Scientific Computation[END_REF][START_REF] Fornberg | A practical guide to pseudospectral methods[END_REF][START_REF] Gottlieb | Theory and applications of spectral methods[END_REF] for the spatial discretization of the SAV model (1.4). We refer the reader to the book of Trefethen [START_REF] Trefethen | Software, Environments, and Tools[END_REF] for details of the implementation of such scheme in MATLAB. We emphasize that our paper presents numerical simulation in dimension d = 1. However, the method can be adapted to higher dimensions. Our convergence and error analysis holds in dimensions 1 ≤ d ≤ 3. Thereby, for the sake of clarity, we here give the details of the space discretization for d = 1. We denote by X N the space spanned by the trigonometric functions up to degree N/2

X N := span{e ikx/L : -N/2 ≤ k ≤ N/2 -1}.
For the time discretisation of the SAV system (1.4) we apply a Crank-Nicholson discretisation with time step τ such that t k = kτ for k ∈ N. At each grid point we thereby approximate the time derivative by

∂ t u(t k+1 , x) ≈ u(t k+1 , x) -u(t k , x) τ .
Let us give the details of the approximation in dimension d = 1, where the domain is defined by Ω = [-π, π] with a mesh size h. In this case the collocation points are x a = 2πa N where a ∈ B with B := {-P, . . . , P -1} if N = 2P is even, {-P, . . . , P } if N = 2P + 1 is odd.

We denote by U k (x a ) the approximation of u(t k , x a ). The Fourier pseudo-spectral discretization is given by

U k (x a ) = p∈B ûk p exp (2iπap/N )
with the Fourier coefficients defined by

ûk p = 1 N b∈B U k (x b ) exp (-2iπbp/N ) .
We approximate the Laplacian by the Fourier differentiation matrix D (2) which for j, l = 0, . . . , N -1 takes the form

D (2) jl =    1 4 (-1) j+1 N + (-1) j+l+1 2 sin 2 (j-l)π N , if j = l -(N -1)(N -2)

12

, otherwise.

However, to avoid the need to use the symmetric matrix D (2) in the previous form and gain in computational time, it can be preferable to use the method proposed in [START_REF] Antoine | Scalar Auxiliary Variable/Lagrange multiplier based pseudospectral schemes for the dynamics of nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF] that uses the fact that the previous differentiation matrix is diagonal in Fourier space. Therefore, inverting this matrix has a very low cost. In our work, since we use the Hamiltonian system (1.2) to analyze the properties of the SAV scheme we will use the previously defined differentiation matrix D (2) .

For the N collocation points x a , we define the interpolation operation I N by

(I N u)(x) = p∈B ũp e 2iπxp/N .
We have the following interpolation error (see Section 5.8.1 in [START_REF] Canuto | Spectral methods, Scientific Computation[END_REF]):

Lemma 2.1 (Interpolation error) For any u ∈ C(0, T ; H m p (Ω) with d ≤ 3, we have

I N u -u H l p (Ω)) ≤ CN l-m |u| m , 0 ≤ l ≤ m, I N ∂ t u -∂ t u H l p (Ω)) ≤ CN l-m |∂ t u| m , 0 ≤ l ≤ m.

The fully discrete SAV scheme

Applying the time discretization described in the previous section, for k = 0 → N T , the semi-discrete model of (1.4) reads

       p k+1 -p k τ = -∆q k+1/2 + r k+1/2 gk+1/2 1 , q k+1 -q k τ = ∆p k+1/2 -r k+1/2 gk+1/2 2 , r k+1 -r k = 1 2 gk+1/2 1 , q k+1 -q k + gk+1/2 2 , p k+1 -p k , (2.1)
where φ k+1/2 = φ k+1 + φ k /2 and gk+1/2 i is a second order extrapolation of g i at time t = t k+1/2 . Denoting by capital letters the vectors of unknowns P k , Q k that are approximations at each collocation nodes of the continuous (in space) unknowns p k , q k , the fully discrete space-time scheme then takes the form

       P k+1 -P k τ = -D (2) Q k+1/2 + R k+1/2 Gk+1/2 1 , Q k+1 -Q k τ = D (2) P k+1/2 -R k+1/2 Gk+1/2 2 , r k+1 -r k = 1 2 Gk+1/2 1 , Q k+1 -Q k + Gk+1/2 2 , P k+1 -P k , (2.2)
where G1 , G2 are the vectors associated to the functions g1 and g2 .

Let us now present two algorithms for the efficient solution of the fully discrete SAV system (2.2). The two methods are equivalent and reduce the problem to the solving of two linear systems that involves only real variables.

Algorithm 1

The algorithm below was originally proposed for solving the fully discrete SAV system arising in gradient flows [START_REF] Shen | The scalar auxiliary variable (SAV) approach for gradient flows[END_REF].

Let us give the procedure on how to solve the system (2.2). First, we need to replace r k+1 in the first two equations using the third equation. This yields that

   P k+1 -P k = -τ D (2) (Q k+1 +Q k ) 2 + τ r k + 1 4 Gk+1/2 1 , Q k+1 -Q k + Gk+1/2 2 , P k+1 -P k Gk+1/2 1 , Q k+1 -Q k = τ D (2) (P k+1 -P k ) 2 -τ r k + 1 4 Gk+1/2 1 , Q k+1 -Q k + Gk+1/2 2 , P k+1 -P k Gk+1/2 2 .
Next we set

Z k = P k Q k , Gk+1/2 = Gk+1/2 2 Gk+1/2 1 , Bk+1/2 = - Gk+1/2 1 Gk+1/2 2 .
This allows us to rewrite the system into a matrix form

AZ k+1 + τ 4 Gk+1/2 , Z k+1 Bk+1/2 = C k , (2.3)
where

A = I τ 2 D (2) -τ 2 D (2) I , and
C k = I -τ 2 D (2) τ 2 D (2) I Z k -τ r k Bk+1/2 + τ 4 Gk+1/2 , Z k Bk+1/2 ,
with I the identity matrix. Multiplying (2.3) by A -1 and taking the discrete inner product with Gk+1/2 , we finally obtain

Gk+1/2 , Z k+1 = Gk+1/2 , A -1 C k 1 + τ 4 Gk+1/2 , A -1 Bk+1/2 . (2.4)
Then, knowing Gk+1/2 , Z k+1 , Z k+1 is computed using (2.3) and r k+1 is calculated from the third equation of (2.2). Therefore, solving the fully discrete SAV model (2.2) reduces to solving the linear system constituted by the equations (2.4) and (2.3).

Algorithm 2

Below we describe a second algorithm recently proposed in [START_REF] Antoine | Scalar Auxiliary Variable/Lagrange multiplier based pseudospectral schemes for the dynamics of nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF] for the numerical solution of the fully discrete NLS SAV scheme (2.2). Rewriting the scheme in its matrix form we have

Z k+1 -Z k τ = -JZ k+1/2 -r k+1/2 Bk+1/2 , r k+1 -r k = 1 2 Gk+1/2 1 , Q k+1 -Q k + Gk+1/2 2 , P k+1 -P k , (2.5)
where

J = 0 D (2)
-D (2) 0 .

Using the decomposition

Z k+1/2 = Z k+1/2 1 + r k+1/2 Z k+1/2 2 , (2.6)
and adding 2 τ Z k on both sides of the first equation of (2.5), we furthermore obtain that

2 τ Z k+1/2 + r k+1/2 Z k+1/2 2 = 2 τ Z k -J Z k+1/2 1 + r k+1/2 Z k+1/2 2 -r k+1/2 Bk+1/2 . (2.7)
Applying the same decomposition to the second equation of (2.5) and adding 2r k on both sides, we get

2r k+1/2 = 2r k + Gk+1/2 , Z k+1/2 1 + r k+1/2 Z k+1/2 2 -Z k . (2.8)
Hence, denoting by I the identity matrix, we first solve the equation (2.7) using the system

2 τ I + J Z k+1/2 1 = 2 τ Z k 2 τ I + J Z k+1/2 2 = -Bk+1/2 . (2.9)
Then we compute r k+1/2 by solving equation (2.8) which yields that

r k+1/2 = 2r k + Gk+1/2 , Z k+1/2 1 -Z k 2 -Gk+1/2 , Z k+1/2 2 .
From the decomposition (2.6) we get Z k+1/2 from which we compute Z k+1 and r k+1 .

Remark 2.2 Since in practice the computation and storage of the invert of a non-diagonal matrix has to be avoided, Algorithm 2 is a preferable choice. Indeed, the main step in Algorithm 2 lies in solving two decoupled linear equations (2.9). To do so, standard tools of linear systems can be applied such as matrix-free preconditioned Krylov solvers. We refer to Appendix C in [START_REF] Shen | Algorithms, analysis and applications[END_REF] for a description of iterative solvers of linear system and preconditioning.

Even though the previous remark already highlights the main advantage of Algorithm 2, we emphasize that the inversion of the main matrix in Algorithm 1 can be carried out efficiently.

Remark 2.3 Referring to [START_REF] Shen | The scalar auxiliary variable (SAV) approach for gradient flows[END_REF], we remark that the inversion of the matrix A in the first algorithm and the matrix 2 τ + J in the second Algorithm can be carried out efficiently using the Sherman-Morrison-Woodbury formula [START_REF] Golub | Matrix computations[END_REF]

(A + U V T ) -1 = A -1 -A -1 U (I + V T A -1 U ) -1 V T A -1 ,
where A is a n × n and U, V are n × k matrices, and I is the k × k identity matrix.

Remark 2.4

Referring to [START_REF] Fu | A structure-preserving algorithm for the fractional nonlinear Schrödinger equation based on the sav approach[END_REF], a fast solver for solving the linear system (2.3) (in Algorithm 1) and (2.9) (in Algorithm 2) exists. It uses the fact that the differentiation matrix D (2) can be decomposed into D (2) = F -1 ΛF where F and F -1 are the corresponding matrices for the discrete Fourier transformation and Λ is a diagonal matrix with eigenvalues of D (2) as its entries. Therefore, the matrix A from Equation (2.3) admits the decomposition

A = F -1 M F, with M = I τ 2 Λ -τ 2 Λ I .
We note that a similar decomposition exists for the matrix [ 2 τ I + J] in equation (2.9). Thanks to the above decomposition, the inverse of the matrix A can be computed explicitly in an efficient manner since

A -1 = F -1 M -1 F, and M -1 = M T (I + τ 2 4 Λ 2 ) -1 0 0 (I + τ 2 4 Λ 2 ) -1 ,
where

(I + τ 2 4 Λ 2
) is a diagonal matrix, such that its inverse is fast to compute.

Conservation properties and inequalities

In this section we outline the conserved quantities of the SAV method. It is well known that due to its design the SAV scheme preserves a modified version of the underlying Hamiltonian. In addition, to the conservation of energy, there is a wide variety of properties in the continuous equation which is feasible to preserve also on the numerical (discrete) level, we refer to Bao and Cai [START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF] as well as Antoine et al. [START_REF] Antoine | Computational methods for the dynamics of the nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF]: i) time-reversibility or symmetry, i.e. the system is unchanged when τ → -τ , ii) gauge-invariance, i.e. if the potential V is changed such that V → V + α with α a real constant then the density |u| 2 remains unchanged, iii) conservation of mass, i.e. u(t) L 2 (Ω) = u(0) L 2 (Ω) , and the Hamiltonian energy, i.e. H(t) = H(0), iv) preservation of the dispersion relation

ω(k) = |k| 2 2 + f (|A| 2 ),
for the plane wave solutions u(t, x) = Ae ik•x-ωt . Proving analytically that the SAV scheme for the NLS equation meets the points ii) and iv) (over long time scales) is up to our knowledge not possible with current techniques. However, the other points can be verified for a large number of nonlinearities. Here, we briefly recall the proofs of the conservation properties and refer to [START_REF] Fu | A structure-preserving algorithm for the fractional nonlinear Schrödinger equation based on the sav approach[END_REF][START_REF] Antoine | Scalar Auxiliary Variable/Lagrange multiplier based pseudospectral schemes for the dynamics of nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF], where they have been first set in context of nonlinear Schrödinger equations and our Theorems 3.1 and 3.2 are found by a combination of the results from [START_REF] Fu | A structure-preserving algorithm for the fractional nonlinear Schrödinger equation based on the sav approach[END_REF] and [START_REF] Antoine | Scalar Auxiliary Variable/Lagrange multiplier based pseudospectral schemes for the dynamics of nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF].

Theorem 3.1 (Conservation of the modified discrete energy) The scheme (2.1) is associated to the discrete modified Hamiltonian

Hk+1 = 1 2 Q k+1 2 1 + P k+1 2 1 + r k+1 2 , (3.1)
and conserves the modified Hamiltonian energy through time i.e.

Hk+1 = Hk . (3.2)
Proof. Taking the inner product with Q k+1 -Q k for the first equation of (2.2) and for the second with -P k+1 -P k , then summing the results we get

0 = 1 2 Q k+1 2 1 -Q k 2 1 + P k+1 2 1 -P k 2 1 +r k+1/2 Gk+1/2 1 , Q k+1 -Q k + Gk+1/2 2 , P k+1 -P k , where |•| 1 = ∇• 0 is the H 1 -seminorm.
Then, multiplying the third equation of (2.2) by 2R k+1/2 and using the result in the previous equation, we obtain

0 = 1 2 Q k+1 2 1 -Q k 2 1 + P k+1 2 1 -P k 2 1 + r k+1 2 -r k 2 ,
from which we can conclude both (3.1) and (3.2).

The SAV scheme also preserves the mass up to an error of order O(τ 3 ), where the latter error is introduced by the second-order extrapolation. Theorem 3.2 (Conservation of the L 2 norm) The scheme (1.4) conserves the L 2 norm of the solution up to an order O(τ 3 ) i.e.

U k+1 2 0 = U k 2 0 + O(τ 3 ), (3.3)
with

U k = P k + iQ k .
Proof. Taking the inner product of first equation of (2.2) with 2P k+1/2 , the second equation with 2Q k+1/2 , and summing the two we get

1 τ P k+1 2 0 -P k 2 0 + Q k+1 2 0 -Q k 2 0 = 2r k+1/2 - Gk+1/2 2 , Q k+1/2 + Gk+1/2 1 , P k+1/2 . Since Gk+1/2 i is a second-order approximation of G k+1/2 i
, we can write

1 τ U k+1 2 0 -U k 2 0 = 2r k+1/2 -G k+1/2 2 , Q k+1/2 + G k+1/2 1 , P k+1/2 + O(τ 2 ).
Then, we find that

G k+1/2 1 = 1 E 1 P k+1/2 , Q k+1/2 + E c ∂E 1 P k+1/2 , Q k+1/2 ∂Q k+1/2 = V (x)Q k+1/2 + f P k+1/2 2 , Q k+1/2 2 Q k+1/2 ,
and

G k+1/2 2 = 1 E 1 P k+1/2 , Q k+1/2 + E c ∂E 1 P k+1/2 , Q k+1/2 ∂P k+1/2 = V (x)P k+1/2 + f P k+1/2 2 , Q k+1/2 2 P k+1/2 ,
from which we easily obtain

-G k+1/2 2 , Q k+1/2 + G k+1/2 1 , P k+1/2 = 0.
Consequently, we obtain (3.3).

To derive H 2 -bound for the solution of the SAV scheme, we use the following proposition. The proof of this technical result can be found in Lemma 2.3 in [START_REF] Shen | Convergence and error analysis for the scalar auxiliary variable (SAV) schemes to gradient flows[END_REF].

Proposition 3.3 (Bound for ∇G k+1/2 i 0
) Assume that the functions g i (i=1,2) satisfy (1.5) and let U 1 ≤ M for some constant M > 0. Then there exists 0 ≤ σ < 1 such that

∇G k+1/2 i ≤ C(M ) 1 + ∇∆P k+1/2 2σ 0 + ∇∆Q k+1/2 2σ 0 . (3.4)
We have the following result on the H 2 -norm of P k+1/2 and Q k+1/2 . Proposition 3.4 (H 2 -bound on the numerical solution) The solution {P k+1 , Q k+1 } of (2.1) satisfies

max k=1,...,N T -1 ∆P k+1 2 0 + ∆Q k+1 2 0 ≤ CT + ∆P 0 2 0 + ∆Q 0 2 0 . (3.5)
Proof.

First, we multiply the first equation of (2.1) by ∆ 2 (Q k+1/2 ), the second equation by ∆ 2 (P k+1/2 ) and integrate over Ω. Then, by summing the two, we obtain, after integration by parts, that

∇∆Q k+1/2 2 0 + ∇∆P k+1/2 2 0 = r k+1/2 ∇ Gk+ 1 2 1 , ∇∆Q k+1/2 ) + r k+1/2 ∇ Gk+ 1 2 2 , ∇∆P k+1/2 .
From the conservation of the modified Hamiltonian (3.1)-(3.2) and assuming a finite initial Hamiltonian, we have

r k+1/2 ∇ Gk+ 1 2 1 , ∇∆Q k+1/2 + r k+1/2 ∇ Gk+ 1 2 2 , ∇∆P k+1/2 ≤ C 2 ∇ Gk+ 1 2 1 2 0 + ∇∆Q k+1/2 2 0 + ∇ Gk+ 1 2 2 2 0 + ∇∆P k+1/2 2 0 .
Then, from the result of Proposition 3.3, for any > 0, we have

∇ Gk+ 1 2 1 2 0 + ∇ Gk+ 1 2 2 2 0 ≤ ∇∆Q k+1/2 2 0 + ∇∆P k+1/2 2 0 + C( ).
Therefore, combining the two previous inequalities, we obtain

∇∆Q k+1/2 2 0 + ∇∆P k+1/2 2 0 ≤ C. (3.6)
Secondly, by multiplying the first equation of (2.1) with ∆ 2 (P k+1/2 ), the second equation with ∆ 2 (Q k+1/2 ), integrating over Ω, and summing the two, we obtain after integration by parts that

∆P k+1 2 0 -∆P k 2 0 + ∆Q k+1 2 0 -∆Q k 2 0 = τ r k+ 1 2 ∇ Gk+ 1 2 1 , ∇∆Q k+1/2 -τ r k+ 1 2 ∇ Gk+ 1 2 2 , ∇∆P k+1/2 .
Then, combining the result of Proposition 3.3 and the inequality (3.6), we have

∆P k+1 2 0 -∆P k 2 0 + ∆Q k+1 2 0 -∆Q k 2 0 ≤ τ C,
and summing from k = 0 → N T , we obtain (3.5).

Remark 3.5 From the fact that H 2 (Ω) ⊆ L ∞ (Ω) for d ≤ 3, we can conclude from the previous proposition that for k = 1, . . . N T -1,

P k+1 L ∞ + Q k+1 L ∞ ≤ C. (3.7)
Next, we present the stability inequality that will be useful in the convergence analysis.

P k+1 2 0 + Q k+1 2 0 + τ 2 N T -1 k=0 P k+1 -P k τ 2 0 + Q k+1 -Q k τ 2 0 ≤ C(τ, H 0 , N T ). (3.8)
Proof.

Multiplying the first equation with 2τ P k+1 , integrating over Ω and using 2(a -b)a = a 2 -b 2 + (a -b) 2 , we obtain

P k+1 2 0 + τ 2 P k+1 -P k τ 2 0 -P k 2 0 = -2τ ∇Q k+1/2 , ∇P k+1 + 2τ r k+1/2 Gk+1/2 1 , P k+1 .
Using the Cauchy-Schwartz inequality and (3.1)-(3.2), we obtain

-2τ ∇Q k+1/2 , ∇P k+1 ≤ 2τ ∇Q k+1/2 0 ∇P k+1 0 ≤ 4τ H 0 .
Then, from the conservation of the Hamiltonian (3.1)-(3.2), and the conservation of the L 2 -norm of the solution (3.3), we obtain using the Cauchy-Schwartz inequality

r k+1/2 Gk+1/2 1 , P k+1 ≤ C Gk+1/2 1 0 P k+1 0 ≤ C G k+1/2 1 0 + O(τ 2 ).
Since from Proposition 3.4 and (3.7), we have that G k+1/2 i 0 ≤ C with i = 1, 2, for a large number of nonlinearities. Therefore, combining the previous inequalities for the right-hand side of (3), we obtain

P k+1 2 0 + τ 2 P k+1 -P k τ 2 0 ≤ Cτ + P k 2 0 .
The same can be found for the second equation by repeating the same calculations. Summing for k = 0 → N T -1, we find (3.8).

4 Convergence analysis

Notations

To study the convergence of the scheme, we introduce the following notation: For k = 0, . . . , N T -1 we set

U (t, x) := t -t k τ U k+1 + t k+1 -t τ U k , t ∈ (t k , t k+1 ],
and

∂U ∂t := U k+1 -U k τ t ∈ (t k , t k+1 ].
We also define

U + := U k+1 , U -:= U k , and
U -U + = (t -t k+1 ) ∂U ∂t , U -U -= (t -t k ) ∂U ∂t t ∈ (t k , t k+1 ], k ≥ 0.
In addition, we take analogous definitions for P and Q: For k = 0, . . . , K T -1 we set

P (t, x) := t -t k τ P k+1 + t k+1 -t τ P k , t ∈ (t k , t k+1 ],
∂P ∂t := P k+1 -P k τ t ∈ (t k , t k+1 ],

P + := P k+1 , P -:= P k , and

P -P + = (t -t k+1 ) ∂P ∂t , and
P -P -= (t -t k ) ∂P ∂t t ∈ (t k , t k+1 ], k ≥ 0.

Convergence theorem

Now we are in the position to establish time convergence for the semi discrete SAV scheme (2.1).

Theorem 4.1 (Convergence) Let {p, q} be a pair of functions such that

   p(t, x) ∈ L 2 [0, T ]; H 1 (Ω) H 1 [0, T ]; H 1 (Ω) q(t, x) ∈ L 2 [0, T ]; H 1 (Ω) H 1 [0, T ]; H 1 (Ω) .
Then for τ → 0 we can extract a subsequence of solutions of (2.1), such that

P, P ± → p strongly in L 2 [0, T ]; L 2 (Ω) , (4.1)
Q, Q ± → q strongly in L 2 [0, T ]; L 2 (Ω) , (4.2)
P, P ± p weakly in L 2 [0, T ]; H 1 (Ω) , (4.3)
Q, Q ± q weakly in L 2 [0, T ]; H 1 (Ω) , (4.4) ∂P ∂t ∂p ∂t weakly in L 2 [0, T ]; H 1 (Ω) , (4.5
)

∂Q ∂t ∂q ∂t weakly in L 2 [0, T ]; H 1 (Ω) , (4.6)
r k+1 r(t) = E 1 [t] + E c weak-star in L ∞ (0, T ) . (4.7)
The limit {p, q} satisfies the nonlinear Schrödinger model (1.2) in the following weak sense

   T 0 ∂p ∂t , η dt = T 0 Ω ∇q∇η + V (x)q + ∂F (|p| 2 ,|q| 2 ) ∂q η dx dt T 0 ∂q ∂t , η dt = T 0 Ω -∇p∇η -V (x)p + ∂F (|p| 2 ,|q| 2 ) ∂p η dx dt, (4.8)
for all η ∈ L 2 [0, T ]; H 1 (Ω) .

Proof.

Step 1: Weak and strong convergences. First, the weak convergences (4.3), (4.4), (4.5) and (4.6) follow from the assumption that the initial Hamiltonian energy is bounded and the stability inequality (3.8).

Then, the weak-star convergence (4.7) also holds true by the conservation of the modified Hamiltonian and the boundedness of the initial state.

From the compact embedding H 1 (Ω) ⊂ L 2 (Ω) ≡ L 2 (Ω) , we can apply the Lions-Aubin Lemma [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF] to find both convergences (4.1) and (4.2).

Step 2: Limit system. Let us work on the first equation of the discrete system. We use a test function η ∈ L 2 [0, T ]; H 1 (Ω) and analyze the convergence of the terms separately. First, from the weak convergence (4.4), we have

T 0 Ω ∇ Q + + Q - 2 ∇η dx dt → T 0 Ω ∇q∇η dx dt.
Secondly, from the fact that Gk+1/2 1 is a second-order approximation of G k+1/2 1

, we have

T 0 Ω r k+1/2 Gk+1/2 1 η dx dt = T 0 Ω r k+1/2 G k+1/2 1 η dx dt + T 0 Ω r k+1/2 O(τ 2 )η dx dt.
From the inequality

E 1 (U k+1/2 ) -E 1 (u) ≤ C U k+1/2 2 L 1 (Ω)
, and the fact that P ± , Q ± p, q weak-star in L ∞ 0, T ; H 1 (Ω) , which follows from the conservation of both the Hamiltonian energy and the L 2 norm, we have

E 1 (U k+1/2 ) E 1 (u) weak-star in L ∞ (0, T ) .
The same holds true for δp using similar arguments. Then, using also the strong convergences (4.1) and (4.2), together with the weak-star convergence (4.7), we obtain

T 0 Ω r k+1/2 G k+1/2 1 η dx dt + T 0 Ω r k+1/2 O(τ 2 )η dx dt → T 0 Ω r(t)g 1 (t)η dx dt.
Finally, for any η ∈ H 1 ([0, T ]; H 1 (Ω)), by integration by parts we have Gathering the previous convergences, we have

(p(T ), η(T )) -(p(0), η(0)) - T 0 p, ∂η ∂t dt = T 0 Ω ∇q∇η + V (x)q + ∂F (|p| 2 , |q| 2 ) ∂q η dx dt. Since ∇q + V (x)q + ∂F (|p| 2 ,|q| 2 )
∂q ∈ L 2 (Ω) which follow from the conservation of the Hamiltonian energy, we know that p ∈ H 1 [0, T ]; H -1 (Ω) . Finally, we find the first equation of the limit system (4.8) and the same arguments can be applied to the second equation. This yields the result.

Error analysis

In this section we analyse the difference between the exact and modified Hamiltonian, and establish a bound on

H[p(t k ), q(t k )] -H[P k , Q k ] .
In addition we prove second-order convergence of the fully discrete SAV scheme (2.2) approximating the solution of the nonlinear Schrödinger equation (1.1). We introduce the following notation to study the error

e k u = θ k u + ρ k u , (5.1)
where

θ k u = U k -(I N u)(t k , x), ρ k u = (I N u)(t k , x) -u(t k , x
). For our convergence result we assume that the solution u of (1.1) is sufficiently smooth satisfying

∂ ttt u L ∞ (0,T ;H 1 (Ω)) + u L ∞ (0,T ;H 2 (Ω)) ≤ C.
(5.2)

We define the different truncation errors by

T k+ 1 2 u = u k+1 -u k ∆t -∂ t u(t k+ 1 2 ), T k+ 1 2 u = u k+ 1 2 -u(t k+ 1 2 ) = u k+1 + u k 2 -u(t k+ 1 2 ).
We commence with two important lemma that will be useful in the global error analysis.

Lemma 5.1 (Boundedness of nonlinear functions) If (p, q) is a solution of (1.4) satisfying (5.2), we have for i = 1, 2

|g i (p, q)| , ∂g i ∂p , ∂g i ∂q , ∂ 2 g i ∂p∂q , ∂ 2 g i ∂p 2 , ∂ 2 g i ∂q 2 ≤ C.
Proof. This result is found by a combination of the fact that u ∈ L ∞ (0, T ; H 2 (Ω)), Remark 3.7, and assumption (1.5).

Remark 5.2 From Lemma 5.1, and the hypothesis (5.2), we know that

|∂ ttt r| ≤ C ∂ ttt p 2 0 + ∂ ttt q 2 0 .
We have the following Lemma on the norm of the truncation errors (see Lemma 4.7 in [START_REF] Wang | Error estimates for second-order SAV finite element method to phase field crystal model[END_REF] for example).

Lemma 5.3 (Truncation errors) For α = -1, 0, 1, 2, we have

T k+ 1 2 ψ 2 H α (Ω) ≤ τ 3 t k+1 t k ∂ ttt ψ(s) 2 H α (Ω) ds, T k+ 1 2 ψ 2 H α (Ω) ≤ τ 3 t k+1 t k ∂ ttt ψ(s) 2 H α (Ω) ds,
Theorem 5.4 (Error analysis) Assume that the solution of (1.2) satisfies (5.2) with initial condition u 0 ∈ H 3 (Ω). Then the discrete solution {P k+1 , Q k+1 } of the fully discrete SAV scheme (2.2) satisfies the error estimate

1 2 ∇e k+1 q 2 0 + 1 2 ∇e k+1 p 2 0 + e k+1 r 2 ≤ C exp [1 -Cτ ] -1 t k+1 τ 4 + N -4 ,
where the constant C depends on the smoothness of the solution (5.2).

Proof.

Step 1: Error equations. We begin by evaluating the model (1.4) at time

t k+1/2      ∂ t p(t k+1/2 ) = -∆q(t k+1/2 ) + r(t k+1/2 )g 1 (t k+1/2 ), ∂ t q(t k+1/2 ) = ∆p(t k+1/2 ) -r(t k+1/2 )g 2 (t k+1/2 ), dr dt (t k+1/2 ) = 1 2 g 1 (t k+1/2 ), ∂ t q(t k+1/2 ) + g 2 (t k+1/2 ), ∂ t p(t k+1/2 ) .
Subtracting the above equations from (2.2) yields

               e k+1 p -e k p τ + T k+1/2 p = -∆ e k+1/2 q + T k+1/2 q + R k+1/2 Gk+1/2 1 -r(t k+1/2 )g 1 (t k+1/2 ), e k+1 q -e k q τ + T k+1/2 q = ∆ e k+1/2 p + T k+1/2 p -R k+1/2 Gk+1/2 2 + r(t k+1/2 )g 2 (t k+1/2 ), e k+1 r -e k r τ + T k+1/2 r = 1 2 Gk+1/2 1 , Q k+1 -Q k τ + Gk+1/2 2 , P k+1 -P k τ -g 1 (t k+1/2 ), ∂ t q(t k+1/2 ) -g 2 (t k+1/2 ), ∂ t p(t k+1/2
) .

(5.3)

We introduce the error e k+1/2 g,1

= Gk+1/2 1 -g 1 (t k+1/2 ).
The rightmost terms of the two first equations of (5.3) can be replaced by

R k+1/2 Gk+1/2 1 -r(t k+1/2 )g 1 (t k+1/2 ) = Gk+1/2 1 e k+1/2 r + T k+1/2 r + r(t k+1/2 )e k+1/2 g,1 , (5.4)
and

R k+1/2 Gk+1/2 2 -r(t k+1/2 )g 2 (t k+1/2 ) = Gk+1/2 2 e k+1/2 r + T k+1/2 r + r(t k+1/2 )e k+1/2 g,2
.

(5.5)

Similarly, we have

1 2 Gk+1/2 1 , Q k+1 -Q k τ -g 1 (t k+1/2 ), ∂ t q(t k+1/2 ) = 1 2 Gk+1/2 1 , e k+1
q -e k q τ + T k+1/2 q + e k+1/2 g,1

, ∂ t q(t k+1/2 ) ,

(5.6) and 1 2

Gk+1/2

2 , P k+1 -P k τ -g 2 (t k+1/2 ), ∂ t p(t k+1/2 ) = 1 2 Gk+1/2 2 , e k+1 p -e k p τ + T k+1/2 p + e k+1/2 g,2
, ∂ t p(t k+1/2 ) .

(5.7)

Plugging (5.4), (5.5),(5.6), and (5.7) into (5.3), we thus obtain , ∂ t p(t k+1/2 ) .

                   e k+1 p -e k p τ + T k+1/2 p = -∆ e k+1/2 q + T k+1/2 q + Gk+1/2 1 e k+1/2 r + T k+1/2 r + r(t k+1/2 )e k+1/2 g,1 , e k+1 q -e k q τ + T k+1/2 q = ∆ e k+1/2 p + T k+1/2 p - Gk+1/2 2 e k+1/2 r + T k+1/2 r -r(t k+1/2 )e k+1/2 g,2 , e k+1 r -e k r τ + T k+1/2 r = 1 2 Gk+1/2 1 , e k+1 q -e k q τ + T k+1/2 q + e k+1/2 g,1 , ∂ t q(t k+1/2 ) + Gk+1/2
Using the decomposition of the error (5.1), we furthermore obtain , ∂ t p(t k+1/2 ) .

                                 θ k+1 p -θ k p τ + ∆ θ k+1 q +θ k q 2 = - ρ k+1 p -ρ k p τ -∆ ρ k+1/2 q + T k+1/2 q + Gk+1/2 1 e k+1/2 r + T k+1/2 r +r(t k+1/2 )e k+1/2 g,1 -T k+1/2 p , θ k+1 q -θ k q τ -∆ θ k+1 p +θ k p 2 = - ρ k+1 q -ρ k q τ + ∆ ρ k+1/2 p + T k+1/2 p - Gk+1/2 2 e k+1/2 r + T k+1/2 r -r(t k+1/2 )e k+1/2 g,2 -T k+1/2 q , e k+1 r -e k r τ + T k+1/2 r = 1 2 Gk+1/2 1 , e k+1 q -e k q τ + T k+1/2 q + e k+1/2 g,1 , ∂ t q(t k+1/2 ) + Gk+1/2
(5.8)

Step 2. Error estimate formula. We use the following notations to make the results more compact

D 1 τ θ k+1 p = θ k+1 p -θ k p τ , D 1 θ k+1 p = θ k+1 p -θ k p .
Taking the inner product of the first equation of the system (5.8) with -D 1 θ k+1 q and the second with D 1 θ k+1 p , and summing the results, we also have,

1 2 D 1 ∇θ k+1 q 2 0 + 1 2 D 1 ∇θ k+1 p 2 0 = D 1 τ ρ k+1 p , D 1 θ k+1 q -D 1 τ ρ k+1 q , D 1 θ k+1 p -∇ρ k+1/2 q , ∇D 1 θ k+1 q -∇ρ k+1/2 p , ∇D 1 θ k+1 p -∇T k+1/2 q , ∇D 1 θ k+1 q -∇T K+1/2 p , ∇D 1 θ k+1 p - Gk+1/2 1 e k+1/2 r + T k+1/2 r , D 1 θ k+1 q - Gk+1/2 2 e k+1/2 r + T k+1/2 r , D 1 θ k+1 p -r(t k+1/2 )e k+1/2 g,1 -T k+1/2 p , D 1 θ k+1 q -r(t k+1/2 )e k+1/2 g,2 + T k+1/2 p , D 1 θ k+1 p .
(5.9)

Multiplying the third equation of (5.8) by 2τ e k+1/2 r , we have

D 1 e k+1 r 2 + 2τ T k+1/2 r e k+1/2 r -τ e k+1/2 r Gk+1/2 1 , D 1 τ ρ k+1 q + T k+1/2 q + Gk+1/2 2 , D 1 τ ρ k+1 p + T k+1/2 p + e k+1/2 g,1 , ∂ t q(t k+1/2 ) + e k+1/2 g,2 , ∂ t p(t k+1/2 ) = e k+1/2 r Gk+1/2 1 , D 1 θ k+1 q + Gk+1/2 2 , D 1 θ k+1 p .
(5.10) Using (5.10) in (5.9), we have

1 2 D 1 ∇θ k+1 q 2 0 + 1 2 D 1 ∇θ k+1 p 2 0 + D 1 e k+1 r 2 = D 1 τ ρ k+1 p , D 1 θ k+1 q -D 1 τ ρ k+1 q , D 1 θ k+1 p -∇ ρ k+1/2 q + T k+1/2 q , ∇D 1 θ k+1 q -∇ ρ k+1/2 p + T k+1/2 p , ∇D 1 θ k+1 p -T k+1/2 r Gk+1/2 1 , D 1 θ k+1 q + Gk+1/2 2 , D 1 θ k+1 p -2τ T k+1/2 r e k+1/2 r + τ e k+1/2 r Gk+1/2 1 , D 1 τ ρ k+1 q + T k+1/2 q + Gk+1/2 2 , D 1 τ ρ k+1 p + T k+1/2 p + e k+1/2 g,1 , ∂ t q(t k+1/2 ) + e k+1/2 g,2 , ∂ t p(t k+1/2 ) -r(t k+1/2 )e k+1/2 g,1 -T k+1/2 p , D 1 θ k+1 q -r(t k+1/2 )e k+1/2 g,2 + T k+1/2 p , D 1 θ k+1 p .
(5.11)

Step 3. Inequalities for the terms on the right-hand side of (5.11). Now, we bound the right-hand side of (5.11). Using Lemma 2.1, Lemma 5.3 and Young's inequality we have

D 1 τ ρ k+1 p , D 1 θ k+1 q ≤ 4 D 1 τ ρ k+1 p 2 0 + 1 16 D 1 θ k+1 q 2 0 ≤ CN -6 + 1 16 D 1 θ k+1 q 2 0 , -D 1 τ ρ k+1 q , D 1 θ k+1 p ≤ 4 D 1 τ ρ k+1 q 2 0 + 1 16 D 1 θ k+1 p 2 0 ≤ CN -6 + 1 16 D 1 θ k+1 p 2 0 .
Then, from Theorem (3.1), we have

-∇ρ k+1/2 q , ∇D 1 θ k+1 q -∇ρ k+1/2 p , ∇D 1 θ k+1 p ≤ ∇ρ k+1/2 q 2 0 ∇D 1 θ k+1 q 2 0 + ∇ρ k+1/2 p 2 0 ∇D 1 θ k+1 p 2 0 ≤ CN -4 , and -∇T k+1/2 q , ∇D 1 θ k+1 q -∇T K+1/2 p , ∇D 1 θ k+1 p ≤ ∇T k+1/2 q 2 0 ∇D 1 θ k+1 q 2 0 + ∇T K+1/2 p 2 0 ∇D 1 θ k+1 p 2 0 ≤ Cτ 4 ,
For the rest of the terms on the right-hand side of (5.11), we use Lemma 5.3, and Proposition 3.6 together with Lemma 5.1, and Remark 5.2, to obtain

-T k+1/2 r Gk+1/2 1 , D 1 θ k+1 q + Gk+1/2 2 , D 1 θ k+1 p ≤ 4 T k+1/2 r 2 Gk+1/2 1 2 0 + Gk+1/2 1 2 0 + 1 16 D 1 θ k+1 q 2 0 + D 1 θ k+1 p 2 0 ≤ Cτ 4 + 1 16 D 1 θ k+1 q 2 0 + D 1 θ k+1 p 2 0 , (5.12)
, ∂ t q(t k+1/2 ) + e k+1/2 g,2 , ∂ t p(t k+1/2 ) ≤ τ 2 ∂ t q(t k+1/2 ) 2 0 e k+1/2 g,1 2
0

+ e k+1 r 2 + e k r 2 + τ 2 ∂ t p(t k+1/2 ) 2 0 e k+1/2 g,2 2
+ T k+1/2 q 2 0 + 1 16 D 1 θ k+1 q 2 0 + D 1 θ k+1 p 2 0 ≤ C e k+1/2 g,1 2 0 + e k+1/2 g,2 2 0 + 2τ 4 + 1 16 D 1 θ k+1 q 2 0 + D 1 θ k+1 p 2 0 . (5.13)
Step 4. Estimating the terms in the inequalities (5.12)-(5.13). First, we aim to emiminate the terms D 1 θ k+1 p 2 0 and D 1 θ k+1 q 2 0 in the above inequalities. Taking the inner product of the first equation of (5.8) with 2τ θ k+1 p , we obtain

D 1 τ θ k+1 p , 2τ θ k+1 p = 2τ ∇θ k+1/2 q , ∇θ k+1 p -2τ D 1 τ ρ k+1 p , θ k+1 p + 2τ ∇ ρ k+1/2 q + T k+1/2 q , ∇θ k+1 p + 2τ e k+1/2 r + T k+1/2 r Gk+1/2 1 , θ k+1 p + 2τ r(t k+1/2 )e k+1/2 g,1 -T k+1/2 p , θ k+1 p .
Knowing that . Using the notations S(p, q) = E 1 (p, q) + C, and

D 1 τ θ k+1 p , 2τ
N 1 (p, q) = δ δq E 1 (p, q), N 2 (p, q) = δ δp E 1 (p, q)
we have that e k+1/2 g,1

= G 1 (P k+1/2 , Q k+1/2 ) -g 1 (p(t k+1/2 ), q(t k+1/2 ))

= N 1 (P k+1/2 , Q k+1/2 ) S(P k+1/2 , Q k+1/2 ) -N 1 (p(t k+1/2 ), q(t k+1/2 )) S(p(t k+1/2 ), q(t k+1/2 ))

= N 1 (P k+1/2 , Q k+1/2 ) S(P k+1/2 , Q k+1/2 ) -N 1 (P k+1/2 , Q k+1/2 ) S(p(t k+1/2 ), q(t k+1/2 ))

+ N 1 (P k+1/2 , Q k+1/2 ) S(p(t k+1/2 ), q(t k+1/2 ))

-N 1 (p(t k+1/2 ), q(t k+1/2 )) S(p(t k+1/2 ), q(t k+1/2 ))

= N 1 (P k+1/2 , Q k+1/2 ) E 1 (p(t k+1/2 ), q(t k+1/2 )) -E 1 (P k+1/2 , Q k+1/2 ) S(P k+1/2 , Q k+1/2 )S(p(t k+1/2 ), q(t k+1/2 )) S(P k+1/2 , Q k+1/2 ) + S(p(t k+1/2 ), q(t k+1/2 ))

+ N 1 (P k+1/2 , Q k+1/2 ) -N 1 (p(t k+1/2
), q(t k+1/2 )) S(p(t k+1/2 ), q(t k+1/2 )) .

From the smoothness assumption (5.2), Lemma 5.1, and Remark 3.5, we have e k+1/2 g,1

2 0 ≤ C P k+1/2 -p(t k+1/2 ) 2 0 + Q k+1/2 -q(t k+1/2 ) 2 0 .
Then, using the notation (5.1) and Lemma 5. Therefore, by the use of Gronwall's Lemma, we can conclude that

1 2 ∇θ k+1 q 2 0 + 1 2 ∇θ k+1 p 2 0 + e k+1 r 2 ≤ C exp [1 -Cτ ] -1 t k+1 τ 4 + N -4 .

Numerical experiments

In this section we numerically confirm our theoretical convergence result given in Theorem 5.4 and illustrate the long time energy conservation of the SAV method. In the following, the numerical results have been obtained from Algorithm 2. However, we want to emphasize that Algorithm 1 leads to the same results and has a comparable computational cost for d = 1.

Our numerical findings suggest the favorable energy preservation of the SAV method compared to classical splitting methods in certain applications such as for non-linearities with non-integer exponents which arise for instance in context of optical dark and power law solitons with surface plasmonic interactions [START_REF] Crutcher | Derivation of the Effective Nonlinear Schrödinger Equations for Dark and Power Law Spatial Plasmon-Polariton Solitons Using Nano Self-Focusing[END_REF]. For the comparison we use the classical first order Lie and second order Strang splitting which are known for their near energy preservation over long times, see, e.g., [START_REF] Faou | Geometric numerical integration and Schrödinger equations[END_REF].

In the numerical examples we plot the deviation of the exact Hamiltonian and the modified Hamiltonian H, i.e., e H = H(t k ) -Hk , the error between the exact Hamiltonian and the discrete nonmodified Hamiltonian e H = H(u(t k )) -H(U k ) , as well as the L 2 error e u = U k -u(t k ) L 2 (Ω) . We choose the potential V = 0 in the Schrödinger equation (1.1).

First test case: cubic nonlinearity

In a first example we consider the nonlinear Schrödinger equation (1.1) with a cubic nonlinearity i.e.

f (|u| 2 ) = β |u| 2
on the spatial domain Ω = [-32, 32]. In Figure 1 we choose a mesh size h = 1/32 and approximate the soliton solution [START_REF] Bao | Numerical methods and comparison for computing dark and bright solitons in the nonlinear Schrödinger equation[END_REF][START_REF] Antoine | Scalar Auxiliary Variable/Lagrange multiplier based pseudospectral schemes for the dynamics of nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF]

u(x, t) = a √ -β sech (a(x -vt)) exp ivx -0.5(v 2 -a 2 )t ,
with the parameters a = 1, β = -1 and v = 1 up to T = 10. Figure 1 numerically confirms the secondorder convergence of the SAV method. The numerical findings also suggest that the error constant of the Strang splitting method is slightly better than the one of the SAV method in this example. In Figure 2 we simulate the solitary wave

u(t, x) = √ 2e it cosh(x)
on the domain [-π 0.11 , π 0.11 ] with N = 256 collocation points and time step size τ = 0.01. We illustrate the evolution of the errors e H , e u and e H over long times, i.e., up to T = 1000. Our numerical findings confirm the conservation of the modified Hamiltonian by the SAV method, see Figure 2. We also observe that the SAV method preserves well the exact energy and L 2 norm over long times. Even though, the error e H of the Strang splitting seems favorable in this example, we have to stress that the modified Hamiltonian is closer to the value of the real Hamiltonian (see error e H on Figure 2).

Second test case: cubic nonlinearity with non-smooth initial condition

In this example we analyse the error behaviour of the SAV scheme in case of non-smooth initial data. For this purpose we solve the NLS equation with cubic nonlinearity with initial data of various Figure 3 shows the convergence behaviour of the SAV scheme, and the two splitting methods for the initial data of different regularity. We find that if α < 3, the SAV method does not maintain its second order convergence rate and for α = 2, the SAV scheme reduces to first order. Decreasing the regularity of the initial condition even more, the convergence worsens and becomes less than order 1. A similar order reduction is observed for the splitting schemes, however, for the latter the error starts to oscillate for α < 3. Again, the error e H is favorable for the Strang splitting for all α. However, the modified Hamiltonian is closer to the real Hamiltonian (see Figure 4 for α = 2 3 ).

Third test case: non-integer exponent

In this example we consider the periodic nonlinear Schrödinger equation (1.1) with nonlinearities with non-integer exponents ( [START_REF] Crutcher | Derivation of the Effective Nonlinear Schrödinger Equations for Dark and Power Law Spatial Plasmon-Polariton Solitons Using Nano Self-Focusing[END_REF]) where the Hamiltonian takes the form

f (|u| 2 ) = β |u| 4/γ , γ > 0
H(u) = Ω 1 2 |∇u| 2 + β γ (4 + 2γ) |u| 4 γ +2 dx.
We carry out simulations for various exponents γ = 2, 8/3, 4, 8 up to time T = 10 with smooth initial value u(0, x) = sin(x) ∈ C ∞ ([-π, π]).

The error e H for different exponents γ is plotted in Figure 5. Our numerical findings suggest that as γ increases the splitting methods suffer from sever order reduction. This loss of convergence of splitting methods was also observed in [START_REF] Ostermann | Low regularity exponential-type integrators for semilinear Schrödinger equations[END_REF]. The SAV method, on the other hand, retains its second order energy convergence for non-integer exponents.

Computing ground states

We use the SAV scheme to simulate ground states of Bose-Einstein condensates, however, unlike the work of Antoine et al. [START_REF] Antoine | Scalar Auxiliary Variable/Lagrange multiplier based pseudospectral schemes for the dynamics of nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF], we propose here to use a different strategy. Indeed, in [START_REF] Antoine | Scalar Auxiliary Variable/Lagrange multiplier based pseudospectral schemes for the dynamics of nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF], the authors use the SAV scheme presented in Section 2 and observe the capacity of the scheme to preserve the initial mass and Hamiltonian for various strengths of the nonlinearity. In the present work, we propose to use a different method and compare the numerical results with reference methods that are designed to simulate the stationary states of the NLS equation for large nonlinearities.

To do so, we reformulate the problem into the solving of a gradient flow equation to compute these stationary states: this method is known has the gradient flow with discrete normalization method [START_REF] Bao | Ground states and dynamics of multicomponent Bose-Einstein condensates[END_REF][START_REF] Bao | Computing the ground state solution of Bose-Einstein condensates by a normalized gradient flow[END_REF].

The SAV scheme is well adapted to this formulation since its original purpose was the simulation of the gradient flow equations. Details of the reformulation and the adaptation of the SAV scheme to the case can be found in Appendix A. We here present numerical results obtained choosing d = 1, V (x) = x 2 /2, β = 400, and u 0 (x) = exp(-x 2 /2) π 1/4 / u 0 0 . We validate our results with the Backward Euler PseudoSpectral (BEPS in short) scheme implemented in the GPELab code [START_REF] Antoine | Gpelab, a matlab toolbox to solve gross-pitaevskii equations i: Computation of stationary solutions[END_REF][START_REF] Gpelab | a matlab toolbox to solve gross-pitaevskii equations ii: Dynamics and stochastic simulations[END_REF].

We denote by E(u) the energy associated to the renormalized system and Ẽ(u) its modified SAV energy (see Appendix A for details).

Figure 6 (left) compares the stationary states obtained with the two schemes for h = 1/8 in space. We clearly see that they both reach the same steady state. Figure 6 (right) depicts the evolution of the energy during the simulation. We observe that the SAV scheme preserves the monotonic decay of the energy. The steady state reached at the end of the simulation has an energy Ẽ(φ) ≈ 22.90. However, using the solution φ g obtained for the SAV scheme and computing the "real" energy, we obtain E(φ) ≈ 21.36 which is the value obtained with the BEPS scheme. We numerically evaluate the order of convergence in space of the SAV scheme for the simulation of ground states. We choose our reference solution to be the result of a simulation with h = 1/32. Then, we vary h from 1/2 to 1/16. Figure 7 shows that the scheme remains second order convergent in space as predicted by our error analysis.

A Gradient flow with discrete normalization for computing ground state

A common method to compute stationary states of the NLS equation (1.1) with a cubic nonlinearity is to write u(t, x) = φ(x) exp(-iµt),

where µ is defined as the chemical potential of the condensate µ(φ) = For the proof of the existence of such state and other mathematical properties we refer the reader to [START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF].

In the following, we adapt the Scalar Auxiliary Variable method to compute the stationary solutions of Equation (1.1). Therefore, endowing the equation with the normalization constraint, and using the projected gradient method [6], the complete system reads ∂ t φ = 1 2 ∆φ -V (x)u -β |u| 2 u(t, x), φ 2 L 2 (Ω) = 1.

Our SAV scheme can be easily adapted to this case, leading to the discrete system

       φ + -φ k τ = 1 2 D (2) φ k+1/2 -r k+1/2 Gk+1/2 , r k+1 -r k = 1 2 Gk+1/2 , φ + -φ k φ k+1 = φ + φ + 2 L (Ω) ,
with φ k+1/2 = φ + +φ k 2 , Gk+1/2 a second order approximation of δE 1 [t k+1/2 ] δφ k+1/2 . We precise that the associated modified SAV energy is

Ẽ(φ) = Ω 1 2
|∇φ| dx + r(t) 2 < +∞.

  where 1 ≤ p ≤ +∞ and m ∈ N, denote the standard Lebesgue and Sobolev spaces equipped with the corresponding norms || • || m,p , || • || m and semi-norms | • | m,p , | • | m .

Proposition 3 . 6 (

 36 Stability inequality) The solution of (2.1) satisfies the stability inequality max k=0,...,N T -1

  P (T ), η(T )) -(P (0), η(0)) .Hence, from the regularity of η and the convergence (4.1) , we obtain as τ → 0 and ∀η ∈ H 1 ([0, T ]; H 1 (Ω)).

2 + e k r 2 ≤ 5 + τ e k+1 r 2 + τ e k r 2 , 2 + e k r 2 ≤ 2 + e k r 2 , 2 + e k r 2 ≤ 2 + e k r 2 ,

 2252222222222 Cτ Cτ N -6 + τ 4 + e k+1 r Cτ N -6 + τ 4 + e k+1 r

2 q- 4 ++ τ 4 + C ∇θ k+1 p 2 0 ,D 1 + CN - 4 ++ CN - 4 +

 24420144 Let us bound the terms on the right-hand side of(5.14). Using Lemma 2.1 we find that 2τ ∇θ k+1/Cτ 4 + ∇θ k+1 p where we have used the Poincaré inequality to obtain the last inequality. Plugging the previous inequalities into (5.14), we obtain Cτ 4 + C e the inner product of the second equation of (5.8) with 2τ θ k+1 q and repeating the same steps as before, we obtainD Cτ 4 + C e

+ τ 3 + N - 4 ≤ C ∇θ k+1 p 2 0+ τ 3 + N - 4 .+ τ 3 + N - 4 .+ D 1 e k+1 r 2 ≤ τ C ∇θ k+1 p 2 0 2 + e k r 2 + C τ 4 +

 342343422224 Step 6. Discrete Gronwall Lemma. The above two estimates together with (5.15) and (5.16) imply N -4 .

Figure 1 :

 1 Figure 1: Error e u (left) and e H (right) versus step size τ at time T = 10.

Figure 2 :

 2 Figure 2: Left Figure: error e u (left) through time. Right Figure: e H (blue, red, yellow) and e H (purple) through time.

Figure 3 :

 3 Figure 3: Error e H versus step size τ for the nonlinear Schrödinger equation starting from different initial conditions in H α (α = 2 3 top-left, α = 2 top-right, α = 3 bottom left and α = 5 bottom-right). The dotted lines represent order τ (green) and τ 2 (purple), respectively.

Figure 4 :

 4 Figure 4: Error e H (blue, red, yellow) through time and e H (purple) through time for initial condition in H 2 3 .

Figure 5 :

 5 Figure 5: Error e H versus step size τ for the nonlinear Schrödinger equation with different non-integer exponents (γ = 2 top-left, γ = 8 3 top-right, γ = 4 bottom left and γ = 8 bottom-right). The dotted lines represent order τ (green) and τ 2 (purple), respectively.

Figure 6 :

 6 Figure 6: (Left) Comparison of stationary solutions of the NLS equation with a large cubic nonlinearity obtained with the SAV scheme (blue) and the BEPS scheme from GPELab (red). (Right) Evolution of the energy Ẽ(u) for the SAV scheme during the simulation.

Ω 1 2 |∇φ| + β |φ| 4 +

 24 V (x) |φ| 2 dx.Therefore, using the previous reformulation in Equation (1.1), we obtainµφ(x) = -1 2 ∆φ(x) + β |φ(x)| 2 φ(x) + V (x)φ(x).

Figure 7 :|φ| 4 +

 74 Figure 7: Error e u versus grid size h for the simulation of ground states with a large cubic nonlinearity.

† The author has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 740623) § The author has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 850941)

Both Algorithm 1 and Algorithm 2 from Section 2 can be applied to compute the solution of the SAV system. Furthermore, using the same calculation as in Section 3, we can easily prove that the scheme dissipates the energy and preserves the normalization constraint.