1 Introduction

In this paper we consider nonlinear Schrödinger equations (NLS) seeking a complex function u(t, x) such that

$$\begin{aligned} \mathrm {i}\partial _t u = - \triangle u + V u + \gamma ( |u|^2 ) u \end{aligned}$$

in a bounded domain \(\mathscr {D} \subset \mathbb {R}^d\), with a homogenous Dirichlet boundary condition on \(\partial \mathscr {D}\) and a given initial value. Here, V(x) is a known real-valued potential and \(\gamma : [0,\infty ) \rightarrow \mathbb {R}\) is a smooth (and possibly nonlinear) function that depends on the unknown density \(|u|^2\). Of particular interest are cubic nonlinearities of the form \(\gamma (|u|^2)u=\kappa |u|^2u\), for some \(\kappa \in \mathbb {R}\). In this case, the equation is called Gross–Pitaevskii equation. It has applications in optics [1, 2], fluid dynamics [3, 4] and, most importantly, in quantum physics, where it models for example the dynamics of Bose-Einstein condensates in a magnetic trapping potential [5,6,7]. Another relevant class are saturated nonlinearities, such as \(\gamma (|u|^2)=\kappa |u|^2 (1 + \alpha |u|^2)^{-1}\) for some \(\alpha \ge 0\), which appear in the context of nonlinear optical wave propagation in layered metallic structures [8, 9] or the propagation of light beams in plasmas [10]. In order to discretize nonlinear Schrödinger equations in time, splitting methods and exponential integrators yield typically highly efficient solution schemes that can be easily combined with a spectral discretization in space (cf. [11,12,13,14,15,16,17,18] and the references therein). If the exact solution to the NLS admits high regularity, such discretization schemes typically show a remarkably good performance. However, if the regularity of the solution is strongly reduced, either by rough potentials V (e.g. disorder potentials or optical lattices) or by rough initial values u(0) (e.g. when effects close to phase transitions are studied), then the performance of these methods can drop dramatically. Here we refer exemplarily to the recent numerical experiments reported in [19,20,21]. To overcome this issue, Ostermann and Schratz proposed new low-regularity time-integrators [19, 20] which improve the convergence in low regularity regimes significantly. However, the approach still relies on a Fourier discretization in space, which is not an optimal choice due to the loss of spectral convergence for non-smooth solutions. Practically, the usage of a (low order) finite element space discretization is often desirable in order to account for spatial low regularity. In the following we will only discuss approaches that can be easily combined with finite elements in space, meaning that we put ourselves into the situation that we assume that the solution to the NLS does not admit much smoothness.

Nonlinear Schrödinger equations come with important physical invariants, where the mass and the energy are considered as two of the most crucial ones. When solving a NLS numerically it is therefore of great importance to also reproduce this conservation on the discrete level. This aspect was emphasized by various numerical studies [21, 22], where it was also found that the complexity of the physical setup (or low-regularity) can stress this issue even further .

For the subclass of power law nonlinearities of the form \(\gamma (|u|^2)= \sum _{k=1}^K \kappa _k |u|^{2\sigma _k}\) for \(\sigma _k \ge 0\) and \(\alpha _k \in \mathbb {R}\), a mass and energy conserving relaxation scheme was proposed and analyzed by Besse [23, 24]. Thanks to its properties, the scheme shows a very good performance in realistic physical setups [21]. Despite the large variety of different numerical approaches for solving the time-dependent NLS (cf. [16,17,18, 25,26,27,28,29,30,31,32,33,34] and the references therein) the literature knows however only one time discretization that conserves both mass and energy simultaneously for arbitrary (smooth) nonlinearities. This discretization, which was first mathematically studied by Sanz-Serna [31] and which is long-known in the physics community, is a Crank-Nicolson-type (CN) approach where the nonlinearity is approximated by a suitable difference quotient involving the primitive integral of \(\gamma \). This is also the time discretization that we shall consider in this paper. Here we note it was analytically and numerically demonstrated that this method is applicable and reliable in low-regularity regimes [21, 35].

A combination of the method with a finite difference space discretization was proposed and analyzed by Bao and Cai [36, 37]. Combining the Crank-Nicolson time discretization with a P1 finite element discretization in space, the first a priori error estimates for the arising method were obtained by Sanz-Serna 1984 [31] for cubic nonlinearities. He considers the case \(d=1\) and derives optimal \(L^{\infty }(L^2)\)-error estimates under the coupling constraint \(\tau _{}\lesssim h\), where \(\tau _{}\) denotes the time step size of the Crank-Nicolson method and h the mesh size of the finite element discretization in space. In 1991, Akrivis et al. [25] improved this result by showing optimal convergence rates in \(L^{\infty }(L^2)\) in dimension \(d=1,2,3\) and under the relaxed coupling constraint \(\tau _{} \lesssim h^{d/4}\). Finally, in 2017 [35], the \(L^{\infty }(L^2)\)-error estimates could be improved yet another time by showing that the coupling constraint can be fully removed. Furthermore, general nonlinearities could be considered, the influence of potentials could be taken into account and even convergence under weak regularity assumptions could be proved (with reduced convergence rates). However, so far, optimal error estimates in \(L^{\infty }(H^1)\) for this particular CN-discretization are still open in the literature.

One reason for this absence of \(H^1\)-results could be related to the techniques used for the error analysis in previous works (cf. [25, 28,29,30,31, 34]) which is based on the following steps: 1. Appropriate truncation of the nonlinearity to obtain a problem with bounded growth. 2. Analyzing the scheme with truncation in the FE space and deriving corresponding \(L^{\infty }(L^2)\)- and/or \(L^{\infty }(H^1)\)-error estimates. 3. Using inverse estimates in the finite element space to show that the truncated approximations are uniformly bounded in \(L^{\infty }(L^{\infty })\) by a term of the form \(C (1 + h^{-s} ( \tau _{}^2 + h^{p}))\), with appropriate powers \(p>s>0\) that depend on the considered space discretization, regularity and space dimension. 4. Concluding that if \(\tau _{}\) and h are coupled in an appropriate way, then the truncated approximations are all uniformly bounded by a constant C and hence coincide with a solution to the scheme without truncation.

This strategy does not only have the disadvantage that it produces unnecessary coupling conditions, but also that it becomes impractically technical when considering \(L^{\infty }(H^1)\)-error estimates for the Crank-Nicolson FEM. This is because it requires a suitable truncation of the primitive integral of \(\gamma \) that is on the one hand consistent with the energy conservation and on the other hand allows for uniform bounds of the approximations in \(L^{\infty }(W^{1,\infty })\). However, thanks to the new techniques developed in [33] and the CN error analysis suggested in [35] in the context of \(L^{\infty }(L^2)\)-error estimates, the truncation step is no longer necessary and the desired \(L^{\infty }(L^{\infty })\)-bounds can be derived with elliptic regularity theory. With this, it is now possible to obtain \(L^{\infty }(H^1)\) estimates in a direct way, not only in the finite element setting, but also in the semi-discrete Hilbert space setting.

In this paper we will therefore build upon the results from [33, 35] to fill the gap in the literature and prove optimal \(L^{\infty }(H^1)\)-error estimates for the energy-conservative Crank-Nicolson approach without coupling constraints and for a general class of nonlinearities. The paper is structured as follows. In Sect. 2 we present the notation and the analytical assumptions on the problem. In Sect. 3 we present the time–discrete Crank-Nicolson method, we recall its well-posedness and optimal error estimates in \(L^{\infty }(L^2)\). Furthermore, we present and prove the new error estimate in \(L^{\infty }(H^1)\). The paper continues with the fully-discrete setting presented in Sect. 4, where the time discretization is combined with a finite element discretization in space. We recall what is known about this discretization and finally prove corresponding \(L^{\infty }(H^1)\)-error estimates, which is the main result of this paper. The paper concludes with a note on how to efficiently implement the method and two numerical experiments to confirm the convergence rates and to illustrate a setting in which it makes computational sense to use the CN-FEM instead of for example a spectral method.

2 Notation and assumptions

We start with introducing the analytical setting of this work. Throughout the paper we assume that \({\mathscr {D}} \subset \mathbb {R}^d\) (for \(d=2,3\)) is a convex bounded domain with polyhedral boundary. On \(\mathscr {D}\), the Sobolev space of complex-valued, weakly differentiable functions with a zero trace on \(\partial \mathscr {D}\) and \(L^2\)-integratable partial derivatives is as usual denoted by \(H^1_0(\mathscr {D}):=H^1_0(\mathscr {D},\mathbb {C})\). The potential \(V \in L^{\infty }(\mathscr {D};\mathbb {R})\) is assumed to be real and nonnegative. Indirectly, we also assume that V is sufficiently smooth so that it is compatible with the regularity assumptions for u listed below (see [35] for a discussion on this aspect). The (possibly nonlinear) function

$$\begin{aligned}\gamma : [0,\infty ) \rightarrow [0,\infty )\end{aligned}$$

is assumed to be two times differentiable, fulfills \(\gamma (0)=0\) and its growth can be characterized with

$$\begin{aligned} |\gamma (|v|^2)v - \gamma (|w|^2)w| \le L(K) |v - w| \qquad \hbox {for all } v,w \in \mathbb {C}\hbox { with } |v|,|w|\le K \end{aligned}$$

where L is a function with the following growth properties

$$\begin{aligned} 0 \le L(s) \le C s^q \quad \hbox {for } s \ge 0 \qquad \hbox {and } {\left\{ \begin{array}{ll} q \in [0,\infty ) &{}\quad \hbox {for } d=2,\\ q \in [0,4) &{}\quad \hbox {for } d=3. \end{array}\right. } \end{aligned}$$

Note that in [35] the admissible growth condition in 3d requires \(q \in [0,2)\), which is however a typo and should be, as above, \(q \in [0,4)\) (cf. [38, Proposition 3.2.5 and Remark 3.2.7] for the original result). Examples for nonlinearities that fulfill these assumptions are mentioned in the introduction. The most common and physically relevant choices covered by our setting are power law nonlinearities \(\gamma (\rho )=\kappa \rho ^{q}\) for \(\kappa \ge 0\) and \(0\le q<\infty \) in 2d and \(0\le q < 4\) in 3d. Other physically relevant nonlinearities that fulfill the conditions are saturated nonlinearities appearing in the modeling of optical wave propagation such as \(\gamma (\rho )=\kappa \rho (1+ \alpha \rho )^{-1}\) for \(\alpha ,\kappa \ge 0\).

The above assumptions cover the regime of so-called defocussing (positive) nonlinearities and guarantees that the NLS and its Crank-Nicolson discretization are well-posed. For focussing (negative) nonlinearities, i.e. \(\gamma : [0,\infty ) \mapsto (-\infty ,0]\), the well-posedness (of both the continuous and discrete models) can no longer be guaranteed without making additional technical assumptions. Typically, effects such as finite time blow ups can occur in this regime. To avoid constantly having to invoke a saving clause we restrict our attention to the defocussing case. We do however point out that under the assumptions that the NLS and the Crank-Nicolson discretizations are well-posed (without blow-up in the time interval [0, T]) then all our error estimates hold without changes.

For the initial value we assume that \(u^0 \in H^1_0(\mathscr {D}) \cap H^2(\mathscr {D})\) and, without loss of generality, that it has a normalized mass, i.e. \(\int _{\mathscr {D}} |u^0(x)|^2 \,dx =1\). With this, the considered nonlinear Schrödinger equation (NLS) reads as follows. For a maximum time \(T>0\) and an initial value \(u^0\), we seek

$$\begin{aligned} u \in L^{\infty }([0,T],H^1_0({\mathscr {D}})) \quad \hbox {and} \quad \partial _t u \in L^{\infty }([0,T],H^{-1}({\mathscr {D}})) \end{aligned}$$

such that \(u(\cdot ,0)=u^0 \) and

$$\begin{aligned} \mathrm {i}\partial _t u = - \varDelta u + V \, u + \gamma ( |u|^2 )\, u \end{aligned}$$
(2.1)

in the sense of distributions. Problem (2.1) admits at least one solution, that is even unique for repulsive cubic nonlinearities in 1d and 2d (cf. [38] in general and [35, Remark 2.1] for precise references). We assume that the solution admits the following additional regularity, which is

$$\begin{aligned}&u_{ttt} \in {L^2(0,T;H^1(\mathscr {D}))} \nonumber \\ \hbox {and} \quad&u \in {L^{\infty }(0,T; W^{1,\infty } ({\mathscr {D}}))}, \quad \partial _t^{(k)} u \in L^2(0,T;H^2(\mathscr {D})) \quad \hbox {for } 0\le k \le 2, \end{aligned}$$
(2.2)

where we note that any solution with such increased regularity must be unique (cf. [35, Lemma 3.1]). In the rest of the paper u hence always refers to this uniquely characterized solution.

It is well known that solutions to the NLS (2.1) preserve the mass, i.e.

$$\begin{aligned} \int _{\mathscr {D}} |u(t,x)|^2 \,dx = \int _{\mathscr {D}} |u^0(x)|^2 \, dx = 1 \end{aligned}$$

and the energy, i.e.

$$\begin{aligned} E[ u(t) ] = E[u^0], \quad \hbox {where } E[u]:= \frac{1}{2} \int _{\mathscr {D}} |\nabla u(x)|^2 + V(x) \, |u(x)|^2 + \varGamma (|u(x)|^2) \, dx, \end{aligned}$$

with \(\varGamma (\rho ) := \int _0^{\rho } \gamma (r)\,dr\).

For brevity, we shall denote the \(L^2\)-norm of a function \(v \in L^2(\mathscr {D}):=L^2(\mathscr {D},\mathbb {C})\) by \(\Vert v \Vert \). The \(L^2\)-inner product is denoted by \(\langle v, w \rangle =\int _{\mathscr {D}} v(x) \, \overline{w(x)} \, dx\). Here, \(\overline{w}\) denotes the complex conjugate of w.

Throughout the paper we will use the notation \(A \lesssim B\), to abbreviate \(A \le C B\), where C is a constant that only depends on u, T, d, \(\mathscr {D}\), V and \(\gamma \), but not on the discretization.

Remark 2.1

In the analysis we restrict our attention to homogeneous Dirichlet boundary conditions. Typically these boundary conditions can be motivated by physical reasoning. For example in the context of Bose Einstein condensates, the magnetic potential V is a trapping potential that becomes very quickly very large and hence traps the condensate in a bounded region. Mathematically this leads to an exponential decay of the solution u to zero (in moderate distances from the origin of the coordinate system) and hence justifies to truncate the computational domain to a simple geometric object on which the problem is solved with zero boundary conditions. A typical alternative found in the literature are periodic boundary conditions which are e.g. favorable for spectral methods. Both the formulation of the Crank-Nicolson method and its error analysis can be easily generalized to that case.

3 Time-discrete Crank-Nicolson scheme

In this section we will state the semi-discrete Crank-Nicolson scheme, recall its well-posedness and available stability bounds, and then use these results to prove optimal \(L^{\infty }(H^1)\)-error estimates in the Hilbert space setting. For that, let T denote the final time of computation, N the number of time-steps, and \(\tau = T/N\) the time step size. By \(t_n\) we shall mean \(t_n=n\tau \). The exact solution at time \(t_n\) shall be denoted by \(u^n:=u(t_n,\cdot )\). We also introduce a short hand notation for discrete time derivatives which is \(D_\tau u^n := (u^{n+1}-u^n)/\tau \) and analogously \(D_\tau u_{\tau _{}}^{n} := ( u_{\tau _{}}^{n+1} - u_{\tau _{}}^{n})/\tau \).

3.1 Method formulation and main result

With the notation above, the semi-discrete Crank-Nicolson approximation \(u_{\tau _{}}^{n+1} \in H^1_0(\mathscr {D})\) to \(u^{n+1}\) is given recursively as the solution (in the sense of distributions) to the equation

$$\begin{aligned} \mathrm {i}D_\tau u_{\tau _{}}^{n} = - \varDelta u_{\tau _{}}^{n+\frac{1}{2}} + V u_{\tau _{}}^{n+\frac{1}{2}} + \frac{\varGamma (|u_{\tau _{}}^{n+1}|^2)-\varGamma (|u_{\tau _{}}^{n}|^2)}{|u_{\tau _{}}^{n+1}|^2-|u_{\tau _{}}^{n}|^2} \, u_{\tau _{}}^{n+\frac{1}{2}}, \end{aligned}$$
(3.1)

where \(u_{\tau _{}}^{n+\frac{1}{2}}:=(u_{\tau _{}}^{n}+u_{\tau _{}}^{n+1})/2\). The initial value is selected as \(u_{\tau _{}}^{0}=u^0\). It is easily seen that the discretization conserves both mass and energy, i.e.

$$\begin{aligned} \int _{\mathscr {D}} |u_{\tau _{}}^{n}|^2 \, dx = \int _{\mathscr {D}} |u^0|^2 \, dx \quad \hbox {and} \quad E[u_{\tau _{}}^{n}] = E[u^0] \quad \hbox {for all } n\ge 0. \end{aligned}$$

The scheme (3.1) is well-posed and admits a set of a priori error estimates. The properties are summarized in the following theorem that is proved in [35, Theorem 4.1].

Theorem 3.1

Under the general assumptions of this paper, there exists a constant \(C(u)>0\) and a solution \(u_{\tau _{}}^{n}\in H^1_0(\mathscr {D})\) to the semi-discrete Crank-Nicolson scheme (3.1) that is uniquely characterized by the property that

$$\begin{aligned} \sup _{0\le n \le N} \left( \Vert u_{\tau _{}}^{n} \Vert _{L^{\infty }(\mathscr {D})} + \Vert u_{\tau _{}}^{n} \Vert _{H^{2}(\mathscr {D})}\right) \le C(u), \end{aligned}$$
(3.2)

and the a priori estimate for the \(L^2\)-error

$$\begin{aligned} \sup _{0\le n \le N} \Vert u_{\tau _{}}^{n} - u^n \Vert \lesssim \tau _{}^2, \end{aligned}$$

where u is the (unique) exact solution with the regularity property (2.2).

Our main result on optimal error estimates in the \(L^{\infty }(H^1)\) reads as follows.

Theorem 3.2

(Optimal \(H^1\)-error estimates for the semi-discrete method) Consider the setting of Theorem 3.1, then the \(L^{\infty }(H^1)\)-error converges with optimal order in \(\tau \), i.e.

$$\begin{aligned} \sup _{0\le n \le N} \Vert u_{\tau _{}}^{n} - u^n \Vert _{H^{1}(\mathscr {D})} \lesssim \tau _{}^2. \end{aligned}$$

The theorem is proved in Sect. 3.2 below.

3.2 Proof of Theorem 3.2

In this section we will prove Theorem 3.2. Let us introduce some notation that is used throughout the proofs. We recall \(D_\tau e^n = ( e^{n+1} - e^{n})/\tau \). Furthermore, we let \(e^{n+1/2}:=(e^{n+1} + e^n)/2\) and \(u^{n+1/2}:=(u^{n+1} + u^n)/2\). For time derivatives at fixed time \(t^n\), we also write \(\partial _t u^n:= \partial _t u(t^n,\cdot )\).

We begin by establishing a differential equation for the time discrete error \(e^n = u^n-u^n_\tau \). This is stated in the following lemma.

Lemma 3.1

(Consistency error) The error \(e^n = u^n-u_\tau ^n\) fulfills the identity

$$\begin{aligned} \mathrm {i}D_\tau e^n+\varDelta e^{n+1/2}-Ve^{n+1/2}-e^n_\gamma = T^n, \end{aligned}$$
(3.3)

where the the consistency error \(T^n\) is given by

$$\begin{aligned} T^n:= & {} \mathrm {i}\,( D_\tau u^n - \partial _t u(t_{n+1/2}))+\varDelta ( u^{n+1/2}-u(t_{n+1/2})) -V(u^{n+1/2}-u(t_{n+1/2})) \nonumber \\&-\,\bigg (\gamma (\xi ^n) -\gamma (|u(t_{n+1/2})|^2) \bigg ). \end{aligned}$$
(3.4)

Here, \(e^n_\gamma := \gamma (\xi ^n)u^{n+1/2}-\gamma (\xi ^n_\tau )u^{n+1/2}_\tau \) for some bounded functions \(\xi ^n,\xi ^n_\tau \in L^{\infty }(\mathscr {D})\) with the properties that

$$\begin{aligned}&\xi ^n(x) \in [\min (|u^n|^2,|u^{n+1}|^2),\max (|u^n|^2,|u^{n+1}|^2)] \quad \hbox {and} \\&\quad \xi ^n_\tau (x) \in [\min (|u^n_\tau |^2,|u^{n+1}_\tau |^2),\max (|u^n_\tau |^2,|u^{n+1}_\tau |^2)] \end{aligned}$$

for almost all \(x\in \mathscr {D}\).

Proof

It is easily verified that exact solution fulfills

$$\begin{aligned} \mathrm {i}D_\tau u^n + \varDelta u^{n+1/2}- Vu^{n+1/2} -\frac{\varGamma (|u^{n+1}|^2)-\varGamma (|u^n|^2)}{|u^{n+1}|^2-|u^n|^2} = T^n. \end{aligned}$$
(3.5)

By the regularity assumptions we can apply Taylor expansion arguments to \(T^n\) to see:

$$\begin{aligned} \sum _{k=0}^{N}\Vert T^k\Vert ^2 \le C\tau ^3 \end{aligned}$$
(3.6)

The argument that proves (3.6) is elaborated in “Appendix A”, where it becomes also visible how the regularity assumptions enter explicitly in the estimate. Next, subtracting (3.5) from (3.1) we find that \(e^n = u^n-u_\tau ^n\) satisfies:

$$\begin{aligned}\mathrm {i}D_\tau e^n+\varDelta e^{n+1/2}-Ve^{n+1/2}-e^n_\gamma = T^n \end{aligned}$$

where \(e^n_\gamma \) denotes the error coming from the nonlinear term, defined by

$$\begin{aligned} e_\gamma ^n= \frac{\varGamma (|u^{n+1}|^2)-\varGamma (|u^n|^2)}{|u^{n+1}|^2-|u^n|^2}u^{n+1/2}-\frac{\varGamma (|u_\tau ^{n+1}|^2)-\varGamma (|u^n_\tau |^2)}{|u^{n+1}_\tau |^2-|u^n_\tau |^2} u_\tau ^{n+1/2}. \end{aligned}$$

Recalling the definition of \(\varGamma \) we have:

$$\begin{aligned} \frac{\varGamma (|u^{n+1}|^2)-\varGamma (|u^n|^2)}{|u^{n+1}|^2-|u^n|^2} = \frac{1}{|u^{n+1}|^2-|u^n|^2}\int _{|u^n|^2}^{|u^{n+1}|^2} \gamma (r)\,dr =:\gamma (\xi ^n) , \end{aligned}$$

likewise

$$\begin{aligned} \frac{\varGamma (|u^{n+1}_\tau |^2)-\varGamma (|u^n_\tau |^2)}{|u^{n+1}_\tau |^2-|u^n_\tau |^2} = \frac{1}{|u^{n+1}_\tau |^2-|u^n_\tau |^2}\int _{|u^n_\tau |^2}^{|u^{n+1}_\tau |^2} \gamma (r)\,dr =: \gamma (\xi ^n_\tau ). \end{aligned}$$

The expression for \(e^n_\gamma \) is thus simplified to

$$\begin{aligned} e^n_\gamma := \gamma (\xi ^n)u^{n+1/2}-\gamma (\xi ^n_\tau )u^{n+1/2}_\tau , \end{aligned}$$

where \(\xi ^n\) is a function taking values between \(|u^n|^2\) and \(|u^{n+1}|^2\) and \(\xi ^n_\tau \) takes values between \(|u^n_\tau |^2\) and \(|u_{\tau _{}}^{n+1}|^2\). \(\square \)

The differential equation in Lemma 3.1 is now used to derive a recurrence formula for the \(H^1\)-norm of the error. Multiplying (3.3) by \(D_\tau e^n\), integrating and taking the real part yields:

$$\begin{aligned} \frac{\Vert \nabla e^{n+1}\Vert ^2-\Vert \nabla e^n\Vert ^2}{2\tau } = \underbrace{ \text {Re}( \langle e_\gamma ^n, D_\tau e^n\rangle )}_{\mathrm{I}} \ \underbrace{-\,\text {Re}(\langle T^n,D_\tau e^n\rangle ).}_{\mathrm{II}} \end{aligned}$$
(3.7)

The idea is to bound the terms I and II in such a way that Grönwall’s inequality can be used. We proceed to bound term I. Multiplying the error PDE (3.3) by \(e^n_\gamma \) results in:

$$\begin{aligned} \mathrm {i}\langle D_\tau e^n, e_\gamma ^n\rangle = \langle \nabla e^{n+1/2},\nabla e_\gamma ^n\rangle + \langle Ve^{n+1/2},e^n_\gamma \rangle + \Vert e_\gamma ^n\Vert ^2 + \langle T^n,e_\gamma ^n\rangle \end{aligned}$$

and consequently

$$\begin{aligned} |\mathrm{I}|= & {} |\text {Re}( \langle D_\tau e^n, e_\gamma ^n\rangle )| \nonumber \\\le & {} |\text {Im}(\langle \nabla e^{n+1/2}, \nabla e_\gamma ^n)\rangle |+|\text {Im}\langle Ve^{n+1/2},e^n_\gamma \rangle |+ |\text {Im}(\langle T^n, e_\gamma ^n\rangle )| \nonumber \\\lesssim & {} \Vert \nabla e^{n+1/2} \Vert ^2+ \Vert \nabla e^n_\gamma \Vert ^2+ \Vert V\Vert _\infty (\Vert e^{n+1/2}\Vert ^2+ \Vert e^n_\gamma \Vert ^2)+\tau ^4+\Vert e^n_\gamma \Vert ^2+\Vert T^n\Vert ^2 \nonumber \\\lesssim & {} \Vert \nabla e^{n+1}\Vert ^2+\Vert \nabla e^{n}\Vert ^2+\Vert \nabla e^n_\gamma \Vert ^2+\Vert e^n_\gamma \Vert ^2+\Vert T^n\Vert ^2+\tau ^4. \end{aligned}$$
(3.8)

In order to use Grönwall’s inequality we need to bound \(\Vert e^n_\gamma \Vert \) and \(\Vert \nabla e^n_\gamma \Vert \) in terms of \(\Vert e^n\Vert \),\(\Vert \nabla e^n\Vert \) and terms of \(\mathscr {O}(\tau ^2)\). These bounds are formulated in the two following lemmas.

Lemma 3.2

Given the optimal \(L^2\)-convergence of Theorem 3.1 and the uniform bounds (3.2), the error coming from the nonlinear term behaves as \(\tau ^2\), i.e. \(\Vert e^n_\gamma \Vert \lesssim \tau ^2 \).

Proof

We introduce the function f by

$$\begin{aligned} f(a,b) = \frac{1}{b-a}\int _{a}^{b} \gamma (r)\,dr. \end{aligned}$$
(3.9)

The partial derivatives of f with respect to the i’th variable are denoted by \(\partial _{i}f\). A standard application of the mean value theorem yields that

$$\begin{aligned} e^n_\gamma= & {} (\gamma (\xi ^n)-\gamma (\xi ^n_\tau ))u^{n+1/2} + \gamma (\xi ^n_\tau )(u^{n+1/2}-u^{n+1/2}_\tau ) \nonumber \\= & {} (f(|u^n|^2,|u^{n+1}|^2)-f(|u^n_\tau |^2,|u^{n+1}_\tau |^2) ) u^{n+1/2}+\gamma (\xi ^n_\tau )e^{n+1/2} \nonumber \\= & {} \big (\partial _1 f(\eta ^n,\eta ^{n+1})(|u^n|^2-|u^n_\tau |^2)+\partial _2 f(\eta ^n,\eta ^{n+1})(|u^{n+1}| ^2-|u^{n+1}_\tau |^2) \big )u^{n+1/2}\nonumber \\&\qquad +\gamma (\xi ^n_\tau )e^{n+1/2}\nonumber \end{aligned}$$

for some functions \(\eta ^{n}\) and \(\eta ^{n=1}\) with

$$\begin{aligned} \min \{ |u^n(x)|^2 , |u^n_{\tau }(x)|^2 \} \le \eta ^n(x) \le \max \{|u^{n}(x)|^2,|u^n_{\tau }(x)|^2 \}. \end{aligned}$$
(3.10)

A quick sanity check shows that the partial derivatives of f are bounded by the derivative of \(\gamma \):

$$\begin{aligned} |\partial _1 f(a,b)|= & {} | \frac{1}{b-a}(\gamma (c)-\gamma (a))| = |\gamma '(\theta ^-) \frac{c-a}{b-a}| \nonumber \\ |\partial _2 f(a,b)|= & {} |\frac{1}{b-a}(\gamma (b)-\gamma (c)) | =| \gamma '(\theta ^+) \frac{b-c}{b-a}| , \end{aligned}$$
(3.11)

where c and \(\theta \) lie somehwere between a and b. Hence

$$\begin{aligned} \Vert e^n_\gamma \Vert\le & {} {} \Vert \partial _1 f(\eta ^n,\eta ^{n+1})\Vert _{L^\infty } \Vert u^{n+1/2}\Vert _{L^\infty }\Vert |u^n|^2-|u^n_\tau |^2\Vert \\&+\, \Vert \partial _2 f(\eta ^n,\eta ^{n+1})\Vert _{L^\infty } \Vert u^{n+1/2}\Vert _{L^\infty }\Vert |u^{n+1}|^2-|u^{n+1}_\tau |^2\Vert \\&+\, \Vert \gamma (\xi ^n_\tau )\Vert _{L^\infty } \Vert e^{n+1/2}\Vert . \end{aligned}$$

Since \(\Vert |u^n|^2-|u^n_\tau |^2 \Vert \le \Vert | u^{n}|+|u^{n+1}\Vert | _{L^{\infty }}\Vert | e^{n}\Vert |\) it now follows with the max-norm bounds for the exact solution and the semi-discrete solution that \( \Vert e^n_\gamma \Vert \lesssim \Vert e^{n+1}\Vert +\Vert e^{n} \Vert \lesssim \tau ^2\).\(\square \)

Lemma 3.3

Given Theorem 3.1, the gradient of the error coming from the nonlinear term is bounded as \(\Vert \nabla e^n_\gamma \Vert \lesssim \Vert \nabla e^{n+1} \Vert +\Vert \nabla e^n\Vert +\tau ^2.\)

Proof

The steps are much the same as in the previous lemma, with the exception that we need to use max-norm estimates for the gradient of the exact solution, which are however not available for the numerical approximation. We begin with a suitable error splitting of the form

$$\begin{aligned} \nabla e^n_\gamma= & {} \nabla [\gamma (\xi ^n)u^{n+1/2}-\gamma (\xi ^n_\tau ) u^{n+1/2}_\tau ] \\= & {} \nabla \gamma (\xi ^n)u^{n+1/2}+\gamma (\xi ^n)\nabla u^{n+1/2} -\nabla \gamma (\xi ^n_\tau )u^{n+1/2}_\tau -\gamma (\xi ^n_\tau )\nabla u^{n+1/2}_\tau \\= & {} u^{n+1/2}_\tau \nabla (\gamma (\xi ^n)-\gamma (\xi ^n_\tau ))+\nabla \gamma (\xi ^n)(u^{n+1/2}-u^{n+1/2}_\tau )\\&\quad + (\gamma (\xi ^n)-\gamma (\xi ^n_\tau ))\nabla u^{n+1/2} + \gamma (\xi ^n_\tau )\nabla (u^{n+1/2}-u^{n+1/2}_\tau ). \end{aligned}$$

Using the previous lemma and the max-norm bounds for the exact solution we may conclude,

$$\begin{aligned}&\Vert \nabla e^n_\gamma \Vert \le \Vert u^{n+1/2}_\tau \Vert _{L^\infty } \Vert \nabla (\gamma (\xi ^n)-\gamma (\xi ^n))\Vert + \Vert \nabla \gamma (\xi ^n)\Vert _{L^\infty }\Vert e^{n+1/2}\Vert \nonumber \\&\quad +\Vert \nabla u^{n+1/2}\Vert _{L^\infty }\tau ^2 +\Vert \gamma (\xi ^n_\tau )\Vert _{L^\infty }\Vert \nabla e^{n+1/2}\Vert . \end{aligned}$$
(3.12)

What is left to bound is the term \(\Vert \nabla (\gamma (\xi ^n)-\gamma (\xi ^n))\Vert \). We consider its dependence on \(|u^n|^2,|u^{n+1}|^2,|u^n_\tau |^2\) and \(|u^{n+1}_\tau |^2\) and find:

$$\begin{aligned} \nabla (\gamma (\xi ^n)-\gamma (\xi ^n_\tau ))= & {} \nabla (f(|u^n|^2,|u^{n+1}|^2)-f(|u^n_\tau |^2,|u^{n+1}_\tau |^2)) \nonumber \\= & {} \partial _1 f^n \nabla |u^n|^2+\partial _2 f^n \nabla |u^{n+1}|^2 - \partial _1 f_{\tau }^n \nabla |u^n_\tau |^2- \partial _2 f_{\tau }^n \nabla |u^{n+1}_\tau |^2 \nonumber \\= & {} (\partial _1 f^n -\partial _1 f^n_{\tau }) \nabla |u^n|^2 + \partial _1 f_{\tau }^n \nabla (|u^n|^2-|u^n_\tau |^2) \nonumber \\&\quad + (\partial _2 f^n -\partial _2 f^n_{\tau })\nabla |u^{n+1}|^2+\partial _2 f_{\tau }^n \nabla (|u^{n+1}|^2-|u^{n+1}_\tau |^2), \end{aligned}$$
(3.13)

where \(\partial _i f^n := \partial _i f(|u^n|^2,|u^{n+1}|^2) \text{ and } \partial _i f^n_{\tau }:= f_i(|u^n_\tau |^2,|u^{n+1}_\tau |^2)\). Another application of the mean value theorem yields:

$$\begin{aligned} \partial _1 (f^n - f_{\tau }^n)= & {} \partial _{1,1}f(\theta ^n,\theta ^{n+1}) (|u^n|^2-|u^n_\tau |^2) + \partial _{1,2}f(\theta ^n,\theta ^{n+1}) (|u^{n+1}|^2-|u^{n+1}_\tau |^2), \nonumber \\ \partial _2 (f^n - f_{\tau }^n)= & {} \partial _{2,1}f(\vartheta ^n,\vartheta ^{n+1}) (|u^n|^2-|u^n_\tau |^2) + \partial _{2,2}f(\vartheta ^n,\vartheta ^{n+1})(|u^{n+1}|^2-|u^{n+1}_\tau |^2), \nonumber \\ \end{aligned}$$
(3.14)

for suitable mean value functions that can be bounded pointwise by the maximum of the exact solution and the semi-discrete CN-approximation. The following quick calculations show that the second partial derivatives of f are bounded by the second derivative of \(\gamma \):

$$\begin{aligned} \partial _{1,1}f(a,b)= & {} \frac{1}{(b-a)^2}2(\gamma (c)-\gamma (a))-\frac{\gamma '(a)}{b-a}= \frac{\gamma '(\theta ^-)-\gamma '(a)}{b-a}+C_{\gamma ''}, \nonumber \\ \partial _{1,2}f (a,b)= & {} \frac{1}{(b-a)^2}(\gamma (b)-2\gamma (c)+\gamma (a)) = \gamma ''(c)+C_{\gamma ''},\nonumber \\ \partial _{2,2}f(a,b)= & {} \frac{\gamma '(b)}{b-a}-\frac{1}{(b-a)^2}2(\gamma (b)-\gamma (c))= \frac{\gamma '(b)-\gamma '(\theta ^+)}{b-a}+C_{\gamma ''}, \nonumber \\ \end{aligned}$$
(3.15)

where it was used that \(c = (b+a)/2+ C_{\gamma ''}(a-b)^2\). It thus becomes clear that \(\Vert \partial _{1,1} f+ \partial _{1,2}f+\partial _{2,2}f \Vert _{L^\infty } \le C_{\gamma ''}\). This gives us the following \(L^2\)-bound on (3.13)

$$\begin{aligned} \Vert \nabla (\gamma (\xi ^n)-\gamma (\xi ^n_\tau ))\Vert \le C_{\gamma ''}(\Vert e^n\Vert +\Vert e^{n+1}\Vert )+C_{\gamma '}\Vert \nabla (|u^{n}|^2-|u^{n}_\tau |^2+|u^{n+1}|^2-|u^{n+1}_\tau |^2) \Vert . \end{aligned}$$

We conclude that

$$\begin{aligned} \Vert \nabla e^n_\gamma \Vert \lesssim \Vert \nabla (|u^n|^2-|u^n_\tau |^2)\Vert +\Vert \nabla (|u^{n+1}|^2-|u^{n+1}_\tau |^2)\Vert +\tau ^2. \end{aligned}$$

Here it is noted that \(\nabla (|u^n|^2-|u^n_\tau |^2)\) may be written as

$$\begin{aligned} \nabla (|u^n|^2-|u^n_\tau |^2) = 2\text {Re}( (u^n-u^n_\tau )\nabla \overline{u}^n+u^n_\tau \nabla (\overline{u}^n-\overline{u}^n_\tau )). \end{aligned}$$

Hence, using the max-norm bounds for the gradient of the exact solution we have

$$\begin{aligned} ||\nabla e^n_\gamma || \lesssim ||\nabla e^{n+1}||+||\nabla e^n||+\tau ^2. \end{aligned}$$

\(\square \)

With Lemmas 3.2 and 3.3 we now have the following bound on term I.

Lemma 3.4

For term \(\mathrm{I}\) which is given by (3.8), we have the estimate

$$\begin{aligned} |\mathrm{I}|\, \lesssim \, \Vert \nabla e^{n+1} \Vert ^2+\Vert \nabla e^n \Vert ^2 + \Vert T^n\Vert ^2+\tau ^4. \end{aligned}$$
(3.16)

We can now proceed to bound term II. Here we explicate the Taylor term using (3.4) to see

$$\begin{aligned} \mathrm{II}= & {} - \text {Re}(\langle T^n , D_\tau e^n \rangle ) \\\le & {} \underbrace{ | \langle ( D_\tau u^n - \partial _t u(t_{n+1/2})),D_\tau e^n\rangle |}_{\mathrm{IIa}}+\underbrace{\text {Re}\langle -\varDelta ( u^{n+1/2}-u(t_{n+1/2})),D_\tau e^n\rangle }_{\mathrm{IIb}} \\&+\,\underbrace{|\langle V(u^{n+1/2}-u(t_{n+1/2})),D_\tau e^n\rangle |}_{\mathrm{IIc}} \\&+\,\underbrace{|\langle \bigg ( \frac{\varGamma (|u^{n+1}|^2)-\varGamma (|u^n|^2)}{|u^{n+1}|^2-|u^n|^2}-\gamma (|u(t_{n+1/2})|^2),D_\tau e^n\rangle \bigg )|}_{\mathrm{IId}}. \end{aligned}$$

We start with estimating \(\mathrm{IIa},\mathrm{IIc}\) and \(\mathrm{IId}\), which can be bounded in a similar way.

Step 1, bounding \(\mathrm{IIa}\):

By replacing \(D_\tau e^n\) using (3.3) (i.e. time derivative is replaced by regularity in space) we have

$$\begin{aligned}&{|\langle D_\tau u^n-\partial _t u(t_{n+1/2}) , D_\tau e^n \ \rangle |}\\&\quad \le |\langle D_\tau u^n-\partial _t u(t_{n+1/2}) ,\varDelta e^{n+1/2}\rangle |+|\langle D_\tau u^n-\partial _t u(t_{n+1/2}) ,V e^{n+1/2}\rangle | \\&\qquad +\, | \langle D_\tau u^n-\partial _t u(t_{n+1/2}) , e_\gamma ^n\rangle | + | \langle D_\tau u^n-\partial _t u(t_{n+1/2}) ,T^n \rangle | \\&\quad \overset{(3.6)}{\lesssim } \Vert \nabla e^{n+1}\Vert ^2 + \Vert \nabla e^n\Vert ^2 + \Vert D_\tau u^n-\partial _t u(t_{n+1/2})\Vert _{H^1}^2 + \Vert T^n\Vert ^2 + \tau ^4 . \end{aligned}$$

Here we see how the assumption that \(\partial _{ttt}u\in L^2(0,T;H^1(\mathscr {D}))\) will enter, as it allows us to conclude that \(\sum \Vert D_\tau u^n-\partial _t u(t_{n+1/2})\Vert _{H^1}^2\lesssim \tau ^3\).

Step 2, bounding \(\mathrm{IIc}\):

We use the same idea as for \(\mathrm{IIa}\) to get

$$\begin{aligned}&{|\langle V(u^{n+1/2}-u(t_{n+1/2})),D_\tau e^n\rangle | \le \Vert V\Vert _{L^\infty }|\langle u^{n+1/2}-u(t_{n+1/2}),D_\tau e^n\rangle |} \\&\quad \le \Vert V\Vert _{L^\infty }\big ( |\langle u^{n+1/2}-u(t_{n+1/2}) ,\varDelta e^{n+1/2}\rangle |+\\&\qquad +\,|\langle u^{n+1/2}-u(t_{n+1/2}) ,V e^{n+1/2}\rangle | +| \langle u^{n+1/2}-u(t_{n+1/2}) , e_\gamma ^n\rangle | \\&\qquad + \, | \langle u^{n+1/2}-u(t_{n+1/2}) ,T^n \rangle | \big ) \\&\quad \lesssim \Vert \nabla e^{n+1}\Vert ^2 + \Vert \nabla e^n\Vert ^2 + \Vert u^{n+1/2}-u(t_{n+1/2})\Vert _{H^1}^2+ \Vert T^n\Vert ^2+ \tau ^4. \end{aligned}$$

Step 3, bounding IId:

We start from

$$\begin{aligned}&\big |\big \langle \frac{\varGamma (|u^{n+1}|^2)-\varGamma (|u^n|^2)}{|u^{n+1}|^2-|u^n|^2}-\gamma (|u(t_{n+1/2})|^2),D_\tau e^n\big \rangle \big | \\&\quad = |\langle \gamma (\xi ^n)-\gamma (|u(t_{n+1/2})|^2), D_\tau e^n\rangle | \end{aligned}$$

replace \(D_\tau e^n\) again using (3.3). Furthermore, in virtue of the assumptions it holds that \(\Vert \nabla (\gamma (\xi ^n)-\gamma (|u(t_{n+1/2})|^2))\Vert \lesssim \Vert u^{n+1/2}-u (t_{n+1/2})\Vert _{H^{1}}\). This is made explicit in the appendix. We thus obtain

$$\begin{aligned} |\langle \gamma (\xi ^n)-\gamma (|u(t_{n+1/2})|^2), D_\tau e^n\rangle |= & {} |\langle \gamma (\xi ^n)-\gamma (|u(t_{n+1/2})|^2), -\varDelta e^{n+1/2} + Ve^{n+1/2}+e^n_\gamma +T^n \rangle | \\\le & {} |\langle \nabla (\gamma (\xi ^n)-\gamma (|u(t_{n+1/2})|^2)),\nabla e^{n+1/2}\rangle | +\Vert T^n\Vert ^2 + \mathscr {O}(\tau ^4) \\\lesssim & {} \Vert \nabla e^{n+1}\Vert ^2+\Vert \nabla e^n\Vert ^2+ \Vert T^n\Vert ^2+\Vert u^{n+1/2}-u(t_{n+1/2})\Vert _{H^1}^2+ \tau ^4. \end{aligned}$$

Step 4, bounding \(\mathrm{IIb}\):

The previous technique does not work on this term since replacing the discrete time derivative with regularity in space would give rise to the term \(\nabla \varDelta (u^{n+1/2}-u(t_{n+1/2}))\), which we can not afford. Instead we use summation by parts in time to get the factor \(D_\tau \varDelta (u^{n+1/2}-u(t_{n+1/2}))\), which when integrated against \(e^{n+1/2}\) can be handled. First we recall:

$$\begin{aligned}D_\tau [a^k b^k] = a^kD_\tau b^k + b^{k+1}D_\tau a^k \quad \Longleftrightarrow \quad -a^kD_\tau b^k = D_\tau [a^k] b^{k+1} - D_\tau [a^kb^k]. \end{aligned}$$

Using this on term \(\mathrm{IIb}\) yields:

$$\begin{aligned}&{\langle -\varDelta (u^{n+1/2}-u), D_\tau e^n\rangle {= \langle \nabla (u^{n+1/2}-u),D_\tau \nabla e^n \rangle }} \\&\quad = D_\tau [\langle \nabla (u^{n+1/2}-u(t_{n+1/2})), \nabla e^{n}\rangle ] - \langle D_\tau \nabla (u^{n+1/2}-u(t_{n+1/2})),\nabla e^{n+1}\rangle \\&\quad \le D_\tau [\langle \nabla (u^{n+1/2}-u(t_{n+1/2})), \nabla e^{n}\rangle ] + |\langle D_\tau \nabla (u^{n+1/2}-u(t_{n+1/2})),\nabla e^{n+1}\rangle | \\&\quad \le D_\tau [\langle \nabla (u^{n+1/2}-u(t_{n+1/2})), \nabla e^{n}\rangle ] + \Vert D_\tau (u^{n+1/2}-u(t_{n+1/2}))\Vert _{H^1}^2 + \Vert \nabla e ^{n+1}\Vert ^2. \end{aligned}$$

Collecting the estimates we have the following estimate for term \(\mathrm{II}\).

Lemma 3.5

For term \(\mathrm{II} = -\text {Re}(\langle T^n , D_\tau e^n \rangle )\) it holds the estimate

$$\begin{aligned} \mathrm{II}\le & {} D_\tau [\langle \nabla (u^{n+1/2}-u(t_{n+1/2})), \nabla e^{n}\rangle ] + C \big ( \Vert \nabla e ^{n+1}\Vert ^2 +\Vert \nabla e^n\Vert ^2 + \tau ^4 +\Vert T^n\Vert ^2 \nonumber \\&+\, \Vert D_\tau u^n-\partial _t u(t_{n+1/2})\Vert _{H^1}^2 +\Vert u^{n+1/2}-u(t_{n+1/2})\Vert _{H^1}^2\big ). \nonumber \\ \end{aligned}$$
(3.17)

Here we note the importance of not estimating the absolute value of the first term since it is necessary to use the fact that n of these terms cancel when summed up, i.e. \(\sum _k D_\tau a^k = \frac{1}{\tau }(a^{n+1}-a^0)\). We are now ready to finish the proof of the first main result.

Proof of Theorem 3.2

We pick off where we left (3.7) and find by using Lemmas 3.4 and 3.5:

$$\begin{aligned}&\frac{\Vert \nabla e^{n+1}\Vert ^2 -\Vert \nabla e^n\Vert ^2}{2\tau }\\&\quad \le D_\tau [\langle \nabla (u^{n+1/2}-u(t_{n+1/2})), \nabla e^{n}\rangle ] +C\big (\Vert \nabla e^{n+1}\Vert ^2 +\Vert \nabla e^n\Vert ^2 +\tau ^4 \\&\qquad +\, \Vert D_\tau u^n -\partial _t u(t_{n+1/2})\Vert ^2_{H^1}+ \Vert T^n\Vert ^2 + \Vert u^{n+1/2}-u(t_{n+1/2})\Vert ^2_{H^1} \big ). \end{aligned}$$

Summing this up and using \(e^0=0\) gives

$$\begin{aligned} \frac{\Vert \nabla e^{n+1}\Vert ^2}{2\tau }&\le C\left( \sum _{k=0}^n \Vert \nabla e^{k}\Vert ^2 \right) + \frac{1}{\tau }\langle \nabla (u^{n+3/2}-u(t_{n+3/2})),\nabla e^{n+1} \rangle \\&\quad +\, C\tau ^3 + \sum _{k=0}^{n} \Vert T^k\Vert ^2+\Vert u^{k+1/2}-u(t_{k+1/2}) \Vert ^2_{H^1}+ \Vert D_\tau u^k\\&\quad -\,\partial _t u(t_{k+1/2})\Vert ^2_{H^1} \end{aligned}$$

and therefore, recalling (3.6),

$$\begin{aligned} \Vert \nabla e^{n+1}\Vert ^2&\le C\left( \sum _{k=0}^n \tau \Vert \nabla e^{k}\Vert ^2 \right) + C \tau ^4 +|\langle \nabla (u^{n+3/2}-u(t_{n+3/2})),\nabla e^{n+1} \rangle |. \end{aligned}$$

Young’s inequality with \(\epsilon >0\) is used on the last term:

$$\begin{aligned} |\langle \nabla (u^{n+3/2}-u(t_{n+3/2})),\nabla e^{n+1} \rangle | \le C\left( \frac{\tau ^4}{\epsilon } + \epsilon \Vert \nabla e^{n+1}\Vert ^2\right) . \end{aligned}$$
(3.18)

Which holds since,

$$\begin{aligned} \Vert \nabla (u^{n+3/2}-u(t_{n+3/2}))\Vert \lesssim \tau ^2 \Vert \partial _{tt}u\Vert _{L^\infty (H^1)}, \end{aligned}$$

where we have \( \Vert \partial _{tt}u\Vert _{L^\infty (H^1)} \lesssim \Vert \partial _{tt}u\Vert _{L^2(H^1)} + \Vert \partial _{ttt}u\Vert _{L^2(H^1)}\) by Sobolev embeddings. Finally we arrive at

$$\begin{aligned} \Vert \nabla e^{n+1}\Vert ^2 \le C\left( \sum _{k=0}^n \tau \Vert \nabla e^k\Vert ^2 \right) + C \tau ^4 + \frac{\tau ^4}{\epsilon } + \epsilon \Vert \nabla e^{n+1}\Vert ^2 \end{aligned}$$

and for e.g. \(\epsilon =1/2\) we can absorb \( \epsilon \Vert \nabla e^{n+1}\Vert ^2\) in the left hand side and conclude

$$\begin{aligned} \Vert \nabla e^{n+1}\Vert ^2 \le C \left( \tau ^4 + \sum _{k=0}^n \tau \Vert \nabla e^k\Vert ^2 \right) . \end{aligned}$$

Grönwall’s inequality now yields:

$$\begin{aligned} \Vert \nabla e^{n+1}\Vert \lesssim \tau ^2. \end{aligned}$$
(3.19)

\(\square \)

4 Fully-discrete Crank-Nicolson scheme

We shall now consider the fully-discrete setting that is based on a finite element discretization in space. For that, we let \(S_h\subset H^1_0(\mathscr {D})\) denote the space of P1 Lagrange finite elements on a quasi-uniform simplicial mesh on \(\mathscr {D}\) with mesh size h. In this setting we have by standard finite element theory (cf. [39]) the following inequality for any \(u \in H^1_0(\mathscr {D})\cap H^2(\mathscr {D})\):

$$\begin{aligned} \Vert u-P_h u\Vert +h\Vert \nabla (u-P_h u) \Vert \lesssim h^2\Vert u\Vert _{H^2}. \end{aligned}$$
(4.1)

Here \(P_h: H^1_0(\mathscr {D}) \rightarrow S_h\) denotes (for example) the Ritz-projection into the finite element space. For a given discrete initial value \(u_h^0\in S_h\) the CN-FEM approximation \(u^{n}_{\tau ,h} \in S_h\) to \(u^{n}\) is given by the fully discrete equation

$$\begin{aligned} \mathrm {i}\left\langle \frac{u^{n+1}_{\tau ,h}-u^n_{\tau ,h}}{\tau },v\right\rangle&= \big \langle \nabla u^{n+1/2}_{\tau ,h},\nabla v \big \rangle + \big \langle Vu^{n+1/2}_{\tau ,h},v\big \rangle \nonumber \\&\quad +\, \left\langle \frac{\varGamma (|u^{n+1}_{\tau ,h}|^2)-\varGamma (|u^{n}_{\tau ,h}|^2)}{|u^{n+1}_{\tau ,h}|^2-|u^n_{\tau ,h}|^2}u^{n+1/2}_{\tau ,h},v\right\rangle \end{aligned}$$
(4.2)

for all \(v\in S_h\). The initial value is selected as \( u^0_{\tau ,h}=P_hu^0\). As in the semi-discrete case, the discretization conserves both mass and energy, i.e.

$$\begin{aligned} \int _{\mathscr {D}} | u^n_{\tau ,h}|^2 \, dx = \int _{\mathscr {D}} |u^0_{\tau ,h}|^2 \, dx \quad \hbox {and} \quad E[u^n_{\tau ,h}] = E[u^0_{\tau ,h}] \qquad \hbox {for all } n\ge 0. \end{aligned}$$

The scheme is well-posed and the corresponding approximations converge in the \(L^{\infty }(L^2)\)-norm with optimal order in space and time to the exact solution. A proof of this statement can be easily extracted from [35, Theorem 3.1 and Lemma 5.3]. In particular, we have the following result.

Theorem 4.1

Under the general assumptions of this paper, there exists a solution \(u^n_{\tau ,h}\in H^1_0(\mathscr {D})\) to the fully discrete Crank-Nicolson scheme (4.2) such that the following a priori error estimates hold

$$\begin{aligned} \sup _{0\le n \le N} \Vert u^n_{\tau ,h} - u^n_\tau \Vert _{L^{2}(\mathscr {D})} \lesssim h^2 \quad \hbox {and} \quad \sup _{0\le n \le N} \Vert u^n_{\tau ,h} - u^n \Vert _{L^{2}(\mathscr {D})} \lesssim h^2+\tau ^2. \end{aligned}$$

With this we are ready to state our final theorem.

Theorem 4.2

(Optimal \(H^1\)-error estimates for the fully discrete method) Let \(u_{\tau _{}}^{n}\in H^1_0(\mathscr {D})\) denote the fully-discrete Crank-Nicolson approximations from Theorem 4.1, then it holds

$$\begin{aligned} \sup _{0\le n \le N} \Vert u^n_{\tau ,h} - u^n \Vert _{H^{1}(\mathscr {D})} \lesssim \tau _{}^2+h. \end{aligned}$$

Proof

First, we recall the inverse estimate on quasi-uniform meshes (cf. [39]), i.e. \(\Vert \nabla v_h \Vert \le C h^{-1} \Vert v_h \Vert \) for all \(v_h \in S_h\), which implies

$$\begin{aligned} \Vert \nabla (P_h(u_\tau ^n)-u_{\tau ,h}^n) \Vert \le C h^{-1} \Vert P_h(u_\tau ^n)-u_{\tau ,h}^n \Vert . \end{aligned}$$
(4.3)

With this, the \(H^1\) convergence result (3.19) together with Theorem 4.1 suffice to show optimal \(H^1\)-convergence rates for the fully discrete method. This is made clear by the following splitting.

$$\begin{aligned} \Vert \nabla (u^n-u_{\tau ,h}^n)\Vert&\le \Vert \nabla (u^n-u_\tau ^n)\Vert +\Vert \nabla (u_\tau ^n-P_h(u_\tau ^n))\Vert +\Vert \nabla (P_h(u_\tau ^n)-u_{\tau ,h}^n)\Vert \\&\overset{(4.3)}{\le } \Vert \nabla (u^n-u_\tau ^n)\Vert + Ch+Ch. \end{aligned}$$

Here we have made use of the inequality (4.1), the uniform \(H^2\)-regularity of \(u^n_\tau \), i.e. \(\Vert u^n_\tau \Vert _{H^2} \le C(u)\) (cf. (3.2)) and the optimal \(L^2\)-estimates. In virtue of Theorem 3.2 we may thus conclude:

$$\begin{aligned} \Vert \nabla (u-u^n_{\tau ,h}) \Vert \le C(\tau ^2+h). \end{aligned}$$
(4.4)

\(\square \)

Detailed numerical studies that confirm the optimal convergence rates stated in Theorem 4.1 and Theorem 4.2 are presented in [21, 35].

5 Implementation and numerical examples

In this section we will discuss how the Crank-Nicolson FEM discretization can be efficiently implemented and practically used. Afterwards, we present two numerical experiments. The first one is to confirm the theoretically predicted convergence rates in Theorem 4.1 and the second experiment demonstrates that our approach is fully competitive in low regularity regimes, where we compare it with a time-splitting spectral method.

5.1 Efficient implementation

The Crank Nicolson method (4.2), albeit popular, suffers from the drawback that it requires solving a fully nonlinear system of equations in each time step. Furthermore, this system of equations is often solved through a Newton step, the implementation of which can become complicated and expensive for general nonlinearities. We present here a competitive fixed point solver which makes the method perform on par in terms of computational time with linearized time-stepping methods such as the RE-FEM proposed by C. Besse [23] which was found to be best performing in [21].

To detail the proposed fixed-point iteration, let \(U^n \in \mathbb {R}^N\) denote the vector of nodal values that belongs to the function \(u^n_{\tau ,h} \in S_h\). Introducing the following matrix notation:

$$\begin{aligned}&(M)_{ij} = \langle v_j , v_i \rangle&\quad A_{ij} = \langle \nabla v_j , \nabla v_i \rangle \\&(M_V)_{ij} = \langle V v_j , v_i \rangle&\quad (M_\varGamma ) _{ij}(U^{n+1},U^n) = \langle \gamma (\xi ^{n+1/2}_{\tau ,h}) v_j , v_i \rangle \end{aligned}$$

the equation (4.2) in matrix form becomes :

$$\begin{aligned} \mathrm {i}M\frac{U^{n+1}-U^n}{\tau } =\big ( A+M_V+M_\varGamma (U^{n+1},U^n) \big )\frac{U^{n+1}+U^n}{2}. \end{aligned}$$

Let \(L_1 = M+\mathrm {i}\tau /2(A+M_V)\) and \(L_2 = M-\mathrm {i}\tau /2(A+M_V)\). Our fixed point iteration takes the form:

$$\begin{aligned} U^{n+1}_{i+1} = L_1^{-1}L_2U^n + L_1^{-1}M_\varGamma (U^{n+1}_i,U^n)(U^{n+1}_i+U^n)/2. \end{aligned}$$
(5.1)

Here we note that matrix \(L_1\) does not change with time. Hence, the above iteration can be done efficiently by precomputing the LU-factorization of \(L_1\). After it is precomputed, each time step only involves matrix-vector multiplications, but no longer the solving of a linear system of equations. In fact, the main cost in each time step account for the assembly of the updated mass matrix \(M_\varGamma (U^{n+1}_i,U^n)\) with the densities from the previous time step and the previous iteration. Compared to this, all other costs are essentially negligible. We find that typically, but dependent on \(\tau \), 4–8 iterations are required to reach a tolerance of machine epsilon. To illustrate the efficiency we conclude with two numerical test problems.

5.2 Harmonic potential

First we consider a smooth potential to confirm the expected convergences rates for a type of nonlinearity that complements previous test cases [21, 35]. Here we seek u(x, t) with

$$\begin{aligned} {\left\{ \begin{array}{ll} \mathrm {i}\partial _t u = -\varDelta u + V u + \gamma (|u|^2)u&{}{} \quad \text{ in } \mathscr {D} \times (0,T], \\ u(\cdot ,t) = 0 &{}{}\quad \text{ on } \partial \mathscr {D} \times (0,T], \\ u(\cdot ,0) = u_0&{}{}\quad \text{ in } \mathscr {D} , \end{array}\right. } \end{aligned}$$
(5.2)

where we consider the saturated nonlinearity \(\gamma (r) := r/(1+r)\) (cf. [8,9,10]) . Furthermore, \(\mathscr {D} = [-5,5]^2\) is the computational domain, the maximum time is selected as \(T = 1\) and the trapping potential \(V(x,y) = (\nu _x x)^2+(\nu _y y)^2\). For the time-dependent problem we set the trapping frequencies to \(\nu _x = 2\) and \(\nu _y = 3\). The initial value \(u_0\) is the unique positive ground state with \(\int _{\mathscr {D}} |u_0|^2 =1\) to the problem with \(\nu _x=\nu _y = 1\), i.e. it solves the eigenvalue problem

$$\begin{aligned}\lambda _0 u_0= -\varDelta u_0+Vu_0+\gamma (|u_0|^2)u_0, \end{aligned}$$

with ground state eigenvalue (chemical potential) \(\lambda _0>0\). The \(H^1\)-errors are presented in Table 1. The \(\mathscr {O}(h)\)-convergence is best seen in column \(\tau =2^{-9}\), where initially the convergence is \(\mathscr {O}(h^{1.5})\) but flattens out to \(\mathscr {O}(h^{1.2})\) for the last data point. Since the reference solution is also computed with \(h=0.0125\), it is expected that this last order of convergence is an overestimate. Using the values in row \(h=0.0125\) we estimate the order of convergence with respect to \(\tau \) to be 1.9, hence, confirming the theoretically predicted rates from Theorem 4.1.

Table 1 \(H^1\)-errors for harmonic potential test case (5.2)
Full size table
Fig. 1
figure 1

Plot of the density (\(|u|^2\)) of the reference solution to (5.3) at \(T=1\)

Full size image

5.3 Discontinuous potential

This example illustrates a moderate setting where, due to reduced regularity of the exact solution, finite element based methods are preferable over spectral methods. In the following we compare the Crank-Nicolson approach with a Strang splitting spectral method of order 2 (SP2) [40] which is known to show a very good performance in smooth settings.

In this test problem we seek u(x, t) with

$$\begin{aligned} \left\{ \begin{array}{ll} \mathrm {i}\partial _t u = -\varDelta u + V u + \gamma (|u|^2)u &{}{}\quad \text{ in } \mathscr {D} \times (0,T], \\ u(\cdot ,0) = u_0&{}{}\quad \text{ in } \mathscr {D} , \end{array}\right. \end{aligned}$$
(5.3)

where we consider the saturated nonlinearity \(\gamma (r) := 10r/(1+r)\). For a fair comparison with the SP2, we consider our problem with periodic boundary conditions which are easier to handle by the spectral method. The generalization of the Crank-Nicolson method to periodic boundary conditions is straightforward. Furthermore, \(\mathscr {D} = [-5,5]^2\) is again the computational domain and the maximum time is selected as \(T = 1\). The trapping potential \(V(x,y) = (\nu _x x)^2+(\nu _y y)^2+100(\mathbb {1}_{|x|\ge 1}(x)+\mathbb {1}_{|y|\ge 1}(y))\), with trapping frequencies \(\nu _x = 1\) and \(\nu _y = 3\), is discontinuous and causes a slight loss of regularity. We stress that this is a moderate test case, as illustrated in Fig. 1 most of the dynamics take place within the unit cube where the potential is smooth. The initial value \(u_0\) is the unique positive ground state with \(\int _{\mathscr {D}} |u_0|^2 =1\) and \(V_0(x,y) = x^2+ y^2\), i.e. it solves the eigenvalue problem

$$\begin{aligned}\lambda _0 u_0= -\varDelta u_0+V_0u_0+\gamma (|u_0|^2)u_0. \end{aligned}$$
Table 2 Errors and computational times for the CN-FEM
Full size table
Table 3 Errors and computational times for the SP2
Full size table

The errors and the computational times of the CN-FEM are presented in Table 2 and the errors and computational times of the SP2 in Table 3. The reference solution, \(u_{\text {ref}}\), is computed using the CN-FEM with \(h = 0.0125\) and \(\tau = 2^{-13} \). The implantation was done in Julia. It is important to keep in mind that the SP2 uses mostly inbuilt functions such as the fast Fourier transform from the C subroutine library (FFTW). These functions are heavily optimized and show an extremely good performance. In spite of this we see, comparing the errors of the CN-FEM for \(h = 0.025\) and \(\tau = 2^{-10}\) to those of the SP2 for \(N_{DoF} = 3.2\cdot 10^{6}\) and \(\tau = 2^{-13}\), that they are on par with respect to CPU time relative to accuracy, with a slight computational advantage for the CN-FEM. This advantage becomes clearer, the larger the region of reduced regularity (e.g. in the context of optical lattices or disorder potentials). This justifies the usage of CN-FEM in low regularity regimes. Furthermore, it is clearly seen that the space discretization dominates the error of the SP2, in order for the spectral method to catch up in terms of accuracy with the CN-FEM with \(h=0.0125\), an estimated 10 to 40 million degrees of freedom would be needed and thus the memory cost would become an issue.