Exponential integrators for nonlinear Schrödinger equations with white noise dispersion

This article deals with the numerical integration in time of the nonlinear Schrödinger equation with power law nonlinearity and random dispersion. We introduce a new explicit exponential integrator for this purpose that integrates the noisy part of the equation exactly. We prove that this scheme is of mean-square order 1 and we draw consequences of this fact. We compare our exponential integrator with several other numerical methods from the literature. We finally propose a second exponential integrator, which is implicit and symmetric and, in contrast to the first one, preserves the L2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2$$\end{document}-norm of the solution.


Introduction
We consider the time discretisation of the following nonlinear Schrödinger equation with white noise dispersion idu + c u • dβ + |u| 2σ u dt = 0 where the unknown u = u(x, t), with t ≥ 0 and x ∈ R d , is a complex valued random process, u = d j=1 denotes the Laplacian in R d , c is a real number, σ is a positive real number, and β = β(t) is a real valued standard Brownian motion. This stochastic partial differential equation is understood in the Stratonovich sense, using the • symbol for the Stratonovich product. The existence of a unique global square integrable solution to (1.1) was shown in [14] for σ < 2/d and in [15] for d = 1 and σ = 2, see also [3]. The existence and uniqueness of solutions to the one-dimensional cubic case of the above problem was also studied in [26]. Furthermore, as for the deterministic Schrödinger equation, the L 2 -norm, or mass, of the solution to (1.1) is a conserved quantity. This is not the case for the total energy of the problem.
We now review the literature on the numerical analysis of the nonlinear Schrödinger equation with white noise dispersion (1.1). The early work [18] studies the stability with respect to random dispersive fluctuations of the cubic Schrödinger equation. Furthermore, numerical experiments using a split step Fourier method are presented. The paper [26] presents a Lie-Trotter splitting integrator for the above problem (1.1). The mean-square order of convergence of this explicit numerical method is proven to be at least 1/2 for a truncated Lipschitz nonlinearity [26,Sect. 5 and 6]. Furthermore, [26] conjectures that this splitting scheme should have order one, and supports this conjecture numerically. An analysis of asymptotic preserving properties of the Lie-Trotter splitting is carried out in [16] for a more general nonlinear dispersive equation. Very recently, the authors of [3] studied an implicit Crank-Nicolson scheme for the time integration of (1.1). They show that this scheme preserves the L 2 -norm and has order one of convergence in probability. Finally, the recent preprint [13] examine the multi-symplectic structure of the problem and derive a multi-symplectic integrator which converges with order one in probability.
In the present publication, we will consider exponential integrators for an efficient time discretisation of the nonlinear stochastic Schrödinger equation (1.1). Exponential integrators for the time integration of deterministic semi-linear problems of the formẏ = Ly + N (y), are nowadays widely used and studied, as witnessed by the recent review [22]. Applications of such numerical schemes to the deterministic (nonlinear) Schrödinger equation can be found in, for example, [4][5][6][7][8][9][10]17,21] and references therein. Furthermore, these numerical methods were investigated for stochastic parabolic partial differential equations in, for example, [23][24][25], more recently for the stochastic wave equations in [2,11,12,27], where they are termed stochastic trigonometric methods, and lately to stochastic Schrödinger equations driven by Ito noise in [1].
The main result of this paper is a mean-square convergence result for an explicit and easy to implement exponential integrator for the time discretisation of (1.1). Indeed, we will show in Sect. 3 convergence of mean-square order one for this scheme as well as convergence in probability. Note that the proofs of the results presented here use similar techniques as the one used in [3].
In order to show the above convergence result, we begin the exposition by introducing some notations and recalling useful results in Sect. 2. After that, we present our explicit exponential integrator for the numerical approximation of the above stochastic Schrödinger equation and analyse its convergence in Sect. 3. Various numerical experiments illustrating the main properties of the proposed numerical scheme will be presented in Sect. 4. In the last section, we discuss the preservation of the mass, or L 2 -norm, by symmetric exponential integrators.

Notation and useful results
We denote the classical Lebesgue space of complex functions by L 2 := L 2 (R d , C), endowed with its real vector space structure, and with the scalar product For s ∈ N, we further denote by H s := H s (R d , C) the Sobolev space of functions in L 2 such that their s first derivatives are in L 2 . The Fourier transform of a tempered distribution v is denoted by F (v) or v. With this, H s is the Sobolev space of tempered distributions v such that (1 + |ζ | 2 ) s/2 v ∈ L 2 .
Next, we consider a filtered probability space ( , F , P, {F t } t≥0 ) generated by a one-dimensional standard Brownian motion β = β(t).
With the above definitions in hand, we can write the mild formulation of the stochastic nonlinear Schrödinger equation (1.1) (with the constant c = 1 for ease of presentation) [3,14,26] with the random propagator S(t, r ) expressed in the Fourier variables as for t ≥ r ≥ 0, ζ ∈ R d and v a tempered distribution. We finally collect some results that we will use in the error analysis presented in Sect. 3: • The random propagator S(t, r ) is an isometry in H s for any s, see for example [3].
• There is a constant C such that, for t ≥ 0, h ∈ (0, 1) and r ∈ (t, t + h) and for any F t -measurable function v ∈ L 2 ( , H s+4 ), one has the bounds (see [ • Without much loss of generality, we will truncate the nonlinearity in (1.1) in Sect. 3 and thus recall the following estimates from [3]. Let f be a function from H s to H s , which sends H s+2 to itself and H s+4 to itself, with f (0) = 0, twice continuously differentiable on those spaces, with bounded derivatives of order 1 and 2. Consider Then, there exists a constant C, which depends on f , such that (see [3,Equations (2.30) and (2.44)]) with h, r, t as in the above point (provided that the right-hand side is finite).

Exponential integrator and mean-square error analysis
This section presents an explicit time integrator for (1.1), and further states and proves a mean-square convergence result for this numerical method. As a by-product result, we also obtain convergence in probability of the exponential integrator.

Presentation of the exponential integrator
Let T > 0 be a fixed time horizon and an integer N ≥ 1. We define the step size of the numerical method by h = T /N and denote the discrete times by t n = nh, for n = 0, . . . , N . Looking at the mild solution (2.1) of the problem (1.1) on the interval [t n , t n+1 ], and discretising the integral (by freezing the integrand at the left-end point of this interval), one can iteratively define the following explicit exponential integrator We thus obtain a finite sequence of numerical approximations u n ≈ u(t n ) of the exact solution to the problem (1.1) at the discrete times t n = nh.

Truncated Schrödinger equation
As in [3], we introduce a cut-off function in order to cope with the nonlinear part of the stochastic partial differential equation Observe that, for s > d/2 and σ ∈ N , for a fixed k ∈ N * , f k is a bounded Lipschitz function from H s to H s which sends H s+2 to H s+2 and H s+4 to H s+4 . It is twice differentiable on these spaces, with bounded and continuous derivatives of order 1 and 2. Thus one has a unique global adapted solution u k to the truncated problem in L ∞ ( , C ([0, T ], H s )) if the initial value u 0 ∈ H s , see [3]. Note that, with the assumptions above, u k ∈ L ∞ ( , C ([0, T ], H s+2 )) as soon as u 0 ∈ H s+2 , and u k ∈ L ∞ ( , C ([0, T ], H s+4 )) as soon as u 0 ∈ H s+4 .
The global solution u k ∈ L ∞ ( , C ([0, T ], H s )) to the truncated problem solves (3.2) and the exponential integrator takes the form Note that this method looks like the composition of two methods: the first is the explicit Euler equation applied to the differential equation u = f k (u), the second is the exact solution of the linear stochastic Schrödinger equation.

Main result and convergence analysis
This subsection states and proves the main result of this paper on the mean-square convergence of the exponential integrator applied to the nonlinear Schrödinger equation with white noise dispersion (1.1).
for the discrete times t n = nh. Here, the constant C does not depend on n, h with nh ≤ T but may depend on k, T, Proof For ease of presentation, we will ignore the index k referring to the cut-off in the notations of the numerical and exact solutions as well as in the nonlinear function f k . But we keep in mind that the constants below may depend on this index. We denote by C such a constant, providing it does not depend on n ∈ N nor on h ∈ (0, 1) such that nh ≤ T . In order to later apply a discrete Gronwall-type argument, we first look at the error between the exact and numerical solutions The so-called mean-square error thus reads with the H s norm · 2 s = Re(·, ·) s = · 2 H s . We now proceed with the estimations of the above quantities. Before estimating the mixed terms in (3.4), we start with the first four terms. By isometry of the random propagator S(t, r ), one gets For the second term, we again use the isometry property of the free random propagator and further the fact that the function f is Lipschitz. This gives us Using the isometry property of S(t, r ), Cauchy-Schwarz's inequality, the estimate (2.2) from Sect. 2, and the fact that the exact solution is bounded, we can bound the third term by In order to estimate the fourth term, we use the isometry property of the free random propagator, Cauchy-Schwarz's inequality, the fact that f is Lipschitz, the estimate (2.4) on the time-variations of the exact solution, and the fact that the exact solution is bounded which is recalled in Sect. 2. We then obtain Next, we go on with deriving bounds for the mixed terms present in (3.4). Using Cauchy-Schwarz's inequality, the above bounds for the moments of I 2 and I 3 , and the fact that for all real numbers a, b, we have ab ≤ 1 2 (a 2 + b 2 ), we obtain the bound Similarly, one has and The term containing I 1 and I 2 can be estimated using Cauchy-Schwarz's inequality and the fact that the function f is Lipschitz: The last two terms |E[(I 1 , I 3 ) s ]| and |E[(I 1 , I 4 ) s ]| demand more work. For the first one, we use the isometry of the random propagator S(t, r ) and Cauchy-Schwarz's inequality to get: We next apply the law of total expectation and again Cauchy-Schwarz's inequality in order to get Finally, using (2.3) and the fact that the exact solution is bounded, one obtains the bound In order to estimate the last term |E[(I 1 , I 4 ) s ]|, we use the mild formulation and then Cauchy-Schwarz's inequality and a Taylor expansions It thus remains to bound the above two terms. In order to start to estimate the first term, we insert the mild solution and obtain For the second term, we again use Cauchy-Schwarz's inequality with the fact that S, D f and D 2 f are bounded and the bound (2.5): Altogether, we arrive at Collecting all the above bounds, the mean-square error (3.4) can thus be estimated by and an application of a discrete Gronwall lemma gives the final bound which concludes the proof of the theorem.
Using the above mean-square convergence result and similar arguments as in [19,26] or [3], one can also show that the exponential method (3.1) has order of convergence one in probability.

Numerical experiments
This section presents various numerical experiments for the nonlinear Schrödinger equation with white noise dispersion (1.1). We will use the following numerical schemes: 1. The explicit exponential integrator (3.1); 2. The Lie-Trotter splitting from [26]. Here, Y (h)u * denotes the value at time h of flow associated to the problem i ∂u ∂t + |u| 2σ u = 0 with initial datum u * ; 3. The Strang splitting where again Y (h) is defined as above; 4. The implicit Crank-Nicolson scheme from [3]. Here, we have set u n+1/2 = 1 2 (u n + u n+1 ), We will consider the stochastic partial differential equation (1.1) on the one and two dimensional torus with periodic boundary conditions. The spatial discretisation is done by a pseudospectral method with M Fourier modes in 1d, and M 2 Fourier modes in 2d.

Numerical experiments in 1d
This subsection presents convergence plots for the above mentioned numerical methods applied to the nonlinear Schrödinger equation with white noise dispersion (1.1); space-time evolution plots; experiments illustrating the influence of the power nonlinearity σ supporting a conjecture proposed in [3]; and finally illustrations of the preservation of the L 2 -norm along numerical solutions.

Convergence plots
In order to illustrate the mean-square convergence of the exponential integrator (3.3), we consider problem (1.1) on the interval [0, 2π ] with parameters c = σ = 1. The initial datum is u 0 (x) = e −(x−π) 2 for x ∈ [0, 2π ]. Furthermore, M = 2 8 Fourier modes are used for the spatial discretisation. The mean-square errors Fig. 1a for various values of the time step h = 2 − for = 6, . . . , 17. Here, we simulate the exact solution u(x, t) with the exponential method, with a small time step h exact = 2 −17 . The expected values are approximated by computing averages over M s = 5000 samples. We computed the estimate for the largest standard errors to be 5.78 × 10 −4 for the Crank-Nicolson scheme and around 10 −6 for the other numerical schemes. These estimates are far from optimal but we observed that using a larger number of samples (M s = 10,000) does not improve significantly the behaviour of the convergence plots. This is also the case for the other convergence plots presented below. In Fig. 1a, we observe convergence of order 1 for the exponential integrator. This is in agreement with Theorem 3.1. The orders of convergence for the splitting schemes and for the Crank-Nicolson scheme are also seen to be 1.
Note that the explicit exponential method (3.3), as well as the Lie and Strang splitting methods (4.1)-(4.2) take full advantage of the exact integration of the stochastic linear part of the Schrödinger equation when one uses periodic boundary conditions and hence the spectral properties of the Laplace operator are exactly known. In con-trast, the Crank-Nicolson scheme (4.3) does not integrate the stochastic part of the equation exactly. One can even argue that the error in the identity which is the cornerstone of the error analysis of the linear part of the Crank-Nicolson scheme applied to a Schrödinger equation, is fully under control in the deterministic case (when X = − h 2 ξ 2 , and ξ is the Fourier variable), while it can be much higher in the stochastic case (when X = −c W n 2 ξ 2 ) even for small time steps. Once again, such an error corresponding directly to the stochastic part of the PDE is not present in the three other schemes (3.1), (4.1) and (4.2). This explains, in 1d as well as in higher dimension (see next section for numerical examples in dimension 2), the relatively poor behaviour of the Crank-Nicolson scheme in this situation (see Fig. 1a) even if c is of order 1. On the other, as observed in Fig. 1b, the good convergence behaviour of the Crank-Nicolson is recovered for c = 0.25 (smaller noise intensity parameter). The other parameters are the same as in the previous numerical experiments. Figure 2 shows the evolution of |u n | 2 along one sample of the numerical solution obtained by the exponential integrator (3.1) for the above problem with the discretisation parameters h = 2 −14 and M = 2 8 . This illustrates, in the case σ = 1, the interplay and the balance between the random dispersion and the nonlinearity. In contrast, the qualitative behaviour is different for higher values of σ .

Evolution plots
The article [3] conjectures that the power nonlinearity σ = 4 is the critical case for (1.1) in dimension one. We now present numerical experiments supporting this conjecture. Problem (1.1) with the initial value u 0 (x) = 2.3 · e −14(x−π) 2 is integrated over the time interval [0, 0.05] with discretisation parameters h = 2 −12 and M = 2 14 . Figure 3 shows the space-time evolution for the power nonlinearity σ = 3.9 and σ = 4.  Blow-up of the solution can be observed in the critical case σ = 4 thus numerically confirming the conjecture from [3].
We now perform another numerical test to support the criticality of the exponent σ = 4. We consider the numerical integration of (1.1) with c = 1.0 for several values of σ , namely 3.0, 3.9 and 4.1, starting from the initial datum u 0 (x) = 10e −10x 2 for x ∈ [−20π, 20π ], over the time interval [0, 7.10 −8 ]. We use the Lie-Trotter method (4.1) and the explicit exponential scheme (3.1) to integrate the problem. We use M = 2 19 Fourier modes in space. We run several samples for each method and we plot the numerical H 1 -norm as a function of time. We observe numerical blow up in finite time for some samples. We count the number of samples that blow up in finite time, and the evolution of this number when one refines the time step. Numerical results are presented in Figs. 4 and 5. We observe that, for each method, for σ = 3.0 and σ = 3.9, reducing the time step leads to a smaller proportion of solutions that blow up in finite time. In contrast when σ = 4.1, for each method, the percentage of solutions that blow up in finite time does not decrease when one reduces the time step. This illustrates numerically the criticality of the exponent σ = 4.0 in dimension one.
Of course, this is only a rough result and one can think of more sophisticated techniques such as adaptive mesh refinement techniques to have a better understanding of the behaviour of the solution close to the blow-up.

Preservation of the L 2 -norm
It is known that the L 2 -norm of the solution to the SPDE (1.1) remains constant for all times [3]. Figure 6 illustrates the corresponding preservation properties of the above numerical integrators along one sample path. For this numerical experiment, we consider the parameters c = 1, σ = 1, h = 2 −5 , M = 2 8 and the initial value u 0 (x) = e −10(x−π) 2 for x ∈ [0, 2π ]. Exact preservation of the L 2 -norm for the splitting schemes and for the Crank-Nicolson scheme is observed, whereas a small drift is observed for the exponential integrator (3.1). In Sect. 5, we will propose a symmetric exponential integrator that preserves exactly the L 2 -norm.

Numerical experiments in 2d
This subsection presents convergence plots for (1.1) in two dimensions as well as experiments illustrating the influence of the power nonlinearity σ supporting a conjecture proposed in [3].

Convergence plots
We illustrate the mean-square convergence of the exponential integrator In this figure, we observe convergence of order 1 for the exponential integrator. This is in agreement with Theorem 3.1.  Time

Evolution plots
Let us now consider the following parameters c = 1, h = 2 −11 , M = 2 7 and the initial value 5 · e −14((x−π/2) 2 +(y−π/2) 2 ) · e i10x for (x, y) grators using ideas from [9]. We thus propose the following symmetric exponential method for the numerical discretisation of nonlinear Schrödinger equation with white noise dispersion (1.1) where N (u) = i|u| 2σ u is the nonlinearity. This numerical method preserves the L 2 -norm as seen in the following proposition. Proof This proof is an adaptation of the proof stating conditions for a Runge-Kutta methods to preserve quadratic invariants, see [20, Section IV.2.1] and further [9]. Let us compute the L 2 -norm of u 1 : using the isometry of the random propagator S(t, r ). We next define Y : . Inserting this quantity in the above equation and using the definition of S(t, r ) yields = u 0 2 since the last term in brackets is zero (this follows from the fact that N * = N (Y ) and the fact that the L 2 -norm is a first integral).

Computational cost
The goal of this section is to compare the computational cost of the explicit exponential method (3.1) and its symmetric version (5.1) introduced in this paper to that of classical numerical methods from the literature, for the integration of the Cauchy problem (1.1). We choose to implement the Crank-Nicolson method (4.3), and the Lie and Strang splitting methods (4.1), resp (4.2). We run the five methods over the time interval [0, 0.01] for the Cauchy problem (1.1) with c = 1, σ = 3 and initial datum We see that all numerical methods but the Crank-Nicolson method actually behave rather similarly for small time steps, since the points we obtain are almost aligned. The fact that the Crank-Nicolson method behaves a rather differently to the other four methods can be explained mainly by the fact that it is the only numerical method in this test that does not integrate the linear stochastic part of the equation exactly (as noticed in the analysis of the convergence plots in Sect. 4.1.1 above). Note that if one takes more complicated space geometries or space discretisations (and not equispaced points on a finite interval that allow using FFT) so that the linear part of the equation can no longer be integrated exactly, then the numerical cost of the other four numerical methods may be of the same order as that of the Crank-Nicolson method. ments show that it outperforms methods that do not integrate exactly the linear part of the equation, such as the Crank-Nicolson method, in terms of size of constant errors, for reasonably large noise intensity (Fig. 1). Furthermore, we used this new scheme in Sect. 4.2 to support a conjecture on the critical power to get blow-up in finite time in the nonlinear Schrödinger equation (1.1). Finally, we proposed another exponential integrator (5.1) which is symmetric and has the same numerical order as the one proposed initially. It however is implicit and hence has higher numerical cost.