Numerical comparison of mass-conservative schemes for the Gross-Pitaevskii equation

In this paper we present a numerical comparison of various mass-conservative discretizations for the time-dependent Gross-Pitaevskii equation. We have three main objectives. First, we want to clarify how purely mass-conservative methods perform compared to methods that are additionally energy-conservative or symplectic. Second, we shall compare the accuracy of energy-conservative and symplectic methods among each other. Third, we will investigate if a linearized energy-conserving method suffers from a loss of accuracy compared to an approach which requires to solve a full nonlinear problem in each time-step. In order to obtain a representative comparison, our numerical experiments cover different physically relevant test cases, such as traveling solitons, stationary multi-solitons, Bose-Einstein condensates in an optical lattice and vortex pattern in a rapidly rotating superfluid. We shall also consider a computationally severe test case involving a pseudo Mott insulator. Our space-discretization is based on finite elements throughout the paper. We will also give special attention to possible coupling conditions between time-step sizes and mesh sizes. The main observation of this paper is that mass conservation alone will not lead to a competitive method in complex settings. Furthermore, energy-conserving and symplectic methods are both reliable and accurate, yet, the energy-conservative schemes achieve a visibly higher accuracy in our test cases. Finally, the scheme that performs best throughout our experiments is an energy-conserving relaxation scheme with linear time-stepping proposed by C. Besse (SINUM,42(3):934--952,2004).


Introduction
This article deals with numerical methods for the nonlinear Schrödinger equation (NLSE) with cubic nonlinearity, which reads i∂ t u = − ∆u + V u + β|u| 2 u. ( The solution u(x, t) is a complexed valued wave function, β ∈ R is a given constant typically describing particle interactions, > 0 is a parameter related to the physical regime and V (x) : R d ⊃ Ω → R represents a potential. The problem is completed by an initial value u(x, 0) = u 0 (x). The equation has applications in several fields such as optics [2,20], fluid dynamics, in particular deep water waves [32,33], and quantum physics. Notably in quantum mechanics, it describes the dynamics of a Bose-Einstein condensate [16]. The equation is often referred to as the Gross-Pitaevskii equation.
When selecting a suitable discretization for problem (1), one needs to be aware that the equation possesses several time invariant characteristics. Two of the most important are the mass M and the energy E, which are preserved for all times. Here we stress that the selection of a numerical method that preserves mass and energy on the discrete scale is subject to the choice of a suitable time-discretization, but it is typically not affected by the choice of the spatial discretization. Another important characteristic of the problem is the Hamiltonian structure of the PDE. More precisely, by introducing q := Re(u) and p := Im(u) we observe that the NLSE can be interpreted as an infinite-dimensional Hamiltonian system of the form where H(p, q) = 1 2 E[u] is the Hamiltonian defined through the energy functional. Discretizing the problem in space while staying continuous with respect to the time variable leads to a finite dimensional Hamiltonian system. For instance, if p i and q i represent finite difference values then the Hamiltonian system is of the form ∂ t q i = ∂H h /∂p i and ∂ t p i = −∂H h /∂q i , where H h represents a corresponding discrete Hamiltonian. A second order finite difference method in one dimension would take the form The following symplectic form is then time invariant: It is well-known that symplectic time integrators only allow for small oscillations in the discrete energy, which is why they can be considered almost energy conservative. A central focus when solving the NLSE numerically has been to mimic time invariants, where the above mentioned invariants of mass, energy and symplecticity are considered to be the most important ones. An early comparative study was made by Sanz-Serna and Verwer [28] who compared closely related temporal discretizations based on finite differences in space. It was found that non mass-conservative methods were prone to nonlinear blow-up and performed poorly. Moreover it was found that, among their test cases, the best-performing method was symplectic and mass-conservative. A slight modification of this method leads to an energy-conservative Crank-Nicolson method as shown by Sanz-Serna [26]. However, the conservation of energy comes at the cost of losing symplecticity. In fact, for β = 0 and d > 1 it is impossible to discretely conserve all of the three aforementioned time invariants simultaneously, since energy conservation and symplecticity are mutually exclusive if the scheme is not exact (cf. [35]). Other numerical comparisons, most often of widely different approaches, have been done notably by Chang et al. [14] and Bao et al. [9]. A recent overview of the most common methods was done by Antoine et al. [7].
In the present paper we complement previous comparisons by answering which of the time invariants (3), (4), in addition to mass (2), is most important to conserve and how efficiently this can be done. A further aim is to find competitive methods to use in a low regularity regime. For that reason we will fix the space-discretization as being of finite element type (FEM).
The low regularity regime is motivated by the non-smooth potentials that arise in the context of disorder in quantum systems where one, for instance, encounters random potentials, see [15] and [25]. Another context is the investigation of quantum phase transitions from superfluid BEC to bosonic Mott insulators.
We stress that alternatively to a finite element discretization in space, one could also use a spectral or a finite difference discretization. In particular in high regularity regimes with smooth potentials and regular domains, the exponential convergence of a spectral approach leads to methods that are computationally very efficient. On the downside, they typically perform poorly if the regularity of the exact solution drops. For an overview of spectral methods for the Gross-Pitaevskii equation we refer to Bao et al. [11]. As for finite difference schemes, they are typically less effected by regularity issues. However, their accuracy still reduces observably. An early paper by Griffiths [19] found that finite element methods gave more favorable results than finite differences. A final argument for using FEM in connection with rough potentials is that it can be naturally combined with adaptive mesh refinement strategies, which may be key to efficient methods in low regularity regimes. For instance, rough disorder potentials are closely related to the phenomenon of Anderson localization [5,6,8] (exponential localization of eigenmodes), where often only locally refined meshes are required. For these reasons we restrict our attention to methods compatible with a finite element formulation. Concerning the application of nonlinear Schrödinger equations, we deem exponential time integrators to have low compatibility with FEM discretizations as this requires a mass lumping approach to be efficiently implemented. Since lumping can introduce a severe perturbation of the properties that are to be conserved, we exclude such methods in our survey.
All methods that we consider are mass-conservative and implicit. The necessity of massconserving schemes to obtain reliable approximations was demonstrated in the extensive numerical studies in [28]. Here, we also refer to the numerical experiments in [21] where the Backward Euler method shows a devastating performance due to its lack of mass conservation. If the nonlinearity is treated implicitly, then a nonlinear system of equations must be solved in each time-step. This is the case for the symplectic one-stage Gauss-Legendre Runge-Kutta method that Sans-Serna and Verwer [28] found to perform best and also for the energy conservative Crank-Nicolson method.
Albeit implicit methods, convergence rates were, until recently, always obtained under a (counterintuitive) time-step condition in terms of the mesh size. Recent developments include analysis of linear mass-conservative methods by Wang [31] (2013), Besse [12] (2004), and Zouraris [36] (2001) whose methods we take into account in our comparison. J. Wang proposes an Adams-Bashforth like linearization of the Crank-Nicolson method and was first to prove optimal convergence rates in L 2 -norm without any condition on the time-step size. The method proposed by Besse conserves, in addition to mass, energy. The Two-Step method analyzed by Zouraris was proved to converge under moderate time-step restrictions. The Relaxation method figures in the overview article by Antoine et al. [7]. The authors present errors for a numerical experiment in which it, and the Crank-Nicolson method perform best in terms of temporal accuracy. Other than this, the linear methods have not, due to their recentness, figured in comparative studies and form an interesting complement to the well established nonlinear methods.
Using a finite element based spatial discretization, we test these five mass-conservative methods on problems of increasing difficulty in one and two dimensions. The test cases include solitons with known analytical solution such as the single traveling soliton and a newly derived stationary soliton [4]. Problems of physical relevance involving optical lattices and angular momentum are also considered. We do not resort to the usage of an artificial source function often found in the literature, as it would distort the physical behavior and break the conservation of mass, energy and symplecticity.
The outline of the paper is as follows: in Section 2 the aforementioned five methods are presented along with their properties and a brief evaluation with respect to both theory and implementation. A summarizing table concludes the section. In Section 3 and Section 4 we present numerical experiments in one and two dimensions. Finally, conclusions are presented in Section 5.

General problem formulation and methods
In this section we will state the five numerical methods that are at the core of our numerical study. We present the methods in a FEM formulation in space along with their properties and computational complexity. For a proper mathematical description of the setup, we let Ω ⊂ R d denote a bounded computational domain and let H 1 0 (Ω) = H 1 0 (Ω, C) denote the Sobolev space of weakly differentiable, square-integrable and compactly supported functions over the complex field. By z we denote complex conjugation of a complex number z and by u, v = Ω uv the L 2 -inner product. With this, we consider the following initial value problem in weak form. Given a constant > 0, a real valued bounded potential V ∈ L ∞ (Ω) and an initial value , H −1 (Ω)) such that u(·, 0) = u 0 and such that: for all v ∈ H 1 0 (Ω) and a.e. t ∈ (0, T ]. Remark 2.1. If β ≥ 0 and V ≥ 0, then problem (5) has at least one solution. This solution is also unique if d = 1, 2. For d = 3, a corresponding uniqueness result is still open. For a proof of these results, as well as a more general existence theory in the focusing regime β < 0 we refer to the excellent textbook by Cazenave [13].
In the following, we let S h denote the space of P1 Lagrange finite elements on a quasiuniform simplicial mesh on Ω with mesh size h. We use standard notation for the timediscretization, i.e. for final computational time T = N τ and n = 0, · · · , N we have t n = nτ and u n (x) = u(x, nτ ). To facilitate reading we also define u n+1/2 = (u n+1 + u n )/2.
As for complexity, it is noted that nonlinear methods suffer from an additional drawback. Namely, if the nonlinear system of equations is solved through a Newton step, then it cannot be done in the complex field since there is a non-holomorphic mapping, in this case z → |z|. To make this aspect explicit, let The solution must therefore be split into real and imaginary part. This leads to a system of equations of size (2m) 2 where m is the number of Lagrange elements. Linear methods lead to systems of size m 2 (over C) with 2m 2 degrees of freedom.
The following methods are selected to complement the earlier numerical study by Sanz-Serna and Verwer [28]. Among the methods investigated in the aforementioned comparison paper, we only include the symplectic Gauss-Legendre Runge-Kutta method since it performed best in their study.

Implicit Midpoint Method (IM-FEM)
The first scheme in our comparative study belongs to a general class of Gauss-Legendre Runge-Kutta methods, where we shall only consider the lowest order case (implicit midpoint method with finite element space-discretization -IM-FEM). For realizations of higher order we refer to [24,27]. It has been shown by Sanz-Serna [27] that all Gauss-Legendre Runge-Kutta methods are symplectic, which is also why the following one stage realization is symplectic.
For a given discrete initial value u 0 h ∈ S h the IM-FEM approximation u n h ∈ S h to u n for n = 1, · · · , N is given by the time-discretization for all v ∈ S h . Besides being symplectic, the method is also mass-conservative (i.e. u n h L 2 (Ω) = u 0 h L 2 (Ω) for all n ≥ 0) as is easily seen by testing in (6) with v = u n+1 h + u n h and taking the imaginary part of the equation. We stress again that this method showed the best performance in the numerical comparison by Sanz-Serna and Verwer [28]. In order to solve the nonlinear algebraic equation given by (6), an inner loop of Newton iterations can be applied.
A first proof of convergence in the fully discrete finite element setting was given in 1984 [30] for d = 1. Later, for arbitrary space dimension, optimal L 2 -convergence rates of order O(τ 2 + h 2 ) were proven by Akrivis et al. for V = 0 and under the time-step condition τ = O(h d/4 ). Optimal L ∞ -and H 1 -error estimates for u N − u N h were obtained by Tourigny [29], who recovers again the constraint τ = O(h d/4 ) to ensure convergence. In [21] it was shown that for optimal L 2 -and H 1 -error estimates, the constraint can be relaxed to τ = O(| ln h| −(1+ε)/4 ) for d = 2 and to τ = O(h (1+ε)/4 ) for d = 3. Here, ε > 0 can be arbitrary close to zero. We note that the analysis in [21] also amounts for potentials V ≥ 0, provided that they are smooth enough.

Crank-Nicolson Method (CN-FEM)
A slight modification of the previous method yields the energy-conservative Crank-Nicolson method (CN-FEM). The method has been of widespread use in the physics community and reads: given u 0 h ∈ S h find u n h ∈ S h , n = 1, · · · , N , such that for all v ∈ S h : Note that the essential difference between the methods (6) and (7) is that the IM-FEM is based on an average of the wave-functions in the nonlinearity, whereas the CN-FEM is based on an average of the densities.

Remark 2.2.
A generalization of the IM-FEM to a larger class of nonlinearities is straightforward. To generalize the CN-FEM, antiderivatives of the nonlinearity are needed. For instance, a generalization of the Crank-Nicolson method for nonlinearities of the form γ(|u| 2 )u for some smooth and integrable γ reads as follows in the fully discrete case: The method is mass and energy conservative. A numerical analysis of this method for V = 0 was first done by Sanz-Serna in one space dimension [26], in two and three dimensions by Akrivis et al. [3]. Akrivis et al. obtained optimal L 2 -convergence rates, i.e.
Recently it was shown that such coupling conditions are in fact not required. More precisely, in [22] optimal convergence in L 2 is proved without any time-step conditions for d = 1, 2, 3, for sufficiently smooth V and a general class of nonlinearities including the cubic case γ : z → |z| 2 . Furthermore for general potentials V ∈ L ∞ (Ω) the authors of [22] prove suboptimal convergence rates of order O(h (d+α)/2 + τ ) for some α > 0 under the mild condition h 4−d−α ≤ Cτ 2 . We note that a proof of convergence in H 1 is still open in the literature.
Practically, the CN-FEM in (7) requires the solving of a nonlinear system of equations in each time-step, a suitable Newton step was proposed and analyzed by Akrivis et al.

Relaxation Method (RE-FEM)
A relaxation method was put forward by Besse [12] and treats the nonlinearity explicitly. A remarkable property of this method is that it conserve mass and energy simultaneously, while being linear in each time-step. Introducing ρ τ = |u τ | 2 , the relaxation method in semi-discrete form reads as follows for n ≥ 0. Find u n+1 where the density approximation ρ n+1/2 τ is recursively defined for n ≥ 0 by Here, u n τ denotes the obtained numerical approximation for u(·, t n ) and u 0 τ = u 0 is the original initial value. For the above time-discrete system with V = 0, convergence was proved under the assumption of sufficiently high regularity [12]. More specifically, if then the sequence of time-discrete solutions (u τ , ρ τ ) to (8) to the exact solution (u, |u| 2 ) as τ → 0.
In the fully discrete case, the Relaxation Finite Element Method (RE-FEM) reads as follows in variational form. Given a suitable approximation u 0 h ∈ S h to the exact initial value u 0 , find u n h ∈ S h , n = 1, · · · , N , such that for all v ∈ S h : The following energy like quantity is exactly conserved: Note that if u n h is as assumed piecewise linear, then for energy conservation it is imperative that the discrete density, ρ n h , be accurately represented as a piecewise polynomial of degree two. This needs to be considered when implementing the method.
We are not aware of any analytical results concerning the convergence of the fully discrete method.

Linearized Crank-Nicolson Method (LCN-FEM)
The following method, first proposed and analyzed by Wang [31], is based on replacing the nonlinear term by an Adams-Bashforth linearization. For that we define the short hand notation and for n = 0, we let u 0 h ∈ S h denote the solution to the intermediate first step With this, the Linearized Crank-Nicolson Finite Element Method (LCN-FEM) reads: find Applying this particular Adams-Bashforth linearization can be motived by elliptic regularity theory: apart from the space-discretization, the problem that needs to be solved for each time-step is a linear elliptic problem that admits H 2 -regularity. Exploiting this observation allows to straightforwardly derive uniform L ∞ -bounds for the arising semi-discrete approximations. This a priori control over the growth of the numerical approximations is the key to a convergence theory that does not require time-step constraints. At the cost of this elegant approach, energy-conservation and symplecticity are lost. We stress that Wang [31] was the first to prove optimal L 2 -convergence rates O(τ 2 + h 2 ) without coupling conditions for the mesh size and the time-step size. The analysis was carried out under the assumption that V = 0. The paper also provides a proof of optimale order convergence in H 1 , i.e. convergence with the rate O(τ 2 + h).

Linearized Two-step Method (Two-Step FEM)
In [36], Zouraris proposed the a new type of linearized method to which we will refer as Two- Step FEM. In the initial double step, the method requires to first solve for an intermediate After that, the first full step is taken by solving for With this, the iterations of the two-step method for n = 1, · · · , N − 1 read: In [36], the method was analyzed for the case of adaptive meshes, though we restrict our discussion of the results to the quasi-uniform case as before. For V = 0, optimal convergence rates in the L 2 -and H 1 -norm of order O(τ 2 +h 2 ) and of O(τ 2 +h 2 ) respectively are presented under certain constraints for the time-step size. For d = 1 there is no constraint at all and for higher dimensions the time-step size needs to be bounded in terms of the mesh size by fulfilling the following relations We note that no such coupling constraints became visible in our numerical experiments.

Summarizing table
The following table summarizes the essential features of the considered methods, i.e. mass and energy conservation, symplecticity and linear/nonlinear time-stepping. We recall that all methods have the same optimal convergence rates in space and time, provided that the solution is sufficiently regular. As for the detailed convergence results including the required assumptions we refer to the discussion in the corresponding subsections above. We stress again that the nonlinear CN-FEM is the only method for which convergence for rough potentials is proved. We do not include this fact in the table below. However, we shall indicate in the table if optimal L ∞ (L 2 )-or L ∞ (H 1 )-estimates are available for a certain method, provided sufficient regularity is available, and if the proofs for these convergence rates required a formal coupling condition between mesh size and time-step size. By conditional we mean that a coupling condition was required and by unconditional that it was not required. The precise conditions can be found in the previous subsections.
and an optimal convergence rate in H 1 is of order O(τ 2 + h). As everywhere in the paper, this refers to the case of piecewise linear FEM (simplicial Langrange elements) for the space-discretization.
Finally, we stress that the computational complexity of the linearized methods RE-FEM, LCN-FEM and Two-Step FEM is basically the same (for a fixed space and time resolution). In comparison, the computational complexity of the nonlinear methods IM-FEM and CN-FEM is roughly ten times higher in one dimension. An exemplary comparison of CPU times is given in Table 4 in Section 3.2 . As we will see later, the higher CPU times can be justified for some test cases by a significantly higher accuracy.

Single Soliton
In the first numerical experiment we consider a single traveling soliton u given by It is characterized as the solution to the following NLSE in full space This solution originates from an early paper by Shabat and Zakharov [34]. For any fixed time t, the solution given by (12)  We run all five methods stated in Section 2 for various mesh sizes and time-step sizes. The discrete initial value is the nodal interpolation of u(x, 0) in the finite element space. As we will see, all methods perform very well for the test case, where the Two-Step FEM performs best in terms of error size and CPU times.
To support this claim, a first set of results is depicted in Fig. 1a. Here, the L 2 -error for the density ρ = |u| 2 at time T = 10 is plotted versus the time-step size τ . The nonlinear methods start with substantially smaller errors than the linear methods. Overall the nonlinear methods outperform the linear methods for coarse time steps, where the energyconservative CN-FEM has only slightly smaller errors than the symplectic IM-FEM. For smaller time steps, the linear methods catch up in terms of performance due to their shorter CPU times (making them 10 times faster than the nonlinear methods). We observe that the LCN-FEM initially converges cubically in the L 2 norm of the density, it does however start with the highest errors. At τ = 2 −7 it is on par with the Two-Step FEM and the cubic convergence flattens out to the expected quadratic convergence. Asymptotically, the Two-Step FEM and the LCN-FEM are on par with roughly four times larger errors than the nonlinear methods, the RE-FEM in turn produces about four times larger errors than the Two-Step FEM.
In Fig. 1b the error in the H 1 -norm is plotted versus the time-step size. Again, the Two-Step FEM fares best. Next in accuracy comes the energy conservative CN-FEM and only slightly less accurate is the energy conservative RE-FEM which is on is on par with the symplectic IM-FEM. The LCN-FEM consistently has about 6 times higher errors than the CN-FEM and roughly 32 times that of the Two-Step FEM.
It is noted in Fig. 1b that for approximately τ = 2 −9 the error in H 1 -norm from the space-discretization with h = 2 −10 seems to become dominant in all methods, resulting in polluted convergence rates in τ .
An interesting observation is made for the LCN-FEM, which achieves remarkably low L 2 -errors if the mesh size and the time-step size are coupled to the ratio h/τ = 2. In this case, the method still shows a quadratic convergence, but the errors are significantly smaller. This is illustrated in Fig. 2a. The ratio does not depend on final time and is not observed for the other methods. As seen in Fig. 2b this feature virtually disappears in the error plot of the density, implying that for this ratio the error in phase is minimal. That is, a small change in the spatial discretization may result in considerably smaller L 2 -errors, but not in smaller density errors. Since the density barely changes, the phase must account for the smaller errors.

Two Stationary Solitons
A more difficult test case related to signal propagation in optical fibers and involving two interacting stationary solitons was constructed in [4]. In this example we are looking for the solution u of the following NLSE with a focusing nonlinearity and with initial value As derived in [4], the exact solution is given by The function u(x, t) is depicted in Fig. 3 and consists of two cocentric interacting solitons. With this it is less smooth than the single traveling soliton considered before. However, as in the previous example, we have an exponential decay with |u| → 0 for |x| → ∞.  Due to the exponential decay, we restrict our computations to a finite computational domain of size [−20, 20] × (0, 2] and prescribe a homogeneous Dirichlet boundary condition on both end of the spatial interval. Again, we make a numerical comparison of all five methods, where the initial value is the nodal interpolation of u(x, 0) in S h .
In Table 2 we present convergence of the density in L 2 -norm for the spatial discretization with mesh size h = 40/51200. For the same spatial discretization, H 1 -norm convergence results are presented in Table 3. From these tables we can make an important observation in terms of the convergence rate of the numerical methods. None of the methods converges right away, but the regime of an asymptotic convergence is entered at different time-steps sizes. For instance, the symplectic IM-FEM starts to show the quadratic convergence (in time) from around τ = 2 −8 , the energy-conserving CN-FEM and the Two-Step FEM from around τ = 2 −7 , the energy-conserving RE-FEM from around τ = 2 −6 (and even τ = 2 −5 for the H 1 -error) and the LCN-FEM from around τ = 2 −8 /2 −9 . Due to this early start into the asymptotics, the RE-FEM scheme performs best and achieves the smallest errors. The CN-FEM is close behind, however has a higher computational complexity due to the nonlinear time-stepping.
For very large time-steps, e.g. τ = 1/8, the nonlinear methods (IM-FEM and CN-FEM) become extremely ill-conditioned and all accuracy is lost. This is related to the focusing nonlinearity (i.e. β < 0) and we could not observe similar effects in the defocusing case (i.e. β ≥ 0).   The LCN-FEM experiences large increases in energy for larger time-steps, leading to a large error in the H 1 -norm. As in the previous example it fares remarkably well for some ratio of h and τ (Fig. 4). However, notably and as opposed to the single soliton case this ratio seems to be dependent on the final time as we find for T = 2 that optimally h/τ ≈ 8.2, changing T = 10 we find that optimally h/τ ≈ 10.2. Furthermore, even with an optimal h/τ ratio the method does not beat the energy conservative methods.
The RE-FEM performs the best and the energy conservative methods outperform by far the symplectic IM-FEM.

Condensate in an optical lattice
We conclude the experiments in 1D with a test case that describes a Bose-Einstein condensate in an optical lattice. This test case demonstrates how the LCN-FEM and the Two-Step FEM can develop an instability resulting from a blow-up of the energy. As the energy remains (almost) conserved for IM-FEM, CN-FEM and RE-FEM, these methods do not suffer from such an instability. For the space interval Ω := (−12, 12) we seek the solution u to with Dirichlet boundary condition u(−12, t) = u(12, t) = 0 for t ≥ 0 and initial condition u(·, 0) = u 0 in Ω. The final time is varied. The problem involves a defocusing nonlinearity with β = 1000. The potential describes the combination of a harmonic confining potential together with an optical lattice and is given by where the trapping frequency is γ x = 1 2 .
The initial state u 0 (x) is selected as the ground state associated with the trapping frequency γ x = 1. This ground state is accurately computed using the inverse iteration method for eigenvalue problems with eigenvector nonlinearities presented in [23]. Decreasing the strength of the harmonic potential for t > 0 (as described above) induces a small periodic motion of the condensate which, in our case, is concentrated around the lattice points [−8, −4, 0, 4, 8].
The density of the condensate remains almost constant. However, the setup is very sensitive to small perturbations of the energy and causes instabilities for the two non-symplectic / non-energy-conservative methods. This is stressed by Fig. 5a, where the energy is plotted versus time for the Two-Step FEM with spatial resolution h = 24/800 and time-step sizes τ = 2 −8 , 2 −9 , . . . , 2 −20 . For the LCN-FEM, a similar plot but for a larger final time is presented in Fig. 5b. The spatial discretization is as before h = 24/800 and τ ranges from 2 −9 to 2 −14 . Energy blow-up occurs at some time for both methods of which the LCN-FEM blows up more violently. However, the most striking feature is that refining the time-step barely improves the stability of the Two-Step FEM, i.e. the blow-up is hardly delayed. On the contrary, refining the time-step greatly delays the blow-up of the LCN-FEM. Furthermore, we note that for the LCN-FEM this blow-up is coupled with the mesh size. This is also illustrated in Fig. 5b, as changing mesh size from h = 24/800 (black dashed line) to h = 24/1600 (blue dashed line) while keeping τ = 2 −14 changes the method from being stable at least till t = 20 to energy blow-up occurring at t = 1.1. Some further computations, not included for sake of brevity, show that this coupling behaves rather erratically. Finding the most stable ratio may therefore prove difficult.
Lastly we present errors and CPU-times for the five methods in Table 5. It is clear that the linear methods, with the exception of the Two-Step FEM, achieve lower errors than the nonlinear methods for the same computational times. It is also noted that the RE-FEM and CN-FEM produce almost identical solution, which is explained by the fact that the change in density is small and that for constant density the CN-FEM and the RE-FEM are equivalent. As mentioned, the Two-Step FEM does not perform well compared to the other methods due to its critical instability. Therefore and for the sake of brevity, it will henceforth not be considered in the remaining experiments.

Optical Lattice
In the next experiment we will consider an optical lattice in 2D. This example serves to show that despite being both mass-conservative and symplectic, the IM-FEM has an instability that manifests itself in the density. This instability is even more pronounced in the LCN-FEM, but fully absent in the energy-conservative methods. In the corresponding test problem, we seek u(x, t) with Here, Ω = (−6, 6) 2 is the computational domain and for the maximum time we selected T = 1. The potential V = V o + V ha + V c , consists of three parts: an optical lattice V o , a harmonic confining potential V ha and a discontinuous confining potential frame V c . We have For the dynamics in the time-dependent problem we set the trapping frequencies to γ x = 4 and γ y = 8. The initial value u 0 is the ground state with γ x = γ y = 1, i.e. it solves the eigenvalue problem with ground state eigenvalue (chemical potential) λ 0 . The initial state is accurately computed using the inverse iteration method for eigenvalue problems with eigenvector nonlinearities presented in [23]. The solution to (13) approximately consists of 25 localized oscillators. Again the mass-conservative LCN-FEM is prone to blow up in energy as indicated in Fig. 9b. This results in meaningless solution plots as can be seen in Fig. 6d and Fig. 8d. Fig. 9a shows that there is a pre-asymptotic regime and furthermore that for h = 0.03, the LCN-FEM requires 4 times higher resolution as the other methods to enter this asymptotic regime, namely τ = 2 −13 . This is clearly coupled with h as we find for h = 0.015 that the energy blows up even for τ = 2 −14 .
Perhaps more notably, the symplectic IM-FEM produces a deteriorated solution for coarser time-steps as the error in density blows up, cf. Fig. 6a and Fig. 8a. Since the symplectic method does not allow a blow-up of the energy, the deterioration seems to be already triggered by small fluctuations in the energy. This guess is supported by the observation that the energy-conservative methods perform equally well and never produce completely erroneous density plots. The energy evolution of the symplectic IM-FEM is shown in Fig.  7, where the energy oscillations become clearly visible.

Rotating Bose Einstein Condensate
The next experiment shows how the low errors of the LCN-FEM can be misleading as the method may still produce polluted plots. In this experiment we consider a rotating Bose-Einstein condensate (BEC) which is typically used to investigate superfluidity on an observable scale (cf. [1,17]). Here the superfluid character of the condensate is expressed through density singularities, so called-vortices, and it is of interest to study the dynamical behavior of such vortex patterns. In an appropriate mathematical model, we first compute the initial value as the ground state of a BEC under angular momentum rotation, i.e. we seek an L 2 -normalized eigenfunction u 0 ∈ H 1 0 (Ω) and corresponding minimal eigenvalue λ > 0 such that Here, denotes the axis-oriented angular momentum operator and ω is the angular velocity. In our experiment we selected Ω = (−6, 6) 2 , V 0 (x, y) = 1 2 x 2 + 1 2 y 2 , β = 100 and ω = 0.8. The eigenvalue problem was solved numerically using the Discrete Normalized Gradient Flow method proposed in [10]. The computational mesh is the same as the mesh used for solving the time-dependent problem afterwards. Here we simulate the situation that the stirring potential (represented by ωL) is switched off after the condensate reaches the ground state. Due to the initial rotation the condensate still has an angular momentum and we can study the dynamical change of the vortex pattern. Modifying the trapping potential so that it models a slightly anisotropic trap, we consider the problem to find u such that on Ω with anisotropic potential V (x, y) = 0.45x 2 + 0.55y 2 and the initial condition u(x, 0) = u 0 , where u 0 is the ground state from above. The computations are carried out for maximum time T = 10. As a reference solution in all computations we use an approximation obtained with the RE-FEM with a time-step size τ = 2 · 10 −4 and mesh size h = 12 · 2 −6 . The mesh size will be constant for all computations so that we can focus on the error caused by the time-discretization alone. This experiment further corroborates the general conclusions that in our experiments the energy conservative methods are slightly more accurate than the symplectic method and that the linearized Crank-Nicolson method is prone to cause unphysical effects. In particular, the LCN-FEM can produce erroneous density plots while still maintaining low errors. A striking example of this is Fig. 11 together with the Table 6. The figure shows the density plots of the LCN-FEM, CN-FEM and the RE-FEM for a fine spatial discretization and τ = T /256. The LCN-FEM solution is polluted but yet the table shows that it has the lowest error in both L 2and H 1 -norm as well as the lowest error in terms of L 2 -norm of the density. The reason for this observation is that LCN-FEM captures the rotational phase significantly better than the other methods. As for a comparison between IM-FEM, CN-FEM and RE-FEM we observe that all methods have a similar accuracy, where the errors become smallest for the energyconserving CN-FEM. However, noting the considerably lower computational complexity of the energy-conserving RE-FEM, we can clearly identify it as the best-performing method for the test case.
It is also worth to mention that all methods suffer from a rotational phase shift compared to the exact solution. This is shown exemplarily for the IM-FEM and RE-FEM in Fig. 10. The isolines between the two numerical approximations are hardly distinguishable, where the RE-FEM is only slightly better at respecting the isolines of the exact solution than the IM-FEM. One such slight difference can be seen in Fig. 10 for the rightmost vortex, where the IM-FEM produces a closed loop whereas the RE-FEM does not.

Mott insulator to superfluid
In our final experiment we use the RE-FEM, which was best-performing in our experiments, to present an interesting 2D simplification of the celebrated experiment by Greiner et al. [18]. Before sketching the setup of the physical experiment, we stress that the computational complexity of the nonlinear methods was too high to carry out the computations at the desired resolution. The following numerical experiment therefore differs from the previous experiments in the sense that we will not compare the various methods and we do not have a numerical reference solution. Instead, we evaluate the numerical results by comparing them  with the observations from the physical experiment. As we will see, the energy-conservative RE-FEM successfully captures the expected phase transition. In the experiment by Greiner et al. [18] an optical lattice is used to study quantum phase transitions of a gas of ultra cold bosons from the Mott insulator phase into the superfluid phase. In the Mott insulator phase, individual atoms are trapped in the sites of an optical lattice and the energy of particle interactions dominates over the kinetic energy. When the optical lattice is removed, the particles are free to move and the wave packages start to interfere with each other (phase transition). The aspect of physical interest is to study the subsequent interference pattern that is formed (or not formed) after the bosons are released from the lattice. As observed in [18] the pattern is closely related to the strength of the lattice potential. More precisely, as the strength of the optical lattice is increased, the interference pattern after release from the lattice goes from a regular high-contrast pattern to an incoherent background.
When the optical lattice potential is increased, the particle to particle interactions grow, at a certain threshold these interactions become too large for the gas to be described by a macroscopic wave function as in the Gross-Pitaevskii model. Instead, the (computationally very heavy) Bose-Hubbard model is needed for an accurate description of Mott insulators. In the following we shall assume that the energy of the particle interactions is high, but still below the threshold, so that we can use the Gross-Pitaevskii equation to compute an approximation of a Mott insulator ground state. In any case, after switching off the lattice, the Gross-Pitaevskii equation is the valid model for studying the dynamics of the ultra cold bosons and hence the formation of interference. A simulation faithful to the experiment by Greiner et al. in 3D could well require, of the order, 10 12 degrees of freedom for the spatial discretization. As this is beyond computability, we will present a simplified 2D case with interesting dynamics requiring reasonable computational power.
From the physical setup described in [18], we can extract a configuration and a parameter range in which we expect to observe the phase transition from Mott insulator to superfluid together with a clear and structured interference pattern. As an admissible configuration for our numerical experiment, we set = 1/2, particle interactions with β = 1000 and the following harmonic and optical lattice potentials: The initial state (as an approximation for the Mott insulator) is taken to be the ground state associated with both the harmonic and the optical potential, i.e. u 0 ∈ H 1 0 (Ω) solves (V h + V o + β|u 0 | 2 )u 0 = λ u 0 , where λ denotes the smallest eigenvalue. We solve the timedependent Gross-Pitaevskii equation (1) on [−20, 20] 2 with homogeneous Dirichlet boundary conditions keeping only the harmonic potential, thus V = 2(x 2 1 + x 2 2 ). We solve until final time T = 0.515 with h = 0.01, accounting for 16 · 10 6 spatial degrees of freedom, and the number of time-steps is 2048, i.e. τ = 0.515/2 11 .
As illustrated in Fig. 12a, the condensate starts in a regular pattern of localized and smooth clouds. After release, it quickly disintegrates into a more chaotic interference pattern with sharp spikes, illustrated in Fig. 12b. At t = 0.515 highly localized and smooth clouds reappear in a regular lattice, shown in Fig. 12d. The interference pattern shown in Fig. 12d has precisely the structure of the pattern observed in the physical experiment by Greiner et al. We stress that this is only possible if the interference of the individual wave packages is accurately captured by the numerical scheme over many time steps. As we can see, the energy-conserving RE-FEM is able to reproduce the desired behavior.

Conclusion
In this paper, we compared five different mass conservative time discretizations for the Gross-Pitaevskii equation. Two of the schemes (CN-FEM and RE-FEM) additionally preserve the energy, one scheme (IM-FEM) additionally preserves the symplectic structure and two schemes (LCN-FEM and Two-Step FEM) only preserve the mass. We observed that all schemes perform very well for simple experimental setups, where energy-conservation and symplecticity do not seem to play a crucial role to obtain good results. However, when moving towards more complicated setups with reduced regularity and additional potential terms, the conservation of energy becomes more and more important and the purely mass-conservative methods, LCN-FEM and Two-Step FEM, can fail to produce accurate approximations on reasonable scales. We also observed that in our experiments the energy-conserving methods gave visibly better approximations than the symplectic IM-FEM. None of our experiments gave any indications that convergence could depend on a possible coupling between mesh size and time-step size as often seen in the literature to prove convergence. Still we saw that an unlucky choice of the ratio h/τ can cause an energy blow-up for the LCN-FEM and the Two-Step FEM. This blow-up vanishes slowly for h, τ → 0, hence not contradicting convergence, however it can introduce severe practical constraints. Surprisingly, we also find that the symplectic IM-FEM can produce completely erroneous density plots, whereas the energy-conservative methods do not. In short, the energy-conservative methods appear to be more reliable in complex physical settings. The linear energy-conservative RE-FEM achieves, with the exception of the single soliton test-case, errors at least on par with the nonlinear CN-FEM, thereby making it by far the most computational efficient method.