This work presents a new methodology for computing ground states of Bose-Einstein condensates based on finite element discretizations on two different scales of numerical resolution. In a pre-processing step, a low-dimensional (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 non-linear eigenvalue problem that characterizes the ground state is solved by some suitable iterative solver exclusively in this low-dimensional space, without significant loss of accuracy when compared with the solution of the full fine scale problem. The pre-processing step is independent of the types and numbers of bosons. A post-processing 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 pre-asymptotic effects; being the coarse scale discretization parameter. Numerical experiments indicate that these high rates may still be pessimistic.

Two-Level discretization techniques for ground state computations of Bose-Einstein condensates


Patrick Henning1, Axel Målqvist2 3, Daniel Peterseim4 5


October 9, 2018

Keywords eigenvalue, finite element, Gross-Pitaevskii equation, numerical upscaling, two-grid method, multiscale method

AMS subject classifications 35Q55, 65N15, 65N25, 65N30, 81Q05

1 Introduction

Bose-Einstein condensates (BEC) are formed when a dilute gas of trapped bosons (of the same species) is cooled down to ultra-low 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 time-dependent Gross-Pitaevskii equation (GPE) [26, 31, 37], which is a nonlinear Schrödinger equation given by


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 two-body 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 time-independent GPE

where denotes the dimensionless length unit and where denotes the accordingly rescaled potential (see, e.g., [8] for a derivation of the time-independent 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 gradient-flows [5, 3, 1, 2, 5, 7, 24, 6, 9, 20], methods based on a direct minimization of the energy functional [8, 11], explicit imaginary-time 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 low-dimensional 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. A-priori error estimates for a conservative Crank-Nicolson finite difference (CNFD) method and a semi-implicit finite difference (SIFD) method were derived by Bao and Cai [4].

In this work, we propose a new space discretization strategy that involves a pre-processing step and a post-processing step in standard finite element spaces. The pre-processing step is based on the numerical upscaling procedure suggested by two of the authors [33] for linear eigenvalue problems. In this step, a low-dimensional approximation space is assembled. The assembling is based on some local orthogonal decomposition that incorporates problem-specific information. The constructed space exhibits high approximation properties. The non-linear problem is then solved in this low-dimensional space by some standard iterative scheme (e.g., the ODA [14]) with very low cost per iteration step. The post-processing step is based on the two-grid method suggest by Xu and Zhou [42]. We emphasize that both, pre- and post-processing, involve only the solution of linear elliptic Poisson-type 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 pre-asymptotic 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 pre-processing step can be reused over and over again independent of . Similarly, the data gained by pre-processing 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 Gross-Pitaevskii 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 matrix-valued function with uniform spectral bounds ,

  • is non-negative (almost everywhere).

  • is non-negative.

The weak solution of the GPE minmizes the energy functional given by

Problem 2.1 (Weak formulation of the Gross-Pitaevskii equation).

Find such that a.e. in , , and

It is well-known (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 two-grid 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, non-empty 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


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 low-dimensional subspace of . For this purpose, we introduce a two-grid upscaling discretization that was initially proposed in [34] for the treatment of multiscale problems. The framework has been applied to non-linear 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 co-existence 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ément-type interpolation operator (c.f. [15])


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 low-dimensional coarse space (with favorable approximation properties) and a high-dimensional residual space . The residual or ’fine’ space is the kernel of the interpolation operator restricted to ,


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


A basis of is given by with . With this definition we obtain the splitting


Some favorable properties of the decomposition, in particular its -quasi-orthogonality, are discussed in Section 6.2. The minimization problem in the low-dimensional space reads as follows.

Problem 3.3 (Pre-processed 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 Poisson-type 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 pre-processing 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 pre-processed space .

  • Once the coarse space has been assembled it can also be re-used 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


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


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


This yields a modified coarse space with a local basis


The number of non-zero entries of the corresponding finite element matrices is proportional to (note that we expect non-zero 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 Post-processing

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 post-processing step on the fine grid. The post-processing applies the two-grid 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 (Post-processed approximation).

Find with

for all . Define .

Let us emphasize that this approach is different from [18], where the post-precessing problem has a different structure and where classical finite element spaces are used on both scales.

4 A-priori error estimates

This section presents the a-priori error estimates for the pre-processed/upscaled approximation with and without the post-processing 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 post-processed 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 pre-processed approximation).

Assume that . For and as above, it holds


For sufficiently small (in the sense of Cancès et al. [13]), we also have


The proof is postponed to Section 6.3. ∎

The additional post-processing improves, roughly speaking, the order of accuracy by one.

Theorem 4.2 (Error estimates for the post-processed approximation).

Assume that is sufficiently small. The post-processed approximation and the post-processed eigenvalue satisfy:


The constant behaves roughly like and can be extracted from the proofs in Section 6.4.2.


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 Lipschitz-boundary, , and sufficiently small , the fine scale error satisfies the optimal estimate


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 Lipschitz-boundary and . Under these assumptions our a priori estimates for the post-processed 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 Hartree-Fock equations, since it suits our pre-processing framework. The ODA involves solving a linear eigenvalue problem in each iteration step. However, after pre-processing 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 Hartree-Fock 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

Figure 1: Results for harmonic potential. Left: Errors of pre-processed approximation (), (), and () vs. coarse mesh size . Right: Errors of post-processed approximation (), (), and () vs. coarse mesh size .

In the first experiment, we consider uniform coarse meshes with mesh width parameters of . The fine mesh for the pre- and post-processing 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 post-processed 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 super-close in the sense of


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.14.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 two-level 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 speed-up 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 high-resolution numerical approximation

Figure 2: Results for harmonic potential. Left: (Estimated) errors of pre-processed approximation for fixed values (), () and () vs. fine mesh size . Right: (Estimated) errors of post-processed approximation for fixed values (), () and () vs. fine mesh size . In both plots, the (estimated) error of the standard FEM on the fine mesh () is depicted for reference.

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 post-processing computations. The error between the exact eigenvalue and coarse approximations and is estimated via a high-resolution 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 post-processing.

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 self-trapping 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.

Figure 3: Results for periodic potential. Left: Errors of pre-processed approximation (), (), and () vs. coarse mesh size . Right: Errors of post-processed approximation (), (), and () vs. coarse mesh size .

For the discontinuous potential, the experimental rates (with respect to ) are slightly worse than those ones observed in Section 5.1.1. However, they are still better than the rates indicated by Theorems 4.14.2 and conjectured in (15).

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


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


6.2 Properties of the coarse space

Recall the local approximation properties of the weighted Clément-type interpolation operator defined in (4),


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


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.,


The decomposition of in and is -orthogonal


and -quasi-orthogonal in the sense that


The proof is verbatim the same as in [33]. ∎

The following lemma estimates the error of the best-approximation 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 finite-dimensional 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


Given , define (since ) and let denote the corresponding finite element approximation, i.e.,

With , Galerkin-orthogonality 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 energy-norm with an accuracy of order O.

Lemma 6.5 (Stability and approximability of the reference solution).

Let solve Problem 3.2. Then it holds


Recall . Since is a projection, we have

The -orthogonality of (5.c) further yields