Abstract
This work presents a new methodology for computing ground states of BoseEinstein condensates based on finite element discretizations on two different scales of numerical resolution. In a preprocessing step, a lowdimensional (coarse) generalized finite element space is constructed. It is based on a local orthogonal decomposition of the solution space and exhibits high approximation properties. The nonlinear eigenvalue problem that characterizes the ground state is solved by some suitable iterative solver exclusively in this lowdimensional space, without significant loss of accuracy when compared with the solution of the full fine scale problem. The preprocessing step is independent of the types and numbers of bosons. A postprocessing step further improves the accuracy of the method. We present rigorous a priori error estimates that predict convergence rates for the ground state eigenfunction and for the corresponding eigenvalue without preasymptotic effects; being the coarse scale discretization parameter. Numerical experiments indicate that these high rates may still be pessimistic.
TwoLevel discretization techniques for ground state computations of BoseEinstein condensates
[2em]
Patrick Henning
[2em]
October 9, 2018
Keywords eigenvalue, finite element, GrossPitaevskii equation, numerical upscaling, twogrid method, multiscale method
AMS subject classifications 35Q55, 65N15, 65N25, 65N30, 81Q05
1 Introduction
BoseEinstein condensates (BEC) are formed when a dilute gas of trapped bosons (of the same species) is cooled down to ultralow temperatures close to absolute zero [10, 19, 22, 38]. In this case, nearly all bosons are in the same quantum mechanical state, which means that they loose their identity and become indistinguishable from each other. The BEC therefore behaves like one ’super particle’ where the quantum state can be described by a single collective wave function . The dynamics of a BEC can be modeled by the timedependent GrossPitaevskii equation (GPE) [26, 31, 37], which is a nonlinear Schrödinger equation given by
(1) 
Here, denotes the atomic mass of a single boson, the number of bosons (typically in the span between and ), is the reduced Plank’s constant and is an external trapping potential that confines the system. The nonlinear term in the equation describes the effective twobody interaction between the particles. If the scattering length is positive, the interaction is repulsive, if it is negative the interaction is attractive. For there is no interaction and (1) becomes the Schrödinger equation. The parameter changes according to the considered species of bosons. We only consider the case in this paper. We are mainly interested in the ground state solution of the problem. This stationary state of the BEC is of practical relevance, e.g., in the context of atom lasers [35, 30, 41]. The ansatz , with the unknown chemical potential of the condensate and a proper nondimensionalization , reduces (1) to the timeindependent GPE
where denotes the dimensionless length unit and where denotes the accordingly rescaled potential (see, e.g., [8] for a derivation of the timeindependent GPE). The ground state of the BEC is the lowest energy state of the system and is therefore stable. It minimizes the corresponding energy
amongst all normalized functions. For any normalized minimizer , is the smallest eigenvalue of the GPE. In this paper, we shall focus on the computation of this ground state eigenvalue. Eigenfunctions whose energies are larger than the minimum energy are called excited states of the BEC and are not stable in general but may satisfy relaxed concepts of stability such as metastability (see [36]). Numerical approaches for the computation of ground states of a BEC typically involve an iterative algorithm that starts with a given initial value and diminishes the energy of the density functional in each iteration step. Different methodologies are possible: methods related to normalized gradientflows [5, 3, 1, 2, 5, 7, 24, 6, 9, 20], methods based on a direct minimization of the energy functional [8, 11], explicit imaginarytime marching [32], the DIIS method (direct inversion in the iterated subspace) [40, 16], or the Optimal Damping Algorithm [14, 12]. We emphasize that, in any case, the dimensionality of the underlying space discretization is the crucial factor for computational complexity because it determines the cost per iteration step. The aim of this paper is to present a lowdimensional space discretization that reduces the cost per step and, hence, speeds up the iterative solution procedure considerably. In the literature, there are only a few contributions on rigorous numerical analysis of space discretizations of the GPE. In particular, explicit orders of convergence are widely missing. In [44, 17], Zhou and coworkers proved the convergence of general finite dimensional approximations that were obtained by minimizing the energy density in a finite dimensional subspace of . This justifies, e.g., the direct minimization approach proposed in [8]. The iteration scheme is not specified and not part of the analysis. The results of Zhou were generalized by Cancès, Chakir and Maday [13] allowing explicit convergence rates for finite element approximations and Fourier expansions. Apriori error estimates for a conservative CrankNicolson finite difference (CNFD) method and a semiimplicit finite difference (SIFD) method were derived by Bao and Cai [4].
In this work, we propose a new space discretization strategy that involves a preprocessing step and a postprocessing step in standard finite element spaces. The preprocessing step is based on the numerical upscaling procedure suggested by two of the authors [33] for linear eigenvalue problems. In this step, a lowdimensional approximation space is assembled. The assembling is based on some local orthogonal decomposition that incorporates problemspecific information. The constructed space exhibits high approximation properties. The nonlinear problem is then solved in this lowdimensional space by some standard iterative scheme (e.g., the ODA [14]) with very low cost per iteration step. The postprocessing step is based on the twogrid method suggest by Xu and Zhou [42]. We emphasize that both, pre and postprocessing, involve only the solution of linear elliptic Poissontype problems using standard finite elements. We give a rigorous error analysis for our strategy to show that we can achieve convergence orders of for the computed eigenvalue approximations without any preasymptotic effects. We do not focus on the iterative scheme that is used for solving the discrete minimization problem. The various choices previously mentioned, e.g., the ODA [14] are possible. Our new strategy is particularly beneficial in experimental setups with different types of bosons, because the results of the preprocessing step can be reused over and over again independent of . Similarly, the data gained by preprocessing can be recycled for the computation of excited states. Other applications include setups with potentials that oscillate at a very high frequency (e.g., to investigate Josephson effects [41, 43]). Here, normally very fine grids are required to resolve the oscillations, whereas our strategy still yields good approximations in low dimensional spaces and, hence, reduces the costs within the iteration procedure tremendously.
2 Model problem
Consider the dimensionless GrossPitaevskii equation in some bounded Lipschitz domain where . Since ground state solutions show an extremely fast decay (typically exponential), the restriction to bounded domains and homogeneous Dirichlet condition are physically justified. We seek (in the sense of distributions) the minimal eigenvalue and corresponding normalized eigenfunction with
The underlying data satisfies the following assumptions:

If , the domain is an interval. If (resp. ), has a polygonal (resp. polyhedral) boundary.

The diffusion coefficient is a symmetric matrixvalued function with uniform spectral bounds ,
(2) 
is nonnegative (almost everywhere).

is nonnegative.
The weak solution of the GPE minmizes the energy functional given by
Problem 2.1 (Weak formulation of the GrossPitaevskii equation).
Find such that a.e. in , , and
It is wellknown (see, e.g., [31] and [13]) that there exists a unique solution of Problem 2.1. This solution is continuous in and positive in . The corresponding eigenvalue of the GPE is real, positive, and simple. Observe that the eigenpair satisfies
for all . Moreover, is the smallest amongst all possible eigenvalues and satisfies the a priori bound .
3 Discretization
This section recalls classical finite element discretizations and presents novel twogrid approaches for the numerical solution of Problem 2.1. The existence of a minimizer of the functional in discrete spaces is easily seen. However, uniqueness does not hold in general. We note that unlike claimed in [44] the uniqueness proof given in [31] does not generalize to arbitrary subspaces of the original solution space.
Remark 3.1 (Existence of discrete solutions [13]).
Let denote a finite dimensional, nonempty subspace of , then there exists a minimizer with , , and
If represents a dense family of such subspaces, then any sequence of corresponding minimizers with converges to the unique solution of Problem 2.1.
3.1 Standard Finite Elements
We consider two regular simplicial meshes and of . The finer mesh is obtained from the coarse mesh by regular mesh refinement. The discretization parameters represent the mesh size, i.e., (resp. ) for (resp. ) and (resp. ). For , let
denote the set of piecewise affine functions. Classical conforming finite element spaces are then given by
Note that on the fine discretization scale, a different choice of polynomial degree, e.g., piecewise quadratic functions, is possible. This would be a better choice for smooth data that allows for a regular ground state. Our method and its analysis essentially require the inclusion . The discrete problem on the fine grid reads as follows.
Problem 3.2 (Reference finite element discretization on the fine mesh).
Find with , and
(3) 
The corresponding eigenvalue is given by .
According to Remark 3.1, is not determined uniquely in general. Moreover, is not necessarily the smallest eigenvalue of the corresponding discrete eigenvalue problem. In what follows, refers to an arbitrary solution of Problem 3.2. It will serve as a reference to compare further (cheaper) numerical approximations with. The accuracy of has been studied in [13]. Under the assumption of sufficient regularity, optimal orders of convergence are obtained (cf. (14)).
3.2 Preprocessing motivated by numerical homogenization
The aim of this paper is to accurately approximate the finescale reference solution of Problem 3.2 within some lowdimensional subspace of . For this purpose, we introduce a twogrid upscaling discretization that was initially proposed in [34] for the treatment of multiscale problems. The framework has been applied to nonlinear problems in [27], to linear eigenvalue problems in [33] and in the context of Discontinuous Galerkin [23] and Partition of Unity Methods [28]. This contribution aims to generalize and analyze the methodology to the case of an eigenvalue problem with an additional nonlinearity in the eigenfunction. We emphasize that the coexistence of two difficulties, the nonlinear nature of the eigenproblem itself and the additional nonlinearity in the eigenfunction, requires new essential ideas far beyond simply plugging together existing theories for the isolated difficulties.
Let denote the set of interior vertices in . For we let denote the corresponding nodal basis function with and for all . We define a weighted Clémenttype interpolation operator (c.f. [15])
(4) 
It is easily shown by Friedrichs’ inequality and the Sobolev embedding (for ) that
defines a scalar product in and induces a norm on which is equivalent to the standard norm. By means of the interpolation operator defined in (4), we construct an orthogonal decomposition of the space into a lowdimensional coarse space (with favorable approximation properties) and a highdimensional residual space . The residual or ’fine’ space is the kernel of the interpolation operator restricted to ,
(5.a) 
The coarse space is simply defined as the orthogonal complement of in with respect to . It is characterized via the orthogonal projection onto the fine space given by
By defining , the coarse space is given by
(5.b) 
A basis of is given by with . With this definition we obtain the splitting
(5.c) 
Some favorable properties of the decomposition, in particular its quasiorthogonality, are discussed in Section 6.2. The minimization problem in the lowdimensional space reads as follows.
Problem 3.3 (Preprocessed approximation).
Find with , and
The corresponding eigenvalue in is given by .
Remark 3.4 (Practical aspects of the decomposition).

The assembly of the corresponding finite element matrices requires only the evaluation of , i.e., the solution to one linear Poissontype problem per coarse vertex. This can be done in parallel. Section 3.3 below will show that these linear problems may be restricted to local subdomains centered around the coarse vertices without loss of accuracy. Hence, even in a serial computing setup, the complexity of solving all corrector problems is equivalent (up to factor ) to the cost of solving one linear Poisson problem on the fine mesh.

The preprocessing step is independent of the parameter which characterizes the species of the bosons. Hence, the method becomes considerably cheaper when experiments need to be carried out for different types and numbers of bosons. A similar argument applies to variations on the trapping potential . Provided that this trapping potential is an element of (in practical applications it is usually even harmonic and admits the desired regularity) the bilinear form (and the associated constructions of and ) can be restricted to the second order term without a loss in the expected convergence rates stated in Theorems 4.1 and 4.2 below. The trapping potential may then be varied without affecting the preprocessed space .

Once the coarse space has been assembled it can also be reused in computations of larger eigenvalues (i.e., not only in the ground state solution).
3.3 Sparse approximations of
The construction of the coarse space is based on fine scale equations formulated on the whole domain which makes them expensive to compute. However, [34] shows that decays exponentially fast away from . We specify this feature as follows. Let denote the localization parameter, i.e., a new discretization parameter. We define nodal patches of coarse grid layers centered around the node by
(6)  
There exists depending on the contrast but not on mesh sizes and fast oscillations of such that for all for all vertices and for all , it holds
(7) 
This result motivates the truncation of the computations of the basis functions to local patches . We approximate from (5.a)(5.c) with such that
(8) 
This yields a modified coarse space with a local basis
(9) 
The number of nonzero entries of the corresponding finite element matrices is proportional to (note that we expect nonzero entries without the truncation). Due to the exponential decay, the very weak condition implies that the perturbation of the ideal method due to this truncation is of higher order and forthcoming error estimates in Theorems 4.1 and 4.2 remain valid. We refer to [34] for details and proofs. The modified localization procedure from [29] with improved accuracy and stability properties may also be applied.
3.4 Postprocessing
Although and will turn out to be highly accurate approximations of the unknown solution , the orders of convergence can be improved even further by a simple postprocessing step on the fine grid. The postprocessing applies the twogrid method originally introduced by Xu and Zhou [42] for linear elliptic eigenvalue problems to the present equation by using our upscaled coarse space on the coarse level.
Problem 3.5 (Postprocessed approximation).
Find with
for all . Define .
Let us emphasize that this approach is different from [18], where the postprecessing problem has a different structure and where classical finite element spaces are used on both scales.
4 Apriori error estimates
This section presents the apriori error estimates for the preprocessed/upscaled approximation with and without the postprocessing step. Throughout this section, denotes the solution of Problem 2.1, the solution of reference Problem 3.2, the solution of Problem 3.3 and the postprocessed solution of Problem 3.5. The notation abbreviates with some constant that may depend on the space dimension , , , , , , and interior angles of the triangulations, but not on the mesh sizes and . In particular it is robust against fast oscillations of and .
Theorem 4.1 (Error estimates for the preprocessed approximation).
Assume that . For and as above, it holds
(10) 
For sufficiently small (in the sense of Cancès et al. [13]), we also have
(11) 
Proof.
The proof is postponed to Section 6.3. ∎
The additional postprocessing improves, roughly speaking, the order of accuracy by one.
Theorem 4.2 (Error estimates for the postprocessed approximation).
Assume that is sufficiently small. The postprocessed approximation and the postprocessed eigenvalue satisfy:
(12)  
(13) 
The constant behaves roughly like and can be extracted from the proofs in Section 6.4.2.
Proof.
The proof is postponed to Section 6.4. ∎
Let us emphasize that both theorems remain valid for replaced with its sparse approximation (cf. Section 3.3) for moderate localization parameter .
We shall discuss the behavior of the finescale errors and . Recall from [13] that for a bounded domain with polygonal Lipschitzboundary, , and sufficiently small , the fine scale error satisfies the optimal estimate
(14) 
The proof in [13] is for constant and hyperrectangle but it is easily checked that the estimates remain valid for any bounded domain with polygonal Lipschitzboundary and . Under these assumptions our a priori estimates for the postprocessed approximation of the ground state eigenvalue summarize as follows
Hence, in this regular setting, the choice ensures that the loss of accuracy is negligible when compared to the accuracy of the expensive full fine scale approximation . However, with regard to the numerical experiment in Section 5.1 below, this choice might be pessimistic.
Moreover, note that the fine scale error depends crucially on higher Sobolev regularity of the solution whereas our estimates for the coarse scale error require only minimal regularity that holds under the assumption (a)–(d) in Section 2. Thus, we believe that in a less regular setting, even coarser choices of relative to will balance the discretization errors on the coarse and the fine scale.
5 Numerical experiments
Any numerical approach for the computation of ground states of a BEC involves an iterative algorithm that starts with a given initial value and diminishes the energy of the density functional in each iteration step. In this contribution, we use the Optimal Damping Algorithm (ODA) originally developed by Cancès and Le Bris [14, 12] for the HartreeFock equations, since it suits our preprocessing framework. The ODA involves solving a linear eigenvalue problem in each iteration step. However, after preprocessing these linear eigenvalue problems are very low dimensional and the precomputed basis of can be reused for each of these problems making the iterations extremely cheap. The approximations produced by the ODA are known to rapidly converge to a solution of the discrete minimization problem (see [21] and [12] for a proof in the setting of the HartreeFock equations). All subsequent numerical experiments have been performed using MATLAB.
5.1 Numerical results for harmonic potential
In this section, we choose the smooth experimental setup of [13, Section 4, p. 109 and Fig. 2 (bottom)], i.e., , , , and with homogeneous Dirichlet boundary condition. Our method depends basically on three parameters, the coarse mesh size , the fine mesh size , and the localization parameter (cf. Section 3.3 and [29]). In all computations of this section we couple to the coarse mesh size by choosing . This choice is made such that the error of localization is negligible when compared with the errors committed be the fine scale discretization and the upscaling. All approximations are computed with the ODA method as presented in [21, Section 2] with accuracy parameter .
Comparison with full fine scale approximation
In the first experiment, we consider uniform coarse meshes with mesh width parameters of . The fine mesh for the pre and postprocessing has width and remains fixed. We study the error committed by coarsening from a fine scale to several coarse scales , i.e., we study the distance between the ground state of Problem 3.2 and either the coarse scale approximation of Problem 3.3 (with underlying finescale ) or its postprocessed version of Problem 3.5. Our theoretical results do not allow predictions about the coarsening error. Most likely, this is an artifact of our theory and we conjecture that and its coarse approximations and are in fact superclose in the sense of
(15)  
This assertion is true in the limit . Section 5.1.2 supports numerically the assertion for positive . Figure 1 reports the numerical results. Observe that the experimental rates with respect to displayed in the figures are in fact better than the rates indicated by Theorems 4.1–4.2 and conjectured in (15). The reason could be the high regularity of the underlying (exact) solution . We do not exploit additional regularity in our error analysis. Similar observations have been made for the linear eigenvalue problem; see [33, Remark 3.3] for details and some justification of higher rates under additional regularity assumptions. Our implementation is not yet adequate for a fair comparison with regard to computational complexity and computing times between standard fine scale finite elements and our twolevel techniques. However, to convince the reader of the potential savings in our new approach, let us mention that the number of iterations of the ODA were basically the same for both approaches in all numerical experiments. This statement applies as well to more challenging setups with larger values of (see, e.g., Section 5.2 below) where ODA needs many iterations to fall below some prescribed tolerance. We, hence, conclude that the actual speedup of our approach is truly reflected by the dimension reduction from to up to the overhead induced by slightly denser (but still sparse) finite element matrices on the coarse level.
Comparison with highresolution numerical approximation
In the second experiment we investigate the role of the fine scale parameter . We consider uniform coarse meshes with mesh width parameters and uniform fine meshes for for pre and postprocessing computations. The error between the exact eigenvalue and coarse approximations and is estimated via a highresolution numerical solution on a mesh of width . The results are reported in Figure 2. For the sake of clarity, we show eigenvalue errors only. We conclude that it would have been sufficient to choose to achieve the accuracy of by our coarse approximation scheme with postprocessing.
5.2 Numerical results for discontinuous periodic potential
This section addresses the case of a BEC that is trapped in a periodic potential. Periodic potentials are of special interest since they can be used to explore physical phenomena such as Josephson oscillations and macroscopic quantum selftrapping of the condensate (c.f. [41, 43]). Here we use a potential that describes a periodic array of quantum wells that can be experimentally generated by the interference of overlapping laser beams (c.f. [39]).
Let , , and . Given and , define
and the potential .
Consider the same numerical setup as in Section 5.1.1 (i.e., we draw our attention again to the coarsening error ) with the exception that we were able to reduce the localization parameter without affecting the best convergence rates possible. Figure 3 reports the errors between the finescale reference discretization and our coarse approximations.
6 Proofs of the main results
In this section we are concerned with proving the main theorems.
6.1 Auxiliary results
An application of [13, Theorem 1] shows that and both converge to in , which guarantees stability.
Remark 6.1 (Stability of discrete approximations).
For sufficiently small we have
(16)  
(17) 
The same results hold for replaced by and replaced by for and sufficiently small.
The bound (16) is obvious using and the convergence which guarantees . Estimate (17) directly follows from the definitions of and which gives us .
Remark 6.2 (bound).
The solution of Problem 2.1 is in . This follows from the uniqueness of which shows that it is also the unique solution of the linear elliptic problem
where . Standard theory for linear elliptic problems (c.f. [25, Theorem 8.15, pp. 189–193]) then yields the existence of a constant only depending on , and such that
(18) 
6.2 Properties of the coarse space
Recall the local approximation properties of the weighted Clémenttype interpolation operator defined in (4),
(19) 
for all . Here, is a generic constant that depends only on interior angles of but not on the local mesh size and . Furthermore, for all and for all it holds
(20) 
where and is given by (4).
Lemma 6.3 (Properties of the decomposition).
The decomposition of into and (stated in Section 3.2) is orthogonal, i.e.,
(21) 
The decomposition of in and is orthogonal
(22) 
and quasiorthogonal in the sense that
(23) 
Proof.
The proof is verbatim the same as in [33]. ∎
The following lemma estimates the error of the bestapproximation in the modified coarse space . The lemma is also implicitly required each time that we use the abstract error estimates stated in [13, Theorem 1]. These estimates require a family of finitedimensional spaces that is dense in . This density property is implied by the following lemma.
Lemma 6.4 (Approximation property of ).
For any given with it holds
Proof.
Given , define (since ) and let denote the corresponding finite element approximation, i.e.,
With , Galerkinorthogonality leads to
This, the triangle inequality, and norm equivalences readily yield the assertion.∎
Next, we show that there exists an element in the space that approximates in the energynorm with an accuracy of order O.
Lemma 6.5 (Stability and approximability of the reference solution).
Let solve Problem 3.2. Then it holds