A note on optimal $H^1$-error estimates for Crank-Nicolson approximations to the nonlinear Schr\"odinger equation

In this paper we consider a mass- and energy--conserving Crank-Nicolson time discretization for a general class of nonlinear Schr\"odinger equations. This scheme, which enjoys popularity in the physics community due to its conservation properties, was already subject to several analytical and numerical studies. However, a proof of optimal $L^{\infty}(H^1)$-error estimates is still open, both in the semi-discrete Hilbert space setting, as well as in fully-discrete finite element settings. This paper aims at closing this gap in the literature.


Introduction
In this paper we consider nonlinear Schrödinger equations (NLS) seeking a complex function u(t, x) such that i∂ t u = −△u + V u + γ(|u| 2 )u in a bounded domain D ⊂ R d , with a homogenous Dirichlet boundary condition on ∂D and a given initial value. Here, V (x) is a known real-valued potential and γ : [0, ∞) → 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 γ(|u| 2 )u = κ|u| 2 u, for some κ ∈ R.
In this case, the equation is called Gross-Pitaevskii equation. It has applications in optics [1,13], fluid dynamics [30,31] and, most importantly, in quantum physics, where it models for example the dynamics of Bose-Einstein condensates in a magnetic trapping potential [12,20,23]. Another relevant class are saturated nonlinearities, such as γ(|u| 2 ) = κ|u| 2 (1 + α|u| 2 ) −1 for some α ≥ 0, which appear in the context of nonlinear optical wave propagation in layered metallic structures [11,17] or the propagation of light beams in plasmas [22]. 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 [16,25].
For the subclass of power law nonlinearities of the form γ(|u| 2 ) = K k=1 κ k |u| 2σ k for σ k ≥ 0 and α k ∈ R, a mass and energy conserving relaxation scheme was proposed and analyzed by Besse [7,8]. Thanks to its properties, the scheme shows a very good performance in realistic physical setups [16]. Despite the large variety of different numerical approaches for solving the time-dependent NLS (cf. [2,3,5,14,18,19,21,24,27,26,28,29,32] 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 [24] and which is long-known in the physics community, is a Crank-Nicolson-type approach where the nonlinearity is approximated by a suitable difference quotient involving the primitive integral of γ. This is also the time discretization that we shall consider in this paper.
A combination of the method with a finite difference space discretization was proposed and analyzed by Bao and Cai [4,6]. Combining the Crank-Nicolson time discretization with a P 1 finite element discretization in space, the first a priori error estimates for the arising method were obtained by Sanz-Serna 1984 [24] for cubic nonlinearities. He considers the case d = 1 and derives optimal L ∞ (L 2 )-error estimates under the coupling constraint τ h, where τ 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. [2] improved this result by showing optimal convergence rates in L ∞ (L 2 ) in dimension d = 1, 2, 3 and under the relaxed coupling constraint τ h d/4 . Finally, in 2017 [15], the L ∞ (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 and the influence of potentials could be taken into account. However, so far, optimal error estimates in L ∞ (H 1 ) 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. [2,14,18,19,24,32]) 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 ∞ (L 2 )-and/or L ∞ (H 1 )-error estimates. 3. Using inverse estimates in the finite element space to show that the truncated approximations are uniformly bounded in L ∞ (L ∞ ) by a term of the form C(1 + h −s (τ 2 + h p )), with appropriate powers p > s > 0 that depend on the considered space discretization, regularity and space dimension. 4. Concluding that if τ 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 ∞ (H 1 )-error estimates for the Crank-Nicolson FEM. This is because it requires a suitable truncation of the primitive integral of γ 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 ∞ (W 1,∞ ). However, thanks to the new techniques developed in [29] and the CN error analysis suggested in [15] in the context of L ∞ (L 2 )-error estimates, the truncation step is no longer necessary and the desired L ∞ (L ∞ )-bounds can be derived with elliptic regularity theory. With this, it is now possible to obtain L ∞ (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 [15,29] to fill the gap in the literature and prove optimal L ∞ (H 1 )-error estimates for the Crank-Nicolson approach without coupling constraints and for a general class of nonlinearities. The paper is structured as follows. In Section 2 we present the notation and the analytical assumptions on the problem. In Section 3 we present the time-discrete Crank-Nicolson method, we recall its well-posedness and optimal error estimates in L ∞ (L 2 ). Furthermore, we present and prove the new error estimate in L ∞ (H 1 ). The paper concludes with the fully-discrete setting presented in Section 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 ∞ (H 1 )-error estimates, which is the main result of this paper.

Notation and Assumptions
We start with introducing the analytical setting of this work. Throughout the paper we assume that D ⊂ R d (for d = 2, 3) is a convex bounded domain with polyhedral boundary. On D, the Sobolev space of complex-valued, weakly differentiable functions with a zero trace on ∂D and L 2 -integratable partial derivatives is as usual denoted by H 1 0 (D) := H 1 0 (D, C). The potential V ∈ L ∞ (D; 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 [15] for a discussion on this aspect). The (possibly nonlinear) function is assumed to be smooth, fulfills γ(0) = 0 and its growth can be characterized with where L is a function with the following growth properties Note that in [15] the admissible growth condition in 3d requires q ∈ [0, 2), which is however a typo and should be, as above, q ∈ [0, 4) (cf. [10, 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 γ(ρ) = κρ q for κ ≥ 0 and 0 ≤ q < ∞ in 2d and 0 ≤ q < 4 in 3d.
For the initial value we assume that u 0 ∈ H 1 0 (D) ∩ H 2 (D) and, without loss of generality, that it has a normalized mass, i.e. 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 such that u(·, 0) = u 0 and in the sense of distributions. Problem (1) admits at least one solution, that is even unique for repulsive cubic nonlinearities in 1d and 2d (cf. [10] in general and [15, Remark 2.1] for precise references). We assume that the solution admits the following additional regularity, which is where we note that any solution with such increased regularity must be unique (cf. [15,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 (1) preserve the mass, i.e.
Here, w denotes the complex conjugate of w.
Throughout the paper we will use the notation A B, to abbreviate A ≤ CB, where C is a constant that only depends on u, T , d, D, V and γ, but not on the discretization.

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 ∞ (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 τ = T /N the time step size. By t n we shall mean t n = nτ . The exact solution at time t n shall be denoted by u n := u(t n , ·). We also introduce a short hand notation for discrete time derivatives which is D τ u n := (u n+1 − u n )/τ and analogously D τ u n τ := (u n+1 τ − u n τ )/τ .

Method formulation and main result
With the notation above, the semi-discrete Crank-Nicolson approximation u n+1 τ ∈ H 1 0 (D) to u n+1 is given recursively as the solution (in the sense of distributions) to the equation where u The initial value is selected as u 0 τ = u 0 . It is easily seen that the discretization conserves both mass and energy, i.e.
The scheme (3) and the a priori estimate for the L 2 -error where u is the (unqiue) exact solution with the regularity property (2).
Our main result on optimal error estimates in the L ∞ (H 1 ) reads as follows. The theorem is proved in Section 3.2 below.

Proof of Theorem 3.2
In this section we well prove Theorem 3.2. Let us introduce some notation that is used throughout the proofs. We recall D τ e n = (e n+1 − e n )/τ . 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 ∂ t u n := ∂ t u(t n , ·).
Proof. It is easily verified that exact solution fulfills By the regularity assumptions we can apply Taylor expansion arguments to T n to see: Subtracting (7) from (3) we find that e n = u n − u n τ satisfies: iD τ e n + ∆e n+1/2 − V e n+1/2 − e n γ = T n where e n γ denotes the error coming from the nonlinear term, defined by Recalling the definition of Γ we have: The expression for e n γ is thus simplified to where ξ n is a function taking values between |u n | 2 and |u n+1 | 2 and ξ n τ takes values between |u n τ | 2 and |u n+1 τ | 2 .
The differential equation in Lemma 3.3 is now used to derive a recurrence formula for the H 1 -norm of the error. Multiplying (5) by D τ e n , integrating and taking the real part yields:
In order to use Grönwall's inequality we need to bound e n γ and ∇e n γ in terms of e n , ∇e n and terms of O(τ 2 ). These bounds are formulated in the two following lemmas. Proof. We recall that e n γ = γ(ξ n )u n+1/2 − γ(ξ n τ )u n+1/2 τ and that γ(ξ n τ ) is defined by Since γ is twice differentiable we have: and hence R n τ τ 2 . The same argument, with the small difference that R n is bounded in virtue of the assumptions, can be applied to find Using this and the standard mean value theorem we find that Since |u n | 2 − |u n τ | 2 = (|u n | − |u n τ )(|u n | + |u n τ |)| 2 e n it now follows that e n γ τ 2 .
Lemma 3.5. Given Theorem 3.1, the gradient of the error coming from the nonlinear term is bounded as ∇e n γ ∇e n+1 + ∇e n + τ 2 .
With Lemma 3.4 and 3.5 we now have the following bound on term I. Lemma 3.6. For term I which is given by (10), we have the estimate We can now proceed to bound term II. Here we explicate the Taylor term using (6) to see
We start with estimating IIa, IIc and IId, which can be bounded in a similar way.
Collecting the estimates we have the following estimate for term II.
Lemma 3.7. For term II = −Re( T n , D τ e n ) it holds the estimate We are now ready to finish the proof of the first main result.
Proof of Theorem 3.2. We pick off where we left (9) and find by using Lemma 3.6 and 3.7: ∇e n+1 2 − ∇e n 2 2τ ≤ C ∇e n+1 2 + ∇e n 2 +τ 4 +D τ [ ∇(u n+1/2 −u(t n+1/2 )), ∇e n ]. (20) Summing this up and using e 0 = 0 gives and therefore Young's inequality with ǫ > 0 is used on the last term: Finally we arrive at ∇e n+1 2 ≤ C n k=0 τ ∇e k 2 + Cτ 4 + τ 4 ǫ + ǫ ∇e n+1 2 and for e.g. ǫ = 1/2 we can absorb ǫ ∇e n+1 2 in the left hand side and conclude Grönwall's inequality now yields: 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 ⊂ H 1 0 (D) denote the space of P1 Lagrange finite elements on a quasi-uniform simplicial mesh on D with mesh size h. In this setting we have by standard finite element theory (cf. [9]) the following inequality for any u ∈ H 1 0 (D) ∩ H 2 (D): Here P h : H 1 0 (D) → S h denotes (for example) the Ritz-projection into the finite element space. For a given discrete initial value u 0 h ∈ S h the CN-FEM approximation u n τ,h ∈ S h to u n is given by the fully discrete equation for all v ∈ S h . The initial value is selected as u 0 τ,h = P h u 0 . As in the semi-discrete case, the discretization conserves both mass and energy, i.e.  With this we are ready to state our final theorem. Proof. First, we recall the inverse estimate on quasi-uniform meshes (cf. [9]), i.e. ∇v h ≤ Ch −1 v h for all v h ∈ S h , which implies ∇(P h (u n τ ) − u n τ,h ) ≤ Ch −1 P h (u n τ ) − u n τ,h .
With this, the H 1 convergence result (22) 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.
Here we have made use of the inequality (23), the uniform H 2 -regularity of u n τ , i.e. u n τ H 2 ≤ C(u) (cf. (4)) and the optimal L 2 -estimates. In virtue of Theorem 3.2 we may thus conclude: Detailed numerical studies that confirm the optimal convergence rates stated in Theorem 4.1 and Theorem 4.2 are presented in [15,16].