In this paper we consider nonlinear Schrödinger equations (NLS) seeking a complex function such that
in a bounded domain , with a homogenous Dirichlet boundary condition on and a given initial value. Here, is a known real-valued potential and is a smooth (and possibly nonlinear) function that depends on the unknown density . Of particular interest are cubic nonlinearities of the form , for some . 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 for some , which appear in the context of nonlinear optical wave propagation in layered metallic structures [11, 17] or the propagation of light beams in plasmas .
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 for and , 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 . 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  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 finite element discretization in space, the first a priori error estimates for the arising method were obtained by Sanz-Serna 1984  for cubic nonlinearities. He considers the case and derives optimal -error estimates under the coupling constraint , where denotes the time step size of the Crank-Nicolson method and the mesh size of the finite element discretization in space. In 1991, Akrivis et al.  improved this result by showing optimal convergence rates in in dimension and under the relaxed coupling constraint . Finally, in 2017 , the -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 are still open in the literature.
One reason for this absence of -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 - and/or -error estimates. 3. Using inverse estimates in the finite element space to show that the truncated approximations are uniformly bounded in by a term of the form , with appropriate powers that depend on the considered space discretization, regularity and space dimension. 4. Concluding that if and are coupled in an appropriate way, then the truncated approximations are all uniformly bounded by a constant 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 -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 . However, thanks to the new techniques developed in  and the CN error analysis suggested in  in the context of -error estimates, the truncation step is no longer necessary and the desired -bounds can be derived with elliptic regularity theory. With this, it is now possible to obtain 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 -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 . Furthermore, we present and prove the new error estimate in . 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 -error estimates, which is the main result of this paper.
2 Notation and Assumptions
We start with introducing the analytical setting of this work. Throughout the paper we assume that (for ) is a convex bounded domain with polyhedral boundary. On , the Sobolev space of complex-valued, weakly differentiable functions with a zero trace on and -integratable partial derivatives is as usual denoted by . The potential is assumed to be real and nonnegative. Indirectly, we also assume that is sufficiently smooth so that it is compatible with the regularity assumptions for listed below (see  for a discussion on this aspect). The (possibly nonlinear) function
is assumed to be smooth, fulfills and its growth can be characterized with
where is a function with the following growth properties
Note that in  the admissible growth condition in requires , which is however a typo and should be, as above, (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 for and in and in . Other physically relevant nonlinearities that fulfill the conditions are saturated nonlinearities appearing in the modeling of optical wave propagation such as for .
For the initial value we assume that and, without loss of generality, that it has a normalized mass, i.e. . With this, the considered nonlinear Schrödinger equation (NLS) reads as follows. For a maximum time and an initial value , we seek
such that and
in the sense of distributions. Problem (1) admits at least one solution, that is even unique for repulsive cubic nonlinearities in and (cf.  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 hence always refers to this uniquely characterized solution.
It is well known that solutions to the NLS (1) preserve the mass, i.e.
and the energy, i.e.
For brevity, we shall denote the -norm of a function by . The -inner product is denoted by . Here, denotes the complex conjugate of .
Throughout the paper we will use the notation , to abbreviate , where is a constant that only depends on , , , , and , but not on the discretization.
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 -error estimates in the Hilbert space setting. For that, let denote the final time of computation, the number of time-steps, and the time step size. By we shall mean . The exact solution at time shall be denoted by . We also introduce a short hand notation for discrete time derivatives which is and analogously .
3.1 Method formulation and main result
With the notation above, the semi-discrete Crank–Nicolson approximation to is given recursively as the solution (in the sense of distributions) to the equation
where . The initial value is selected as . It is easily seen that the discretization conserves both mass and energy, i.e.
Our main result on optimal error estimates in the reads as follows.
Theorem 3.2 (Optimal -error estimates for the semi-discrete method).
Consider the setting of Theorem 3.1, then the -error converges with optimal order in , i.e.
The theorem is proved in Section 3.2 below.
3.2 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 . Furthermore, we let and . For time derivatives at fixed time , we also write .
We begin by establishing a differential equation for the time discrete error . This is stated in the following lemma.
The error fulfills the identity
and for some bounded functions with the properties that
for almost all .
It is easily verified that exact solution fulfills
By the regularity assumptions we can apply Taylor expansion arguments to to see:
where denotes the error coming from the nonlinear term, defined by
Recalling the definition of we have:
The expression for is thus simplified to
where is a function taking values between and and takes values between and . ∎
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 (5) by results in:
In order to use Grönwall’s inequality we need to bound and in terms of , and terms of . These bounds are formulated in the two following lemmas.
We recall that and that is defined by
Since is twice differentiable we have:
and hence . The same argument, with the small difference that is bounded in virtue of the assumptions, can be applied to find
Using this and the standard mean value theorem we find that
Since it now follows that ∎
Given Theorem 3.1, the gradient of the error coming from the nonlinear term is bounded as
Starting from the equality
we recall that
where . Using formula (14) with and we find:
where the remainder is bounded as . It is straightforward to check that . Taking the gradient of (3), multiplying by and integrating ( is interpreted weakly, is therefore integrated by parts), gives us:
Recalling it becomes clear that
This means that . Thus,
Using this and the mean value theorem we find:
As the exact same steps are valid for , with the difference that is bounded in virtue of the assumptions, the term becomes:
Since is three times differentiable it follows that:
where the remainder is bounded as . Using the mean value theorem we find:
Thus we are able to conclude:
For term I which is given by (3.2), 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 and IId, which can be bounded in a similar way.
Step 1, bounding IIa: