On optimal convergence rates for discrete minimizers of the Gross-Pitaevskii energy in LOD spaces

In this paper we revisit a two-level discretization based on the Localized Orthogonal Decomposition (LOD). It was originally proposed in [P.Henning, A.M{\aa}lqvist, D.Peterseim. SIAM J. Numer. Anal.52-4:1525-1550, 2014] to compute ground states of Bose-Einstein condensates by finding discrete minimizers of the Gross-Pitaevskii energy functional. The established convergence rates for the method appeared however suboptimal compared to numerical observations and a proof of optimal rates in this setting remained open. In this paper we shall close this gap by proving optimal order error estimates for the $L^2$- and $H^1$-error between the exact ground state and discrete minimizers, as well as error estimates for the ground state energy and the ground state eigenvalue. In particular, the achieved convergence rates for the energy and the eigenvalue are of $6$th order with respect to the mesh size on which the discrete LOD space is based, without making any additional regularity assumptions. These high rates justify the use of very coarse meshes, which significantly reduces the computational effort for finding accurate approximations of ground states. In addition, we include numerical experiments that confirm the optimality of the new theoretical convergence rates, both for smooth and discontinuous potentials.


Introduction
In this paper we consider the stationary Gross-Piteavskii equation (GPE), which can be written as an eigenvalue problem with eigenvector nonlinearity seeking u ∈ H 1 0 (Ω) with u L 2 (Ω) = 1 and λ ∈ R such that −∆u + V u + β|u| 2 u = λu.
Here, Ω ⊂ R d (with d = 1, 2, 3) is a computational domain, V (x) models a trapping potential and β ∈ R ≥0 is a constant. The perhaps most prominent application of the GPE is the modeling of so-called Bose-Einstein condensates [11,19,23,42]. Bose-Einstein condensates (BECs) are formed when a dilute gas of bosons is cooled to ultra-low temperatures so that almost all particles occupy the same quantum state, i.e. they become indistinguishable from each other and behave like one giant "macro particle". The first experimental realization of a BEC goes back to the Nobel prize-winning works by Cornell, Ketterle and Wieman [5,17]. Bose-Einstein condensates are of great interest in the study of macroscopic quantum effects such as superfluidity [2,40] and the closely related phenomenon of superconductivity [24].
In the context of BECs, β ≥ 0 models a repulsion parameter that depends on the number of bosons, as well as their type and mass. In this setting, the GPE has a positive spectrum, where the smallest eigenvalue is simple. This eigenvalue λ is called the ground state eigenvalue of the problem and it describes the chemical potential of a BEC in the lowest energy state. The corresponding ground state eigenfunction u is unique up to its sign. An equivalent characterization of the ground state u is to find a minimizer of an energy functional E (defined in equation (2.1) below) over the space of functions in H 1 0 (Ω) with the normalization constraint u L 2 (Ω) = 1.
Due to the complexity of physical experiments, computing the ground sate numerically is a very important task. Starting from the energy minimizing perspective, a numerical method has two essential components. The first component is an iterative method that allows to find minimizers of the energy E on the constrained manifold. Here there exists a variety of well-established schemes, where we exemplarily refer to methods based on self consistent field iterations (SCF) [13,14,18,43], methods based on discrete Sobolev gradient flows / Riemannian optimization [7,6,4,9,10,15,16,30,25,33,45] or the "J-method" [32,3]. We shall however not further discuss this aspect in our paper. Instead, we shall focus on the second important component of any numerical method for the computation of ground states: the choice of a suitable discrete space in which the energy functional is minimized. Standard choices are Lagrange finite element spaces (which can be often associated with finite difference discretizations) or Fourier spaces. These standard choices were analytically explored in [12]. In theses discrete spaces, each iteration of a scheme for finding a discrete minimizer of E requires the solution of a (typically elliptic) problem. Hence, the complexity of the final method depends naturally on the choice and the dimension of the discrete spaces.
In addition to the aforementioned standard choices, the usage of an alternative discrete space was proposed in [27]. Besides achieving high convergence orders, the proposed spaces have the advantage that they achieve their full potential already under minimal regularity assumptions on u and the potential V . Their construction is based on techniques from numerical homogenization, which build a low-dimensional generalized finite element space using a localized orthogonal decomposition (LOD). The LOD splits an ideal solution space into a low-dimensional approximation space and a high-dimensional remainder space. The low-dimensional space (to which we refer to as the LOD-space) is enriched with problem-specific features and admits a set of quasi-local basis functions that can be efficiently used in practical computations. In particular, due to the very good approximation properties and the quasi-locality of a basis, LOD-spaces can be easily used in a traditional Galerkin approach. Historically, the LOD was first introduced in [38] for elliptic multiscale problems with a rough coefficient. Since then, it has been further developed to suit a range of different problems, e.g., parabolic equations [36,35,37], equations describing wave phenomena, see, e.g., [1,41,28,21,34], non-linear problems [26,44], and more. In the context of nonlinear Schrödinger equations and GPEs, LOD techniques have been suggested in [27,31]. For a review of the LOD we refer to the textbook by Målqvist and Peterseim [39] and the recent survey article on numerical homogenization [3].
In this paper we revisit the LOD for the computation of ground states of the Gross-Pitaevskii equation (as introduced in [27]) in order to prove the higher order convergence rates that were numerically observed by the authors, but which have not yet been theoretically established. To be precise, the a priori error estimates proved in [27] for the GPE predict convergence rates of order O(H 2 ) for the H 1error, of order O(H 3 ) for the L 2 -error and the eigenvalue error and of order O(H 4 ) for the error in energy, where H is the "coarse" mesh size of the LOD-space. As already noted in [27], these rates seemed suboptimal when compared to numerical experiments. In fact, in this contribution we prove the optimal rates, which are of order O(H 3 ) for the H 1 -error, of order O(H 4 ) for the L 2 -error and (with a dramatic improvement) of order O(H 6 ) for the eigenvalue-error and the error in energy. For that we will not require any additional regularity assumptions and the results still hold for rough potentials V ∈ L 2+σ (Ω) with σ > 0.
Besides the obvious advantage of efficiently approximating ground states in the LOD space, our new results have direct implications for the simulation of the dynamics of Bose-Einstein condensates. In such a setting, the dynamics are described by nonlinear Schrödinger equations (NLS), where the initial values are typically ground states of the GPE (with respect to some modified configuration). In [31] it was recently proved that when the NLS is discretized with LOD spaces, then the energy is approximated and conserved with an accuracy of order O(H 6 ), provided that this accuracy can be already guaranteed for the initial value. Hence, this is exactly what we are establishing in this paper, as it justifies that the ground states computed in LOD-spaces can be straightforwardly used as initial values in an NLS without sacrificing the accuracy with which the energy is approximated over time.
Finally we note that in [27], an additional post-processing step on a finer mesh is suggested, which improves the quality of the approximation even further. We shall not focus on the post-processing in this paper. However, we note that the convergence orders can be further improved by such a technique.
The outline of the paper is as follows; in Section 2 the Gross-Pitaevskii equation is presented in more detail, in Section 3 we describe the LOD discretization technique, in Section 4 we present proofs of the optimal convergence rates in the H 1 -and L 2norm as well as for the energy and eigenvalue. Finally, in Section 5 we provide two numerical experiments that confirm the predicted convergence rates, both for smooth and discontinuous potentials.
We consider the Gross-Pitaevskii energy functional E : Here we make the following assumptions on the domain Ω, the trapping potential V and the interaction parameter β: (A1) Ω is a convex domain with polygonal boundary (for d = 1, 2, 3).
Observe that E is two times Fréchet-differentiable, where we denote the corresponding first and second derivatives by E and E respectively.
Under the above assumptions, we can seek the ground state of the energy E. The ground state is defined as a function u ∈ H 1 0 (Ω) that minimizes E(·) under the normalization constraint that Ω |u| 2 dx = 1, i.e.
The ground state is unique up to a sign and it can be equivalently characterized by an eigenvalue problem with eigenvector nonlinearity. The following result summarizes these analytical properties and a corresponding proof can be found in [12]. Furthermore, the unique positive minimizer u can be equivalently characterized by the Gross-Pitaevskii eigenvalue problem (GPE), which seeks the nonnegative eigenfunction u ∈ H 1 0 (Ω) with u = 1 and corresponding smallest eigenvalue λ > 0, such that The smallest eigenvalue λ and the ground state energy E(u) are connected through Note the nontrivial observation in Theorem 2.1 that if λ is the smallest eigenvalue of problem (2.3), then the corresponding eigenfunction u is a ground state of E. This means that minimizing the energy and minimizing the eigenvalue is equivalent.
We also note that (2.2) is a weak formulation that can be written using the Fréchet derivative of E. More precisely, with we observe that any eigenpair (λ, u) of (2.3) satisfies This formulation expresses that the eigenvalue problem should be seen as finding the critical points of the (even, positiv and convex) functional E on the L 2 -sphere S := {φ ∈ H 1 0 (Ω)| φ = 1}. Consequently, we also know that the spectrum of the differential operator is real and unbounded.
In the more general form of the Gross-Pitaevskii eigenvalue problem, where a rotational term is included, the eigenfunctions u are typically complex-valued, but the spectrum remains real. We write |u| 2 u instead of u 3 (which is the same for real functions), to work with the standard notation for the Gross-Pitaevskii equation that remains valid also in complex-valued settings.
Using the equivalence of the energy minimization perspective (2.2) and the GPE (2.3) (respectively (2.5)), we can easily prove additional regularity for the ground state u that will be useful in the convergence analysis. Lemma 2.2. Assume (A1)-(A3). Let u be a minimizer of (2.2) and let λ > 0 denote the corresponding ground state eigenvalue in the sense of (2.3), then u ∈ C 0 (Ω) ∩ H 2 (Ω) and |u| 2 u ∈ H 1 0 (Ω) ∩ H 2 (Ω) and we have the following bounds where C is a generic constant that only depends on Ω and d. Note that u H 1 (Ω) can be bounded again by √ λ.
Proof. In the following we write a b to abbreviate a ≤ C b for a constant C that only depends on Ω and d. In the first step we write If g ∈ L q/2 (Ω), for q > d, it follows that u ∈ C 0 (Ω), see e.g. [22,Theorem 8.30]. To verify that this condition is fulfilled, we define q such that 6q/(12 − q) = 2 + σ, where σ is as in (A2). For d = 1 and d = 2, this gives us q = 3 and consequently also q > d as desired. For d = 3, the condition q > d is equivalent to which is consistent with the assumption σ > 0 in (A2). Without loss of generality, we further assume σ ≤ 1 so that q ≤ 4 for d = 3. For q as introduced above, we need to prove now that g ∈ L q/2 (Ω). For that we use Hölder's inequality and the Sobolev embedding where we exploited u = 1. The term V L 2+σ (Ω) is bounded by (A2) and for the last term we get Consequently, we verified g ∈ L q/2 (Ω) for q > d and we have u ∈ C 0 (Ω). It is worth to mention that u ∈ C 0 (Ω) was already proved in [12], however, we could not verify the result under the stated regularity assumptions which slightly differ from ours. Therefore we decided to explicitly elaborate the argument here. Now that we have u ∈ C 0 (Ω), we can use this result to show that we indeed even have g ∈ L 2 (Ω), as . Using standard elliptic regularity theory for equation (2.8) we conclude that u ∈ H 2 (Ω) with the estimate . Hence, the bound in (2.6) follows.
For (2.7) note that where the multi-index α denotes the partial derivatives. For |α| = 0 the term reduces to (2.9). For the first order derivatives, i.e., for |α| = 1, we get where Hölder's inequality and the Sobolev embedding gives . For the second order derivatives, i.e. for |α| = 2, we have where α 1 , α 2 denotes first order derivatives such that α = α 1 +α 2 . For the first term (2.11) we have using Hölder's inequality and the Sobolev embedding once again . For the second term in (2.11) we get . This gives . Combining all the terms and using u H 1 (Ω) ≤ u H 2 (Ω) we achieve (2.7).

Discretization
Approximations of the ground state can be found by minimizing the energy E over a discrete space V δ subject to the L 2 -normalization constraint. The following abstract approximation result was established in [12, Theorem 1 and Remark 2] and will be a central ingredient for our proofs.
Theorem 3.1. Assume (A1)-(A3) and let u be the ground state of (2.2) (up to sign). Furthermore, let (V δ ) δ>0 be a family of finite-dimensional subspaces of H 1 0 (Ω) with the property that with the property (u, u δ ) ≥ 0. Then there exist (problem-dependent) generic constants C, c 1 , c 2 > 0 such that i.e., u δ is a H 1 -quasi best approximation to u; and it also holds Finally, for the L 2 -error it holds where for ω ∈ H −1 (Ω), the function The results established in [12] also contain an abstract approximation result for the eigenvalue, which is however not optimal for all discrete spaces and needs to be therefore revisited individually depending on the considered discrete setting.
Following [27], our next goal is the construction of a suitable discrete space that allows for very high convergence orders under minimal regularity assumptions. For that we start from a conventional conforming finite element discretization based on a (coarse) quasi-uniform triangulation T H of the domain Ω. Here H denotes the mesh size parameter, i.e. H = diam(T ), for T ∈ T H . With this, let V H be a classical P 1 finite element space on T H with The classical finite element method seeks an approximation u H ∈ V H with u H = 1, and where the corresponding eigenvalue is given by . Since E is weakly lower semi-continuous on H 1 0 (Ω) and bounded from below by zero, discrete minimizers always exist on finite dimensional spaces. However, it is typically not clear if the solution to (3.3) is also unique up to sign, see e.g. the discussion in [12]. Independent of that uniqueness, it can be shown that for any minimizer u H with (u H , u) ≥ 0 (i.e. consistent with the sign of u) it holds where u is the unique (up to sign) ground state and C is a constant independent of H. The space that we construct in the next subsection will have the same dimension as V H , but it will be able to boost the convergence rates to 3rd order for the H 1 -error, 4th order for the L 2 -error and 6th order for eigenvalue and energy error.
3.1. Localized Orthogonal Decomposition. In the following we summarize the technique proposed in [27] to make the paper self-contained. Let P H : H 1 0 (Ω) → V H denote the classical L 2 -projection onto V H and let V f := ker(P H ) = {v ∈ H 1 0 (Ω)| P H (v) = 0}, be the kernel of the projection. Due to H 1 -stability of the L 2 -projection on quasiuniform meshes (cf. [8]), the space V f is a closed subspace of H 1 0 (Ω).
We expect V f to contain fine scale details from H 1 0 (Ω) that are not captured by the coarse space V H . With this definition we obtain the following L 2 -orthogonal splitting Next, we define an inner product a(·, ·) on H 1 0 (Ω), that is based on the linear terms in the GPE (2.3), by a(v, w) := (∇v, ∇w) + (V v, w).

This inner product induces another orthogonal splitting, namely
respect to a(·, ·). Note that, by construction, the V LOD space has the same dimension as V H . This is precisely the space in which we will look for discrete minimizers. This is fixed in the following definition.
Definition 3.2 (LOD ground state approximation). The LOD ground state approximation is obtained by using the V LOD space in (2.2). That is, we seek u LOD ∈ V LOD such that (u LOD , u) ≥ 0, u LOD = 1, and where the corresponding discrete eigenvalue is given by Remark 3.3 (Practical realization). The practical realization of (3.4) requires some remarks. First of all, the condition (u LOD , u) ≥ 0 is never enforced in the numerical method, as it does not matter if u LOD converges to the positive or the negative ground state, the respective approximation properties are identical. Second, for the construction of V LOD it is necessary to solve a set of saddle point problems that determine the basis functions of V LOD . For an efficient realization, this includes a localization step which is analytically justified by an exponential decay of the basis function in units of H (cf. [29,38]). Furthermore, the space V f needs a discrete representation on a finer mesh of size h < H. Comprehensive details on all these implementation details are given in [20].
The influence of these additional approximations (i.e. truncation/localization and fine scale representation of the LOD basis functions) has been studied extensively in the literature and we note that they do not lead to a reduction of the convergence orders of the final method when chosen appropriately. For the Gross-Pitaevskii eigenvalue problem, the influence has been analytically studied in [27].
Further practical tricks for an efficient treatment of the nonlinear term are discussed in [31].
Finally we also stress that aside from the construction of the LOD space, a practical realization also requires an iterative method for finding a discrete minimizer. Here various techniques are possible, where we refer to the introduction for a corresponding literature survey.
In the light of Remark 3.3, all the proofs below are based on the ideal setting of Definition 3.2, i.e. the setting without localization and fine scale representation. This simplifies both the arguments and the notation. By considering localization, an additional term depending on the size of the localization patches would appear in the error estimates. The implications of truncation and fine scale discretization are already thoroughly discussed in [27] for the GPE, which is why we omit this part in our paper. This allows us to keep the presentation short and to focus on the main novelty of this contribution, namely the arguments that actually improve the convergence orders in the proofs.
We have the following main result, which we shall prove in Section 4.
where C depends on V , λ, β, Ω, u L ∞ (Ω) and the mesh regularity of T H .
Compared to the original result proved in [27], the rates for the L 2 and H 1 -error are both improved by one order each, the rate for the energy by two orders and the rate for the eigenvalue is dramatically improved by three orders. However, we also stress that the equation considered in [27] is slightly more general, involving another (potentially discontinuous) coefficient in the kinetic term. Therefore the rates in [27] are optimal for this generalized setting, but only suboptimal for the GPE. Finally, we also note that the rates in (3.5) can be even further improved significantly by using a post-processing technique as in [27].

Proof of Optimal Convergence Rates
In this section we prove the optimal convergence rates of order O(H 3 ) in the H 1norm, O(H 4 ) in the L 2 -norm, and O(H 6 ) for both the energy and the eigenvalue. The error bounds presented below all depend on some constant C, which may depend on the ground state eigenvalue λ, the potential V , the interaction constant β, and the domain Ω, but not on the mesh size H.
In the analysis below we shall make use of the a(·, ·)-orthogonal projection P LOD : The following approximation properties of the LOD space are well known and are explicitly stated and proved in [31, Section 2.1].
The constant C only depends on Ω and the mesh regularity of T H .

4.1.
Convergence in H 1 -norm. We start with the H 1 -error estimate which follows from the abstract H 1 -best approximation result in Theorem 3.1, combined with the approximation properties in Lemma 4.1 and the regularity result in Lemma 2.2. We have: Proof. We note that

4.2.
Convergence in L 2 -norm. Next, we turn our attention to proving the L 2error estimate which reads as follows. The proof is similar to [27,Section 6.3] with several simplifying modifications that keep the proof short.
We analyze the terms on the right hand side one by one. Hölder's inequality and the Sobolev inequality u L 6 (Ω) ≤ C u H 1 (Ω) (for d ≤ 3) give us for the first term For the second term we readily have λψ eLOD ≤ Cλ e LOD . For the third term we use (2.9) and u = 1 to get With u = 1, the last term is estimated as Combining the estimates for all terms yields

4.4.
Convergence of Eigenvalues. The optimal convergence rates for the eigenvalue involve the most substantial changes compared to [27]. In the first step, we prove the following auxiliary result, which established optimal convergence rates in the H −2 -norm.
We conclude with (4.3) that For II 2 we use Lemma 2.2 For the last terms we use Lemma 2.2 again together with the regularity bound ψ ω H 2 (Ω) ≤ C ω to conclude We deduce II ≤ CH 6 ω H 2 (Ω) .
For the third and fourth term in (4.9) we use the error bound in the L 2 -norm from Proposition 4.3, the regularity of ψ ω , and Lemma 2.2 to derive Combining everything, we obtain (ω, u − u LOD ) ≤ CH 6 ω H 2 (Ω) , which gives With this we are ready to prove the final estimate for the eigenvalue, which also finishes the proof of Theorem 3.4. Proposition 4.6. Let λ and λ LOD be the corresponding eigenvalues to u and u LOD in (2.4) and (3.4) respectively. It holds that |λ − λ LOD | ≤ CH 6 .
Proof. We use the definition of the eigenvalues to estimate where the first term is bounded by Conclusion 4.4. It remains to bound the second term. We have For the first term we use Proposition 4.2 to achieve The second term splits into . We note that by Hölder's inequality and Sobolev embeddings we have the bound This relies on the fact that u LOD H 1 (Ω) ≤ C, which can simply be deduced by using the error bound in H 1 -norm Thus, using Lemma 4.5 and Proposition 4.2 we deduce II ≤ CH 6 , which completes the proof.

Numerical Experiments
In this section we perform two numerical experiments to validate the theoretical convergence rates for both a smooth and a discontinuous potential V . The numerical experiments are performed for the ideal method, without localization of the corrections, to eliminate any potential reduction in convergence rate coming from the size of the localization patches. We emphasize that these experiments are only to verify that the derived rates are indeed optimal for the Gross-Pitaevskii equation and in a practical application the localization should always be implemented. For further details regarding localization we refer to [27] where the method was first suggested.

Smooth potential.
For the smooth potential we consider a setup similar to [9]. The computational domain is given by the square Ω = [−6, 6] × [−6, 6], the potential is set to the smooth (harmonic) function V (x, y) = 0.5(x 2 + y 2 ), and the interaction parameter to β = 100.
The resulting minimizing problem (3.4) is solved by using the iterative scheme, referred to as normalized gradient flow, presented in [9].
A reference solution is computed on mesh a fine mesh of size h = 2 −4 . The LOD approximations u LOD are computed for meshes of decreasing size H = 2 1 , 2 0 , 2 −1 , 2 −2 and then compared to the reference solution. The resulting convergence rates for the H 1 -norm, the L 2 -norm, the energy, and the eigenvalue are plotted in Figure 1. The plots clearly confirm the convergence orders predicted by Theorem 3.4.
The potential is set to a checkerboard with squares of size 2 −2 which results in 24 × 24 squares on the domain. The values at the squares alternate between 0 and 1, see Figure 2.
A reference solution is computed on a mesh of size h = 2 −5 . The triangles in the mesh are aligned with the discontinuities in the potential, so that the potential can be computed exactly on each triangle to avoid additional numerical errors. The LOD approximations u LOD are computed for meshes of decreasing size    H = 2 0 , 2 −1 , 2 −2 , 2 −3 and then compared to the reference solution. The resulting convergence rates for the H 1 -norm, the L 2 -norm, the energy, and the eigenvalue are plotted in Figure 3. Again, the plots clearly confirm the predicted convergence orders.