On splitting methods for Schrödinger-Poisson and cubic nonlinear Schrödinger equations

Abstract

We give an error analysis of Strang-type splitting integrators for nonlinear Schrödinger equations. For Schrödinger-Poisson equations with an $H^4$-regular solution, a first-order error bound in the $H^1$ norm is shown and used to derive a second-order error bound in the $L_2$ norm. For the cubic Schrödinger equation with an $H^4$-regular solution, first-order convergence in the $H^2$ norm is used to obtain second-order convergence in the $L_2$ norm. Basic tools in the error analysis are Lie-commutator bounds for estimating the local error and $H^m$-conditional stability for error propagation, where $m=1$ for the Schrödinger-Poisson system and $m=2$ for the cubic Schrödinger equation.

1. Introduction

In this paper we give an error analysis of the Strang splitting time integration method applied to nonlinear Schrödinger equations

$$\begin{equation} i\frac{\partial \psi }{\partial t} = -\Delta \psi + V\psi \,,\qquad x\in {\mathbf{R}}^3,\ t>0, {\tag{1.1}} \end{equation}$$

where

$$\begin{equation} V = V[\psi ] =\pm |\psi |^2 {\tag{1.2}} \end{equation}$$

in the case of the cubic nonlinear Schrödinger equation, and

$$\begin{equation} -\Delta V = \pm |\psi |^2 {\tag{1.3}} \end{equation}$$

in the case of the Schrödinger-Poisson equations. The equations are considered with asymptotic boundary conditions $\lim _{|x|\to \infty } \psi (x,t)=0$ and $\lim _{|x|\to \infty } V(x)=0$. The Poisson equation in (1.3) is thus to be interpreted as giving $V$ by the convolution with the fundamental solution of the negative Laplacian,

$$V = V[\psi ]=\mp \Delta ^{-1}|\psi |^2 := \pm \frac{1}{4\pi |x|}\ast |\psi |^2\,.$$

In both cases, the initial data is given as $\psi (x,0)=\psi _0(x)$ for $x\in {\mathbf{R}}^3$.

The cubic nonlinear Schrödinger equation arises as a model equation from several areas of physics; see, e.g., Sulem and Sulem 20. The one-dimensional problem ($x\in {\mathbf{R}}$) is important in fiber optics; see Agrawal 1. Schrödinger-Poisson equations (1.1), (1.3) (also known as the Hartree equation), and generalizations are basic equations in quantum transport; see, e.g., Brezzi, and Markowich 6 and Illner, Zweifel, and Lange 13. The more elaborate Schrödinger-Poisson system considered there has the same mathematical difficulties as (1.1) with (1.3), so we restrict our attention to this simpler set of equations.

In this paper we study the approximation properties of a semi-discretization in time. The numerical integrator we consider is a Strang-type splitting method, yielding approximations $\psi _n$ to $\psi (t_n)$ with $t_n=n\tau$ for a step size $\tau >0$ via

$$\begin{eqnarray} \psi _{n+1/2}^- &=& e^{\frac{i}{2}\tau \Delta }\psi _n, \\ \psi _{n+1/2}^+ &=& e^{-i\tau V[\psi _{n+1/2}^-]}\,\psi _{n+1/2}^-, {\tag{1.4}}\\ \psi _{n+1} &=& e^{\frac{i}{2}\tau \Delta }\psi _{n+1/2}^+ \,. \end{eqnarray}$$

Here, $e^{it\Delta }$ is the solution operator of the free Schrödinger equation, expressed in terms of Fourier transforms as ${\mathcal{F}}^{-1}e^{-it|\xi |^2}{\mathcal{F}}$ and approximately computed by FFT in a Fourier spectral method, whereas the exponential of $V$ acts as a pointwise multiplication operator. Note that $|\psi _{n+1/2}^+|= |\psi _{n+1/2}^-|$ and hence $V[\psi _{n+1/2}^+]=V[\psi _{n+1/2}^-]$. Method (1.4) is therefore explicit and time-reversible. The method is the composition of the exact flows of the differential equations

$$i\frac{\partial \psi }{\partial t} = -\Delta \psi \qquad \text{and}\qquad i\frac{\partial \psi }{\partial t} = V[\psi ]\psi \,.$$

Such splitting methods are widely used; see, e.g., the early references Strang 19 and Hardin and Tappert 11, the study of the split-step Fourier method for the cubic nonlinear Schrödinger equation by Weideman and Herbst 21 and its use in fiber optics as in Agrawal 1, Section 2.4, the use of splitting methods for the time-dependent Kohn-Sham equations (closely related to the above Schrödinger-Poisson equations) in time-dependent density functional theory by Appel and Gross 2, and the papers by Bao, Mauser, and Stimming 4 on the use in the Schrödinger-Poisson-$X\alpha$ model and by Bao, Jaksch, and Markowich 3 on the numerical solution of the Gross-Pitaevskii equation for Bose-Einstein condensation, which is closely related to the cubic nonlinear Schrödinger equation. We further refer to the review of splitting methods by McLachlan and Quispel 18.

To our knowledge, there is as yet no rigorous convergence result in the literature for the splitting method for the cubic nonlinear Schrödinger equation. We mention, however, the work by Besse, Bidégaray, and Descombes 5, where an error analysis is given for globally Lipschitz-continuous nonlinearities, which is not the case with the cubic nonlinearity considered here. For the Schrödinger-Poisson equation, a first-order $L_2$ error bound over a time interval $[0,T]$ with suitably small $T$ for initial data in the Sobolev space $H^2$ has been shown by Fröhlich 8.

Here, we derive error bounds for the Strang splitting over any given finite time interval that are second-order accurate in the $L_2$ norm under the condition of $H^4$ spatial regularity. This is more stringent than the $H^2$ regularity needed for linear Schrödinger equations with a smooth bounded potential 14. The higher regularity requirement for the nonlinear equations considered here is caused by a term $\Delta ^2\psi$ in the double Lie commutator of $i\Delta$ with the nonlinearity, whereas in the linear case there is a cancellation of higher derivatives that leaves only second-order derivatives. It is also interesting to compare with finite-difference time-stepping methods such as the Crank-Nicholson method or the implicit midpoint rule, for which second-order error bounds involve bounds on the third time derivative of the solution, which would require $H^6$-spatial regularity.

We remark that Weideman and Herbst 21 report an instability phenomenon in the Strang splitting for the cubic Schrödinger equation for certain step sizes, caused by resonances between the linear part, which has its spectrum on the imaginary axis, and the nonlinearity. This instability can lead to an exponential error growth that is stronger than in the error propagation by the equation itself, and can thus impair the long-time behaviour of the method. It should be noted, however, that this potential long-time instability is not at odds with the finite-time stability and convergence results given here.

We restrict our attention in this paper to nonlinear Schrödinger equations (1.1) on the whole space ${\mathbf{R}}^3$. Our arguments would apply similarly to problems with periodic boundary conditions and in lower space dimension, and could be extended to nonlinear Schrödinger equations with other power nonlinearities.

We only study semi-discretization in time but we expect that the results extend to various types of full discretization, uniformly in the spatial discretization parameter. What needs to be checked is the discrete version of the Lie commutator bounds established in this paper for the spatially continuous case. Once such bounds are available, the theory extends to the fully discrete case without further ado. The same remark apparently applies to splitting methods for other nonlinear evolution equations such as the KdV equation, where similarly the scheme of proof given here becomes applicable once the necessary Lie bracket bounds are established.

Throughout the paper, $L_2=L_2({\mathbf{R}}^3)$ denotes the Hilbert space of Lebesgue square integrable functions, and $H^k=H^k({\mathbf{R}}^3)$ is the Sobolev space of $L_2$-functions having all generalized derivatives up to order $k$ in $L_2$. We denote the solution of (1.1) at time $t$ by $\psi (t)=\psi (\cdot ,t)$. The $L_2$ norm is preserved along the solution, and we assume it to be of unit norm: $\|\psi (t)\|_{L_2}=\|\psi _0\|_{L_2}=1$.

The paper is organized as follows. In the first part (Sections 2 to 6) we consider the Schrödinger-Poisson equation (1.1), (1.3) and then, in Sections 7 and 8, we extend the results and techniques to the cubic Schrödinger equation. Sections 2 and 7 state the results of this paper. In Section 3 we give some inequalities related to the nonlinearity in the Schrödinger-Poisson equation. In Section 4 we prove the first-order error bound in the $H^1$ Sobolev norm for solutions in $H^3$, and in Section 5 this is used to show the second-order error bound in $L_2$ for $H^4$-regular solutions. Section 6 proves an $H^2$-regularity result of the numerical solution. Finally, Section 8 outlines the modifications in the proofs needed for the cubic Schrödinger equation.

PART A. SCHRÖDINGER-POISSON EQUATIONS

2. Error bounds for solutions in $H^4$: Statement of results

In this section we formulate error bounds in the $H^1$ and $L_2$ norm and state some related results. According to a result by Illner, Zweifel, and Lange 13, the Schrödinger-Poisson equation (1.1), (1.3) has a global strong solution: $\psi _0\in H^2$ implies $\psi (t)\in H^2$ for all $t\ge 0$. The result can be extended to yield $H^k$ regularity of solutions to $H^k$ initial data for any $k\ge 2$ globally in time. We suppose that the solution $\psi (t)$ to the Schrödinger-Poisson equation (1.1), (1.3) is in $H^4$ for $0\le t \le T$, and set

$$m_k = \max _{0\le t \le T} \| \psi (t) \|_{H^k} \qquad \text{for }\ k\le 4.$$

Our main result concerning the error of the Strang-type splitting scheme (1.4) reads as follows.

Theorem 2.1.

Suppose that the exact solution $\psi (t)$ to the Schrödinger-Poisson equation 1.1, 1.3 is in $H^4$ for $0\le t \le T$. Then, the numerical solution $\psi _n$ given by the splitting scheme 1.4 for the Schrödinger-Poisson equation 1.1, 1.3 with step size $\tau >0$ has a first-order error bound in $H^1$ and a second-order error bound in $L_2$,

$$\begin{array}{l} \| \psi _n - \psi (t_n) \|_{H^1} \le C(m_3,T) \,\tau \\[2.84526pt] \| \psi _n - \psi (t_n) \|_{L_2}\, \le C(m_4,T) \,\tau ^2 \end{array} \qquad \text{ for }\ t_n=n\tau \le T\,.$$

The following auxiliary results are of independent interest. We write the step of the splitting scheme (1.4) briefly as

$$\psi _{n+1} = \Phi _\tau (\psi _n)\,.$$

Proposition 2.2 ($H^1$-conditional $L_2$- and $H^1$-stability).

If $\psi ,\phi \in H^1$ with

$$\|\psi \|_{H^1} \le M_1,\quad \ \|\phi \|_{H^1} \le M_1,$$

then

$$\begin{eqnarray*} &&\| \Phi _\tau (\psi )-\Phi _\tau (\phi ) \|_{L_2} \le e^{c_0\tau } \, \|\psi -\phi \|_{L_2}, \\ &&\| \Phi _\tau (\psi )-\Phi _\tau (\phi ) \|_{H^1} \le e^{c_1\tau } \, \|\psi -\phi \|_{H^1}, \end{eqnarray*}$$

where $c_0,c_1$ only depend on $M_1$.

Note that also the $L_2$-stability estimate depends on bounds in $H^1$. The proof of Theorem 2.1 therefore proceeds by first showing the $H^1$ error bound, which, in particular, establishes the required bound of the $H^1$ norm of numerical solutions. We then are in the position to prove the $L_2$ error bound using the $H^1$-conditional $L_2$-stability.

Proposition 2.3 (Local error in $H^1$).

If $\psi _0\in H^3$ with $\|\psi _0\|_{H^3} \le M_3$, then the error after one step of the method 1.4 is bounded in the $H^1$ norm by

$$\| \psi _1 - \psi (\tau ) \|_{H^1} \le C_3\tau ^2\,,$$

where $C_3$ only depends on $M_3$.

Proposition 2.4 (Local error in $L_2$).

If $\psi _0\in H^4$ with $\|\psi _0\|_{H^4} \le M_4$, then the error after one step of the method 1.4 is bounded in the $L_2$ norm by

$$\| \psi _1 - \psi (\tau ) \|_{L_2} \le C_4\tau ^3\,,$$

where $C_4$ only depends on $M_4$.

Proposition 2.5 ($H^2$ regularity of the numerical solution).

If $\psi _0\in H^2$ and

$$\|\psi _n\|_{H^1} \le M_1\qquad \text{ for all $n$ with $n\tau \le T$},$$

then in fact

$$\|\psi _n\|_{H^2} \le e^{c_2n\tau } \|\psi _0\|_{H^2}\qquad \text{ for } n\tau \le T,$$

where $c_2$ only depends on $M_1$.

3. Some inequalities

Hardy’s inequality (e.g., 15, p. 350)

$$\int _{{\mathbf{R}}^3} \frac{|u(y)|^2}{|y|^2}\, dy \le 4 \int _{{\mathbf{R}}^3} |\nabla u(y)|^2\, dy \qquad (u\in H^1)$$

implies some further inequalities that play an important role in the following.

Lemma 3.1.

For $u\in H^1$ and $v,w\in L_2$,

$$\begin{equation} \| \Delta ^{-1}(uv)w\|_{L_2} \le K_0\,\| u \|_{H^1} \,\|v\|_{L_2} \, \| w \|_{L_2}\,, {\tag{3.1}} \end{equation}$$

and for $u,v\in L_2$ and $w\in H^1$,

$$\begin{equation} \| \Delta ^{-1}(uv)w\|_{L_2} \le K_0\,\| u \|_{L_2} \,\|v\|_{L_2} \, \| w \|_{H^1}\,. {\tag{3.2}} \end{equation}$$

Proof.

(a) Inequality (3.1) is essentially Lemma 3.3 of 13. We have

$$\| \Delta ^{-1}(uv)w\|_{L_2} \le \| \Delta ^{-1}(uv) \|_{L_\infty } \, \|w\|_{L_2}$$

and further, using the Cauchy-Schwarz inequality and Hardy’s inequality,

$$\begin{eqnarray*} \| \Delta ^{-1}(uv) \|_{L_\infty } &=& \sup _x \int _{{\mathbf{R}}^3} \frac{u(x-y)v(x-y)}{4\pi |y|}\, dy \\ &\le & \frac{1}{4\pi }\,\sup _x \Bigl ( \int _{{\mathbf{R}}^3} \frac{|u(x-y)|^2}{|y|^2}\, dy \Bigr )^{1/2} \Bigl ( \int _{{\mathbf{R}}^3} {|v(x-y)|^2}\, dy \Bigr )^{1/2} \\ &\le & \frac{1}{2\pi }\, \| u \|_{H^1}\, \| v \|_{L_2}\,. \end{eqnarray*}$$

(b) For the proof of (3.2) we use a duality argument. Using partial integration and the $L_\infty$ bound of part (a), we obtain

$$\begin{eqnarray*} &&\| \Delta ^{-1}(uv)w\|_{L_2} = \sup _{\|\phi \|_{L_2}=1} \int _{{\mathbf{R}}^3} \Delta ^{-1}(uv)w\phi \, dx = \sup _{\|\phi \|_{L_2}=1} \int _{{\mathbf{R}}^3} uv \,\Delta ^{-1}(w\phi )\, dx \\ &&\qquad \le \| uv \|_{L_1} \, \sup _{\|\phi \|_{L_2}=1} \|\Delta ^{-1}(w\phi )\|_{L_\infty } \le \| u \|_{L_2}\, \| v \|_{L_2}\,\frac{1}{2\pi }\, \| w \|_{H^1}\,, \end{eqnarray*}$$

which yields the result with $K_0=\frac{1}{2\pi }$.

■

With the product rule of derivatives, Lemma 3.1 immediately yields the following bounds.

Lemma 3.2.

We have

$$\begin{equation} \| \Delta ^{-1}(uv)w\|_{H^1} \le K_1\left( \| u \|_{H^1} \,\|v\|_{H^1} \, \| w \|_{L_2} +\| u \|_{H^1} \,\|v\|_{L_2} \, \| w \|_{H^1} \right) {\tag{3.3}} \end{equation}$$

for $u,v,w\in H^1$, and

$$\begin{equation} \| \Delta ^{-1}(uv)w\|_{H^2} \le K_2 \sum _{(k,\ell ,m)} \| u \|_{H^k} \,\|v\|_{H^\ell } \, \| w \|_{H^m} {\tag{3.4}} \end{equation}$$

for $u,v,w\in H^2$, where the sum is over all permutations $(k,\ell ,m)$ of $(0,1,2)$.■

For further inequalities concerning $\Delta ^{-1}(uv)w$ we refer to Castella 7 and Illner, Zweifel, and Lange 13.

4. Proof of the first-order error bound in $H^1$

4.1. $H^1$-conditional stability: Proof of Proposition 2.2

(a) Since $e^{i\tau \Delta }$ preserves both the $L_2$ and the $H^1$ norm, we only need to compare $e^{-i\tau V[\psi ]}\psi$ and $e^{-i\tau V[\phi ]}\phi$, which are the solutions at time $\tau$ of the linear initial value problems

$$\begin{eqnarray*} && i\dot{\theta }= V[\psi ]\theta , \quad \ \theta (0)=\psi , \\ && i\dot{\eta }= V[\phi ]\eta , \quad \ \eta (0)=\phi , \end{eqnarray*}$$

with $H^1$ norms of $\psi$ and $\phi$ bounded by $M_1$. We rewrite the difference of the right-hand sides as

$$\begin{eqnarray*} V[\psi ]\theta - V[\phi ]\eta &=& \Delta ^{-1} |\psi |^2 \cdot \theta - \Delta ^{-1} |\phi |^2 \cdot \eta \\ &=& \Delta ^{-1} \bigl ( |\psi |^2 -|\phi |^2 \bigr )\,\theta + \Delta ^{-1} |\phi |^2 \, (\theta -\eta ) \\ &=& \Delta ^{-1} \bigl ( (\psi -\phi )\overline{\psi }\bigr )\,\theta + \Delta ^{-1} \bigl ( \phi (\overline{\psi }-\overline{\phi })\bigr )\,\theta + \Delta ^{-1} (\phi \overline{\phi }) \, (\theta -\eta ). \end{eqnarray*}$$

(b) By Lemma 3.1, we thus obtain

$$\begin{eqnarray*} \| V[\psi ]\theta - V[\phi ]\eta \|_{L_2} &\le & K_0 \,\|\psi -\phi \|_{L_2} \, \|\psi \|_{H^1}\, \|\theta \|_{L_2} + K_0 \,\|\psi -\phi \|_{L_2} \, \|\phi \|_{H^1}\, \|\theta \|_{L_2} \\ && +\ K_0\, \|\phi \|_{H^1} \, \|\phi \|_{L_2} \, \|\theta -\eta \|_{L_2} \end{eqnarray*}$$

and hence, recalling unit $L^2$ norms of $\phi$ and $\theta$,

$$\| \theta (t)-\eta (t) \|_{L_2} \le \|\psi -\phi \|_{L_2} + 2K_0 M_1 t\,\|\psi -\phi \|_{L_2} + \int _0^t K_0 M_1 \, \|\theta (s)-\eta (s)\|_{L_2}\, ds\,,$$

so that by the Gronwall inequality,

$$\| e^{-i\tau V[\psi ]}\psi - e^{-i\tau V[\phi ]}\phi \|_{L_2} = \| \theta (\tau )-\eta (\tau ) \|_{L_2} \le e^{c_0\tau }\, \|\psi -\phi \|_{L_2}$$

where $c_0$ depends on $M_1$.

(c) We proceed in the same way for the $H^1$ estimate, using now Lemma 3.2 (and recalling unit $L_2$ norms) for the estimate

$$\| V[\psi ]\theta - V[\phi ]\eta \|_{H^1} \le 2K_1 \bigl ( \|\phi \|_{H^1} + \|\theta \|_{H^1} \bigr ) \, \|\psi -\phi \|_{H^1} + 2K_1 \|\phi \|_{H^1} \, \|\theta -\eta \|_{H^1}\,.$$

Next we estimate the $H^1$ norm of $\theta (t)$. By Lemma 3.2 and unit $L_2$ norms, we have

$$\| V[\psi ]\theta \|_{H^1} \le K_1 \left( \|\psi \|_{H^1}^2 + \|\psi \|_{H^1} \, \|\theta \|_{H^1} \right)\,,$$

which yields, using the bound $\|\psi \|_{H^1} \le M_1$,

$$\|\theta (t)\|_{H^1} \le (1+ K_1 M_1 t)\, \|\psi \|_{H^1} + \int _0^t M_1\, \|\theta (s)\|_{H^1} \, ds\,.$$

With the Gronwall inequality we thus obtain

$$\begin{equation} \| e^{-itV[\psi ]}\psi \|_{H^1} = \|\theta (t)\|_{H^1} \le e^{a_1 t} \, \|\psi \|_{H^1}\,, {\tag{4.1}} \end{equation}$$

where $a_1$ only depends on $M_1$. With the above estimate of $\| V[\psi ]\theta - V[\phi ]\eta \|_{H^1}$ this gives

$$\begin{eqnarray*} \| \theta (t)-\eta (t) \|_{H^1} &\le & \|\psi -\phi \|_{H^1} + \int _0^t 2K_1 M_1 (1+e^{a_1 s})\,\|\psi -\phi \|_{H^1} \,ds \\ &&\ +\ \int _0^t 2K_1 M_1\, \|\theta (s)-\eta (s)\|_{H^1}\, ds\,. \end{eqnarray*}$$

Once again by the Gronwall inequality, we finally obtain

$$\| e^{-i\tau V[\psi ]}\psi - e^{-i\tau V[\phi ]}\phi \|_{H^1} = \| \theta (\tau )-\eta (\tau ) \|_{H^1} \le e^{c_1 \tau } \, \|\psi -\phi \|_{H^1}\,,$$

where $c_1$ only depends on $M_1$.■

4.2. Commutator bound

We consider the vector fields on dense subsets of $H^1$,

$$\widehat{T}(\psi ) = i\Delta \psi \,,\qquad \widehat{V}(\psi ) = -iV[\psi ]\psi \,,$$

that appear in (1.1). Their Lie commutator

$$\begin{eqnarray*} &&\!\!\!\!\!\! [\widehat{T},\widehat{V}](\psi ) = \widehat{T}'(\psi ) \widehat{V}(\psi ) - \widehat{V}'(\psi ) \widehat{T}(\psi ) \\ &&\!\!\!\!\!\! = i\Delta \bigl ( -i\Delta ^{-1} (\psi \overline{\psi })\psi \bigr ) + i\Delta ^{-1} (i\Delta \psi \,\overline{\psi })\psi + i\Delta ^{-1} (\psi \overline{i\Delta \psi })\psi + i\Delta ^{-1} (\psi \overline{\psi })i\Delta \psi \\ &&\!\!\!\!\!\! = 2 \Delta ^{-1} (\nabla \psi \cdot \nabla \overline{\psi })\, \psi +2 \Delta ^{-1} (\nabla \psi \,\overline{\psi })\cdot \nabla \psi +2 \Delta ^{-1} (\psi \,\nabla \overline{\psi })\cdot \nabla \psi +2 \Delta ^{-1} (\psi \overline{\Delta \psi })\psi \end{eqnarray*}$$

plays an essential role in the error estimate.

Lemma 4.1.

The commutator is bounded in $H^1$ by

$$\| [\widehat{T},\widehat{V}](\psi ) \|_{H^1} \le C\, \|\psi \|_{H^2} ^3 + C\,\|\psi \|_{H^1}^2\,\|\psi \|_{H^3}\quad \ \text{ for all }\ \psi \in H^3\,.$$

Proof.

The bound follows by applying Lemma 3.2 to the terms in $[\widehat{T},\widehat{V}](\psi )$. We note that the first three terms can be estimated using only the $H^2$ norm, but the last term requires a stronger norm.

■

The estimate of the local error is now obtained with a nonlinear version of the analysis of splitting methods by Jahnke and Lubich 14, similar to Lubich 17; cf. also Kozlov, Kværnø and Owren 16 for another related technique.

4.3. Preparation: Lie derivatives

We use the calculus of Lie derivatives (see, e.g., 9, Sect. III.5 or 12, Sect. IV.1.4). Since this formalism only relies on the differentiability and the semi-group property of the flow, it is applicable in the present infinite-dimensional setting as well as in the finite-dimensional case. For a vector field $F$ on $H^1$, such as $\widehat{T}$ or $\widehat{V}$ or $\widehat{H}=\widehat{T}+\widehat{V}$, we denote by $\varphi _F^t$ the flow at time $t$ of the differential equation $\dot{\psi }=F(\psi )$, that is, $\varphi _F^t(v)$ is the solution at time $t$ of this differential equation with initial value $\psi (0)=v$. We consider the Lie derivative $D_F$ defined by

$$\bigl (D_F G \bigr )(v) = {d\over dt}\Big |_{t=0} G(\varphi _F^t(v)) = G'(v) F(v)$$

for another vector field $G$ on $H^1$ and $v\in H^1$, and we set

$$\Bigl (\exp (tD_F) G\Bigr )(v) = G(\varphi _F^t(v)).$$

In particular, for the identity ${\mathrm{Id}}$, the flow is reproduced as $\exp (tD_F) {\mathrm{Id}}(v) = \varphi _F^t(v)$. We then have the rule

$${d \over dt} \exp (t D_F) G(v) = \Bigl (D_F \exp (t D_F) G\Bigr )(v) = \Bigl (\exp (t D_F) D_F G\Bigr )(v).$$

The commutator $[D_F , D_G]=D_F D_G - D_G D_F$ of the Lie derivatives of two vector fields $F$ and $G$ is the Lie derivative of the commutator of the vector fields in reversed order:

$$[D_F , D_G] = D_{[G,F]}.$$

4.4. Local error: Proof of Proposition 2.3

(a) For notational simplicity we write $D_H$, $D_T$, $D_V$ instead of $D_{\widehat{H}}$, $D_{\widehat{T}}$, $D_{\widehat{V}}$, respectively. We start from the nonlinear variation-of-constants formula

$$\begin{align*} \psi (\tau ) &= \exp (\tau D_{H}) {\mathrm{Id}}\, (\psi _0) = \exp (\tau D_{T}) {\mathrm{Id}} \,(\psi _0) \\ &\qquad +\int _0^\tau \exp ((\tau -s)D_{H}) \, D_{V} \, \exp (s D_T) {\mathrm{Id}}\, (\psi _0) \, ds. \end{align*}$$

Using this formula once more for the expression under the integral, we obtain

$$\begin{align*} \psi (\tau ) &= \exp (\tau D_{T}) {\mathrm{Id}} \,(\psi _0) \\ &\qquad +\int _0^\tau \exp ((\tau -s)D_{T}) \, D_{V} \, \exp (s D_T) {\mathrm{Id}} \,(\psi _0) \, ds + r_1 \end{align*}$$

with the remainder

$$r_1 = \int _0^\tau \int _0^{\tau -s} \exp ((\tau -s-\sigma )D_H) D_V \exp (\sigma D_T) D_V \exp (s D_T) {\mathrm{Id}} \,(\psi _0) \, d\sigma \, ds\nobreakspace .$$

On the other hand, in this notation the numerical solution reads

$$\psi _1 = \exp ({\textstyle {1\over 2}}\tau D_T) \exp (\tau D_V) \exp ({\textstyle {1\over 2}}\tau D_T) {\mathrm{Id}} \,(\psi _0),$$

and Taylor expansion $\exp (\tau D_V)= I + \tau D_V + \tau ^2 \int _0^1 (1-\theta )\exp (\theta \tau D_V) D_V^2\, d\theta$ gives

$$\psi _1 = \exp (\tau D_T) {\mathrm{Id}} \,(\psi _0) + \tau \exp ({\textstyle {1\over 2}}\tau D_T) D_V \exp ({\textstyle {1\over 2}}\tau D_T) {\mathrm{Id}} \,(\psi _0) + r_2$$

with the remainder

$$r_2 = \tau ^2 \int _0^1 (1-\theta )\exp ({\textstyle {1\over 2}}\tau D_T)\exp (\theta \tau D_V) D_V^2\exp ({\textstyle {1\over 2}}\tau D_T) {\mathrm{Id}}\,(\psi _0)\, d\theta \,.$$

(b) The error now becomes

$$\begin{eqnarray} \psi _1-\psi (\tau ) &=& \tau \, \exp ({\textstyle {1\over 2}}\tau D_T) D_V \exp ({\textstyle {1\over 2}}\tau D_T) {\mathrm{Id}} \,(\psi _0) \\ &&- \int _0^\tau \exp ((\tau -s)D_{T}) \, D_{V}\, \exp (s D_T) {\mathrm{Id}}\, (\psi _0) \, ds + (r_2 -r_1), {\tag{4.2}} \end{eqnarray}$$

and hence the principal error term is just the quadrature error of the midpoint rule applied to the integral over $[0,\tau ]$ of the function

$$\begin{equation} f(s)=\exp ((\tau -s)D_{T}) \, D_{V}\, \exp (s D_T) {\mathrm{Id}} \,(\psi _0). {\tag{4.3}} \end{equation}$$

We express the quadrature error in first-order Peano form,

$$\tau \, f({\textstyle {1\over 2}}\tau ) - \int _0^\tau f(s)\, ds = \tau ^2 \int _0^1 \kappa _1(\theta )\, f'(\theta \tau )\, d\theta$$

with the (scalar, bounded) Peano kernel $\kappa _1$ of the midpoint rule. Since

$$\begin{eqnarray*} f'(s) &=& -\exp ((\tau -s)D_{T}) \, [D_T,D_{V}]\, \exp (s D_T) {\mathrm{Id}}\, (\psi _0) \\ &=& \exp ((\tau -s)D_{T}) \, D_{[\widehat{T}, \widehat{V}]}\, \exp (s D_T) {\mathrm{Id}} \,(\psi _0) \\ &=& e^{is\Delta } \, [\widehat{T}, \widehat{V}] ( e^{i(\tau -s)\Delta } \psi _0)\nobreakspace , \end{eqnarray*}$$

the commutator bound of Lemma 4.1 shows that the quadrature error is bounded by

$$\begin{equation} \Bigl \| \tau \, f({\textstyle {1\over 2}}\tau ) - \int _0^\tau f(s)\, ds \Bigr \|_{H^1} \le C\tau ^2 \| \psi _0 \|_{H^2}^3 \,. {\tag{4.4}} \end{equation}$$

(c) Finally, we estimate the remainder terms. For $\|\psi _0\|_{H^1}\le M_1$, we now show that they are bounded by

$$\begin{equation} \| r_1 \|_{H^1} + \| r_2 \|_{H^1} \le C_1 \tau ^2\,, {\tag{4.5}} \end{equation}$$

where $C_1$ only depends on $M_1$. From the definitions we have

$$\begin{eqnarray*} \exp (\rho D_H) D_V \exp (\sigma D_T) D_V \exp (s D_T) {\mathrm{Id}} \,(\psi _0) &=& e^{is\Delta }\widehat{V}'(e^{i\sigma \Delta }\psi (\rho )) e^{i\sigma \Delta } \widehat{V}(\psi (\rho )), \\ \exp ({\textstyle {1\over 2}}\tau D_T)\exp (\theta \tau D_V) D_V^2\exp ({\textstyle {1\over 2}}\tau D_T) {\mathrm{Id}}\,(\psi _0) &=& e^{i\tau \Delta /2} \widehat{V}'(\eta ) \widehat{V}(\eta ), \end{eqnarray*}$$

where $\eta =e^{-i\theta \tau V[\phi ]}\,\phi$ with $\phi =e^{i\tau \Delta /2}\psi _0$, with $\|\eta \|_{H^1} \le e^{a_1\tau }\, \|\psi _0\|_{H^1}$ by (4.1). Since Lemma 3.2 yields the bounds (for $\psi$ of unit $L_2$ norm)

$$\begin{equation} \| \widehat{V}(\psi ) \|_{H_1} \le C \,\|\psi \|_{H^1}^2\quad \text{ and }\quad \| \widehat{V}'(\psi )\phi \|_{H_1} \le C \,\|\psi \|_{H^1}^2\, \|\phi \|_{H^1}\,, {\tag{4.6}} \end{equation}$$

we obtain the bound (4.5).■

4.5. Proof of the $H^1$ error bound of Theorem 2.1

The stated error bound follows from Propositions 2.2 and 2.3 with the standard argument of Lady Windermere’s fan 10, Sect. II.3. Note that the boundedness in $H^1$ required by the stability lemma, is ensured by induction by the $H^1$ error bound.■

5. Proof of the second-order error bound in $L_2$

5.1. Double-commutator bound

Lemma 5.1.

The double commutator of $\widehat{T}$ with $\widehat{V}$ is bounded in $L_2$ by

$$\| [\widehat{T}, [\widehat{T}, \widehat{V}]] (\psi ) \|_{L_2} \le C \,\| \psi \|_{H^4}^3 \qquad \text{for all }\ \psi \in H^4.$$

Proof.

Direct calculation shows that among a plethora of more harmless terms that can be bounded, by Lemmas 3.1 and 3.2, in terms of the $H^3$ or even $H^2$ norms, the double commutator contains also the term $4i\Delta ^{-1} (\psi {\Delta ^2\overline{\psi }})\psi$, which can be bounded in terms of the $H^4$ norm.

■

5.2. Local error in $L_2$: Proof of Proposition 2.4

(a) We return to the error formula (4.2) and write the principal error term in second-order Peano form

$$\tau \, f({\textstyle {1\over 2}}\tau ) - \int _0^\tau f(s)\, ds = \tau ^3 \int _0^1 \kappa _2(\theta )\, f''(\theta \tau )\, d\theta$$

with the Peano kernel $\kappa _2$ of the midpoint rule and $f$ of (4.3). We have

$$\begin{eqnarray*} f''(s) &=& \exp ((\tau -s)D_{T}) \, [D_T,[D_T,D_{V}]]\, \exp (s D_T) {\mathrm{Id}}\, (\psi _0) \\ &=& \exp ((\tau -s)D_{T}) \, D_{[\widehat{T}, [\widehat{T}, \widehat{V}]]}\, \exp (s D_T) {\mathrm{Id}} \,(\psi _0) \\ &=& e^{is\Delta } \, [\widehat{T}, [\widehat{T}, \widehat{V}]] ( e^{i(\tau -s)\Delta } \psi _0)\nobreakspace , \end{eqnarray*}$$

and hence Lemma 5.1 shows that the quadrature error is bounded in $L_2$ by $C\tau ^3 \| \psi _0 \|_{H^4}^3$.

(b) With the function

$$g(s,\sigma ) = \exp ((\tau -s-\sigma )D_T) D_V \exp (\sigma D_T) D_V \exp (s D_T) {\mathrm{Id}} \,(\psi _0)\,$$

the remainder term can be expressed as

$$r_2 -r_1 = \frac{\tau ^2}{2} \, g\Bigl ( \frac{\tau }{2}, 0\Bigr ) - \int _0^\tau \int _0^{\tau -s} g(s,\sigma ) \, d\sigma \, ds +\widetilde{r}_2 - \widetilde{r}_1,$$

where, in the same way as in part (c) of Section 4.4, the remainders can be bounded by

$$\| \widetilde{r}_1 \|_{L_2} + \| \widetilde{r}_2 \|_{L_2} \le \widetilde{C}_2 \tau ^3$$

with $\widetilde{C}_2$ depending only on $\|\psi _0\|_{H^2}$. The other two terms in $r_2-r_1$ form the quadrature error of a first-order two-dimensional quadrature formula, and are therefore bounded by

$$\Bigl \| \frac{\tau ^2}{2} \, g\Bigl ( \frac{\tau }{2}, 0\Bigr ) - \int _0^\tau \int _0^{\tau -s} g(s,\sigma ) \, d\sigma \, ds \Bigr \|_{L_2} \le C \tau ^3 \Bigl ( \max \Big \| \frac{\partial g}{\partial s} \Big \|_{L_2} + \max \Big \| \frac{\partial g}{\partial \sigma } \Big \|_{L_2} \Big )\,,$$

where the maxima are taken over the triangle $0\le s \le \tau$, $0\le \sigma \le \tau -s$. The partial derivatives of $g$,

$$\begin{eqnarray*} \frac{\partial g}{\partial s} (s,\sigma ) &=& \exp ((\tau -s-\sigma )D_T) D_{[\widehat{T}, \widehat{V}]} \exp (\sigma D_T) D_V \exp (s D_T) {\mathrm{Id}} \,(\psi _0) \\ && +\ \exp ((\tau -s-\sigma )D_T) D_V \exp (\sigma D_T) D_{[\widehat{T}, \widehat{V}]} \exp (s D_T) {\mathrm{Id}} \,(\psi _0), \\ \frac{\partial g}{\partial \sigma } (s,\sigma ) &=& \exp ((\tau -s-\sigma )D_T) D_{[\widehat{T}, \widehat{V}]} \exp (\sigma D_T) D_V \exp (s D_T) {\mathrm{Id}} \,(\psi _0)\,, \end{eqnarray*}$$

only contain $\widehat{V}$ and the simple commutator $[\widehat{T}, \widehat{V}]$ and their derivatives. The $L_2$ norms of ${\partial g}/{\partial s}$ and ${\partial g}/{\partial \sigma }$ can therefore be bounded in terms of the $H^2$ norm of $\psi _0$ using (4.6) and the argument of the proof of Lemma 4.1. Together, this shows

$$\| r_2 -r_1 \|_{L_2} \le C_2 \tau ^3\,,$$

where $C_2$ only depends on $\|\psi _0\|_{H^2}$. Recalling the error formula (4.2) and combining the above bound with that of part (a) yields the result of Proposition 2.4.■

5.3. Proof of the $L_2$ error bound of Theorem 2.1

With the $H^2$ regularity of the exact solution, with the $L_2$ bound of the local error of Proposition 2.4, and with the $H^1$-conditional $L_2$-stability of Proposition 2.2 together with the $H^1$ bound of the numerical solution established in Section 4, the result is obtained with the standard argument of Lady Windermere’s fan 10, Sect. II.3.■

6. $H^2$ regularity: Proof of Proposition 2.5

Since $e^{i\tau \Delta }$ preserves the $H^2$ norm, we only need to bound the $H^2$ norm of $e^{-i\tau V[\phi ]}\phi$ for $\phi \in H^2$, which is the solution at time $\tau$ of

$$i\dot{\eta }= V[\phi ]\eta , \quad \ \eta (0)=\phi .$$

By Lemma 3.2 and $\|\eta \|_{L_2}=\|\phi \|_{L_2}=1$,

$$\| V[\phi ]\eta \|_{H^2} \le K_2 \left( \|\phi \|_{H^1}\, \|\eta \|_{H^2} + \|\phi \|_{H^2} \, \|\eta \|_{H^1} + \|\phi \|_{H^2} \, \|\phi \|_{H^1} \right).$$

By (4.1) we have

$$\|\eta (t)\|_{H^1} \le e^{a_1 t} \, \|\phi \|_{H^1}\,,$$

where $a_1$ depends only on $M_1$. For the $H^2$ norm we then obtain

$$\|\eta (t)\|_{H^2} \le \|\phi \|_{H^2} + \int _0^t C_1 \left( \|\phi \|_{H^2} + \|\eta (s)\|_{H^2} \right)\,ds$$

where $C_1$ only depends on $M_1$, and hence once again by the Gronwall inequality,

$$\|\eta (t)\|_{H^2} \le e^{a_2 t} \|\phi \|_{H^2}$$

where again $a_2$ only depends on $M_1$. Combining these estimates yields Proposition 2.5.■

PART B. THE CUBIC NONLINEAR SCHRÖDINGER EQUATION

7. Error bounds for solutions in $H^4$: Statement of results

For the cubic nonlinear Schrödinger equation (1.1), (1.2) with solutions in $H^4$ similar results are obtained. We suppose that the solution $\psi (t)$ to the cubic Schrödinger equation (1.1), (1.2) is in $H^4$ for $0\le t \le T$, and set

$$m_4 = \max _{0\le t \le T} \| \psi (t) \|_{H^4} \,.$$

Theorem 7.1.

The numerical solution $\psi _n$ given by the splitting scheme 1.4 for the cubic nonlinear Schrödinger equation with step size $\tau >0$ has a first-order error bound in $H^2$ and a second-order error bound in $L_2$,

$$\begin{array}{l} \| \psi _n - \psi (t_n) \|_{H^2} \le C(m_4,T)\, \tau \\[2.84526pt] \| \psi _n - \psi (t_n) \|_{L_2} \,\le C(m_4,T)\, \tau ^2 \end{array} \qquad \text{ for } t_n=n\tau \le T\,.$$

We again write the step of the splitting scheme (1.4) briefly as

$$\psi _{n+1} = \Phi _\tau (\psi _n)\,.$$

Proposition 7.2 ($H^2$-conditional $L_2$-, $H^1$- and $H^2$-stability).

If $\psi ,\phi \in H^2$ with

$$\|\psi \|_{H^2} \le M_2,\quad \ \|\phi \|_{H^2} \le M_2,$$

then

$$\begin{eqnarray*} &&\| \Phi _\tau (\psi )-\Phi _\tau (\phi ) \|_{L_2} \le e^{c_0\tau } \, \|\psi -\phi \|_{L_2}, \\ &&\| \Phi _\tau (\psi )-\Phi _\tau (\phi ) \|_{H^1} \le e^{c_1\tau } \, \|\psi -\phi \|_{H^1}, \\ &&\| \Phi _\tau (\psi )-\Phi _\tau (\phi ) \|_{H^2} \le e^{c_2\tau } \, \|\psi -\phi \|_{H^2}, \end{eqnarray*}$$

where $c_0,c_1, c_2$ only depend on $M_2$.

Note that the $L_2$- and $H^1$-stability estimates depend on bounds in $H^2$.

Proposition 7.3 (Local error in $H^2$).

If $\psi _0\in H^4$ with $\|\psi _0\|_{H^4} \le M_4$, then the error after one step of the method 1.4 is bounded in the $H^2$ norm by

$$\| \psi _1 - \psi (\tau ) \|_{H^2} \le C_4\tau ^2\,,$$

where $C_4$ only depends on $M_4$.

Proposition 7.4 (Local error in $L_2$).

If $\psi _0\in H^4$ with $\|\psi _0\|_{H^4} \le M_4$, then the error after one step of the method 1.4 is bounded in the $L_2$ norm by

$$\| \psi _1 - \psi (\tau ) \|_{L_2} \le C_4\tau ^3\,,$$

where $C_4$ only depends on $M_4$.

There is also an analogue of Proposition 2.5 for the cubic nonlinear Schrödinger equation, inferring $H^3$-regularity of the numerical solution from bounds in $H^2$.

For the one-dimensional cubic nonlinear Schrödinger equation we would obtain also $H^1$-conditional stability (essentially because $H^1({\mathbf{R}})\subset L_\infty ({\mathbf{R}})$).

8. Outline of the proofs

The proof of Theorem 7.1 and the above propositions is analogous to the corresponding results for the Schrödinger-Poisson equation. Essentially, the operator $\Delta ^{-1}$ is to be replaced by the identity operator in all formulas. The estimates of Lemma 3.1 need to be replaced by

$$\begin{eqnarray*} \| uvw \|_{L_2} &\le & K_0 \,\| u \|_{H^1}\,\| v \|_{H^1}\, \| w \|_{H^1}, \\ \| uvw \|_{L_2} &\le & K_0 \, \| u \|_{L_2}\,\| v \|_{H^2}\, \| w \|_{H^2}\,. \end{eqnarray*}$$

The first bound follows from the Sobolev embedding $H^1\subset L_6$, and the second bound from the Sobolev embedding $H^2\subset L_\infty$. We then have the further bounds

$$\begin{eqnarray*} \| uvw \|_{H^1} &\le & K_1 \,\| u \|_{H^1}\,\| v \|_{H^2}\, \| w \|_{H^2}, \\ \| uvw \|_{H^2} &\le & K_2 \, \| u \|_{H^2}\,\| v \|_{H^2}\, \| w \|_{H^2}\,. \end{eqnarray*}$$

The commutator bounds now become

$$\begin{array}{rcr} \| [\widehat{T}, \widehat{V}] (\psi ) \|_{H^2} &\le & C \| \psi \|_{H^4}^3, \phantom {\,.} \\[2.84526pt] \| [\widehat{T}, [\widehat{T}, \widehat{V}]] (\psi ) \|_{L_2} &\le & C \| \psi \|_{H^4}^3 \,. \end{array}$$

With these bounds the results follow in the same way as before.