In quantum physics, Bose–Einstein condensates (BEC) are important objects of study that feature various interesting properties including macroscopic quantum effects, superfluidity, and occurrence of quantum vortices (in the presence of a magnetic field). In order to model the steady states of BEC consisting of a collection of bosonic quantum particles, the time-independent Gross–Pitaevskii equation (GPE) is widely used, see [16, 28, 24]
. It can be derived from the many-body Schrödinger equation with a given interaction potential in the limit of a large number of particles by applying a Hartree–Fock ansatz of a symmetric tensor product of a single-particle function (in contrast to a single determinant for fermions). Indeed, the Hartree–Fock ansatz becomes exact in the (dilute) mean-field limit (i.e. under certain assumptions on the interaction between particles such as radial symmetry, repulsive and short-range), see[41, 40] for some rigorous results. The GPE is a nonlinear eigenvalue problem that represents the Euler-Lagrange equation of the Gross–Pitaevskii energy functional, given by
under the following normalization constraint for the single particle functions :
Here, , , is a bounded, connected, and open set with Lipschitz boundary, is a potential function with almost everywhere, and is a constant. We note that the fourth-order term in (causing the associated eigenvalue to be nonlinear if ) results from the interaction of particles. The global minimizer of under the constraint (2) is called the (normalized) ground state of (1).
State-of-the-art of numerical methods for BEC
by solving the corresponding Euler-Lagrange formulation, i.e. the GPE,
or by direct minimization (under the normalization constraint (2)).
In the context of (i), classical spatial discretization approaches such as finite element methods [14, 32, 45, 44], finite difference schemes [22, 44], Fourier methods [11, 19], or (pseudo-) spectral methods [13, 9, 12] can be applied in a classical way. One possibility to deal with the nonlinearity occurring in the eigenvalue problem is to employ the Roothaan iteration scheme , also referred to as self-consistent field (SCF) iteration procedure. Alternatively, Newton’s method  or adaptations to the inverse power method  can be used. The idea of combining the iterative solution of the nonlinear eigenvalue problem and of mesh refinements in finite element discretizations has been addressed in [21, 18, 32, 22].
Among the class of methods (ii) we point to the imaginary time method [11, 4, 3, 7] which, upon employing an imaginary time transformation, , is based on the observation that the time-dependent GPE relaxes to the ground state as time evolves. In addition, upon interpreting the imaginary time method as a steepest descent method for the energy functional from (1), the closely related gradient flow approaches [8, 10, 46, 35, 9, 25, 39, 44] can be derived. Further numerical solution methods in the scope of (ii) include a recently proposed preconditioned conjugate gradient method , direct energy minimization using symmetric properties  (which simplify the constraint minimization problem), or the combination of gradient flows and Riemannian optimization .
The aim of this work is to provide a numerical approximation procedure for the ground state of the Gross–Pitaevskii functional (1), under the constraint (2), which is based on a simultaneous interplay of gradient flow iterations and local adaptive finite element mesh refinements. This idea follows the recent developments on the (adaptive) iterative linearized Galerkin (ILG) methodology [33, 34, 23, 1, 2, 37], whereby adaptive discretizations and iterative nonlinearity solvers are combined in an intertwined way; we also refer to the closely related works [30, 29, 15, 31].
A key building block of the numerical scheme to be presented in this paper concerns the decision of whether local mesh refinements or gradient flow iterations should be given preference on a given finite element discretization space. This is accomplished by monitoring the energy decay resulting from the gradient flow, and by performing a comparison to the energy loss caused by the latest mesh refinement. Then, depending on which effect is currently dominant, we either undertake another gradient flow step, or the mesh is refined adaptively. We emphasize that this is a very natural approach for the given problem since both the (conforming) finite element method and the recently proposed gradient flow method  are both guaranteed to be energy-decreasing (owing to the variational principle). The proposed numerical method thereby generates a sequence of finite element approximations defined on adaptively refined spaces which provide a corresponding sequence of monotonically decreasing energies.
Proceeding along the energy minimization approach in 
, the adaptive mesh-refinement strategy in this work is based on identifying a subset of elements in the mesh for which a local refinement will potentially provide a significant (local) contribution to the total energy decay. To this end, for each element in the mesh, we first apply a local gradient flow step on a locally refined patch; these computations, since local and independent, can be done in parallel and only involve very few degrees of freedom. Then, a Dörfler marking strategy[27, Sec. 4.2] selects the most promising elements for refinement. The numerical tests illustrate optimal convergence rates in the number of unknowns for a variety of examples—both linear or nonlinear models, with smooth or irregular potentials, will be investigated.
In Section 2 we present the framework of the Gross-Pitaevskii equation and motivate its associated gradient flow system. Moreover, Section 3 is devoted to the finite element discretization and the adaptive mesh refinement strategy. Furthermore, Section 4 presents various numerical tests in 2D. Finally, we add some concluding remarks in Section 5.
2. The Gross–Pitaevskii equation and gradient flows
2.1. Nonlinear eigenvalue formulation
We observe that the energy functional from (1) is Fréchet differentiable on the Sobolev space (defined to be the space of all functions with weak gradient in and zero trace along the boundary ); indeed, a simple calculation reveals that
where signifies the dual product in . The Euler-Lagrange formulation of the constrained minimization problem
with signifying the -unit sphere in , is given by
is an eigenfunction of (4) associated with an eigenvalue , then we note that
Given that (almost everywhere in ) and , the Gross–Pitaevskii eigenvalue problem (4) has a unique (-normalized) positive eigenfunction which is the ground state of the Bose–Einstein condensate (1), see [35, Lem. 5.4]; in particular, is an eigenfunction to the minimal (and simple) eigenvalue, signified by of (4), see .
2.2. Continuous gradient flow
The ground state will be determined iteratively. To this end, we employ the projected gradient flow approach proposed in . One of the key ideas is to introduce a weighted energy inner product on : For fixed , we let
Owing to the Riesz representation theorem, for any , it exists a unique such that
If , we notice that , and therefore . Hence, we may consider the linear mapping
for , we have , and is the identity map. Using (7), it is fairly elementary to verify that
i.e. is the orthogonal projection onto the (tangential plane) with respect to the -inner product.
Based on the above definitions, we are now ready to present the gradient flow induced by the inner product from (6). More precisely, we consider a trajectory which, for a given initial value , follows the dynamical system
The existence of a solution has been discussed in [35, Sec. 3.1 & 3.2].
For any , we define to be the Riesz representative of the steepest descend direction at with respect to the -inner product, i.e.
for any . This implies that the gradient flow from (10) follows the orthogonal projection of the steepest descend direction , for , onto the tangential plane .
In particular, it follows that . This means that the gradient flow stays on the sphere for any , i.e., physically speaking, it is mass preserving; cf. [35, Lem. 3.3].
Since is nonnegative for any , the monotonicity property from (iii) implies that there is with . Hence, applying (12), we obtain the identity
Therefore, exploiting (7), for all , we infer that
i.e. is an eigenfunction for the GPE (4) to the eigenvalue .
2.3. Discrete gradient flow
For the purpose of computing an approximation of the continuous gradient flow trajectory from (10), we use the forward Euler discretization method. For a given initial value this yields a sequence of functions which, for , is defined by
Here, is a sequence of positive (discrete) time steps that is assumed uniformly bounded from above and below with bounds and , respectively, such that . The scheme (13) is called discrete gradient flow iteration (GFI).
Provided that the maximal time step is sufficiently small, it can be seen that the GFI scheme (13) yields guranteed energy reduction. More precisely, there exists such that for all it holds that , see [35, Lem. 4.7]. Moreover, if for all then, for any starting value with , the (full) sequence generated by the GFI (13) satisfies for all and converges strongly in to the unique positive ground state , see [35, Theorem 5.1].
3. Adaptive gradient flow finite element discretization
We now focus on the adaptive spatial discretization of the gradient flow iteration scheme (13).
3.1. Finite element discretization
Consider a sequence of conforming and shape-regular partitions of the domain into simplicial elements (i.e. triangles for and tetrahedra for ). Moreover, for a (fixed) polynomial degree and any subset , we introduce the finite element space
with signifying the (local) space of all polynomials of maximal total degree on , . Furthermore, similarly as before, we denote by
the -unit sphere in . In the sequel, we apply the notations and . We further denote by the restriction of the energy functional from (1) to the Galerkin space . Then, due to the compactness of , it exists a minimizer of , i.e. , with . Furthermore, if is a dense family of finite element subspaces of , then any sequence of minimizers of with converges in to the ground state of ; we refer to [19, Theorem 1] or [47, Theorem 3.1].
3.2. Discrete GFI
Let us define the (space) discrete version of the gradient flow iteration (13) on a finite element subspace . For we denote by the unique solution of
cf. (7). Then, for , the space discrete GFI is given by
with a sequence of discrete time steps as in (13).
Consider a fixed mesh and associated approximation space . Let be the sequence generated by the discrete GFI (14) with some initial value . If , with as in Remark 2.1, for all , then the corresponding energies are strictly monotone decreasing, and there exists a limit energy . Furthermore, up to subsequences, we have strongly in , where , with , is a discrete eigenfunction of the corresponding GPE, i.e. there is so that
We refer to [35, Corollary 4.11] for details.
In practical computations, in order to guarantee a positive energy decay in each iteration step, we propose the time step strategy within (14) given by
where, for , we write to denote the output of the discrete GFI (14) based on the time step and on the previous approximation . We observed in several examples that for the choice , i.e. using above, no time correction was needed; we also refer to [35, Remark 4.8] for a discussion of the fixed time step . For that reason, and for the sake of keeping the computational cost minimal, we fix the time step in the local GFI from Algorithm 1 below. We still use, however, the time step strategy for the global GFI in Algorithm 2.
3.3. Local energy decay and adaptive mesh refinements
For any element we consider the open patch comprising of and its immediate face-wise neighbours. Moreover, given , we define the modified patch by uniformly (red) refining the element into a (fixed) number of subelements; here, we assume that the introduction of any hanging nodes in is removed by doing (e.g. green) refinements, see Figure 1.
We consider basis functions of the locally supported space . Furthermore, for any given , we introduce the extended space
Suppose we have found an accurate approximation of the discrete GPE (15), for some . Then, by performing one local discrete GFI-step in we obtain a new local approximation, denoted by , with . We emphasize that has a small dimension, and hence the discrete GFI (14) based on entails hardly any computational cost (for instance, for dimension and polynomial degree , the dimension of the locally refined space , cf. Figure 1, is typically 3 or 4).
for all . The value indicates the potential energy reduction due to a refinement of the element . This observation motivates the energy-based adaptive mesh refinement procedure outlined in Algorithm 1.
|Perform one discrete GFI-step in the low-dimensional space to obtain a potentially improved local approximation .|
3.4. Adaptive strategy
From a practical viewpoint, once the discrete GFI approximation from (14) is close to a solution of (15), on a given finite element space, we expect that any further GFI steps will no longer contribute an essential decay to the energy in (17). In this case, in order to further reduce the energy, we need to enrich the finite element space appropriately. More specifically, for , suppose that we have performed a reasonable number (possibly depending on ) of GFI-iterations (14) in . Consider now a (hierarchically) refined mesh of , for example, obtained by Algorithm 1. Then we may embed the final guess on the previous space into the enriched finite element space in order to obtain an initial guess on the refined mesh :
For each GFI-iteration we monitor two quantities. Firstly, we introduce the increment on each iteration given by
Secondly, we compare to the energy loss as compared to the previous mesh refinement, i.e.
We stop the iteration for as soon as becomes small compared to , i.e. once there is no notable benefit (relatively speaking) in performing any more discrete GFI steps on the current space . Specifically, for , this is expressed by the bound
for some parameter . We implement this procedure in Algorithm 2.
|Mark and adaptively refine the mesh using Algorithm 1 to generate a new mesh .|
4. Numerical Experiments
We apply Algorithm 2 for some numerical computations in two space dimensions, i.e. , with Cartesian coordinates denoted by . In all examples, we choose the initial guess such that for any node in the interior of the initial mesh , where is the appropriate constant to fulfil the norm constraint (2). Moreover, we set , as well as and in Algorithms 1 and 2, respectively. Even if the stopping criterion in Line 12 of Algorithm 2 may not be satisfied, for the purpose of our tests, we stop the computations once the number of degrees of freedom (i.e. the dimension of the finite element space ) exceeds .
4.1. Linear GPE with smooth potential
We consider first the case where in (1), and is a smooth function:
Note that the associated eigenvalue problem (4) is linear for this example. It is known that, for , the energy of the ground state is given by , see, e.g., . We note that the mass of is essentially concentrated at the origin due to the global minimum of at . Therefore, the restriction to the bounded domain , which we use in our computations, has almost no effect on the minimal value of the ground state energy. Figure 2 (left) illustrates that the approximation of the energy converges with optimal rate in terms of the numbers of degrees of freedom. In addition, we see that the mesh is mainly refined around the origin, where the mass of the ground state is concentrated, see Figure 2 (right). These results underline that the proposed adaptive gradient flow procedure effectively detects the local behaviour of the model.
4.2. Linear GPE with singular potential .
We perform another experiment with in (1), i.e. the associated eigenvalue problem (4) is again linear. In contrast to the previous example, however, we consider a potential which features a severe point singularity at the origin . Specifically, the energy functional is given by
with . In order to properly resolve the singularity, we use a higher-order quadrature formula, and begin with an initial mesh that exhibits some a priori refinement at the origin, see Figure 3 (left). This experiment has been conducted already in , where the authors obtained an approximated minimal eigenvalue for the corresponding GPE (4); we compare our results to this reference ground state energy (divided by 2 on account of (5)). We can see from Figure 3 (right) that the proposed Algorithm 2 achieves a sequence of energy approximations for the ground state which decays at an (almost) optimal rate.
4.3. Linear GPE with potential wells
We test a final example with . The potential is given by the sum of four Gaussian bells, see Figure 4 (left). This experiment is borrowed from [42, Experiment 4.2], however, with a (constant) shift such that in the underlying domain . As we can see from Figure 4 (right), the mass of the ground state is mainly concentrated at two adjacent hills in the lower left part of the domain; we note that this is perfectly in line with the results obtained in the paper . Moreover, the energy-based adaptive mesh refinement has properly resolved the two local hills featured in the ground state , see Figure 5.
The authors from  have computed an approximated minimal eigenvalue , whereby this is the value adapted for our reformulation of the problem. Based on this approximation, we observe an optimal rate of convergence for the energy of the ground state in Figure 6 (left). We remark, for this example, that the performance of Algorithm 2 crucially depends on the choice of the parameter in Line 6. Indeed, if (instead of ) is selected, then the numerical results exhibit a considerably less favourable asymptotic convergence regime, see Figure 6 (right). An explanation can be inferred from Figure 7: We see that the choice leads to a significantly higher number of gradient flow steps on each Galerkin space , which seems to be essential for the effective numerical solution of this problem.
4.4. Nonlinear GPE with harmonic confinement potential
We now consider a nonlinear Bose-Einstein condensate, i.e. in (1), and use the smooth potential :
the domain is given by . This experiment was conducted previously in , where an approximation for the energy of the ground state has been documented. Based on this approximation, we test our Algorithm 2. Again, we obtain (almost) optimal convergence for the approximation of the energy of the ground state, see Figure 8.
4.5. Nonlinear GPE with optical lattice potential
As before, we choose and , with an oscillating potential function , see Figure 10 (left) for its contour plot. More precisely, the energy functional is given by
This experiment was also considered in  with an asserted approximation of the ground state energy. Based on the adaptive Algorithm 2 presented in this work, a smaller value for the ground state energy has been computed; we suppose that this (improved) approximation results from the adaptive (and thereby more effective) refinement of the meshes. To be specific, we have obtained the approximation based on an adaptively refined mesh with degrees of freedom; here, the underlined digits are stable (i.e. the computations indicate that they do not change anymore as the iterations continue). In Figure 9 (left) we have depicted the error for the approximations of the ground state energy with respect to our reference value. This plot indicates an asymptotically optimal rate of convergence of Algorithm 2 for the given problem. The suboptimal behaviour in the pre-asymptotic phase is due to the fact that the mesh is not fine enough to resolve the oscillations of the potential sufficiently well.
As in Example 4.3, we have also run this experiment for both values and . In contrast to the previous test, we observe no considerable difference in the performance of the corresponding computations. Indeed, except for the first few mesh refinements, the number of gradient flow steps on each Galerkin space is the same for both cases, see Figure 11.
4.6. Nonlinear energy functional with a nonsymmetric potential
Finally, we revisit the test problem [10, Example 4.3.II]:
with the symmetric domain . In  an approximated value of for the energy of the ground state has been obtained. Again, based on an adaptively refined mesh with degrees of freedom, we have obtained an approximation . Taking this value as a reference ground state energy for the given example, we observe optimal convergence, see Figure 12.
In this work, we have considered a computational procedure for the numerical approximation of the ground state and its associated energy of the Gross-Pitaevskii equation, which applies an effective interplay of a gradient flow iteration method and adaptive mesh refinements. Both of these techniques rely on energy minimization and guaranteed energy reduction. Thereby, they are based on the underlying structure of the problem at hand in a very natural way. Our scheme is fairly simple to implement and, for the test problems presented here, exhibits either optimal or close to optimal convergence rates for the approximation of the ground state energy.
-  M. Amrein and T. P. Wihler, An adaptive Newton-method based on a dynamical systems approach, Commun. Nonlinear Sci. Numer. Simul. 19 (2014), no. 9, 2958–2973.
by same author,
Fully adaptive Newton-Galerkin methods for semilinear elliptic partial differential equations, SIAM J. Sci. Comput. 37 (2015), no. 4, A1637–A1657.
-  X. Antoine and R. Duboscq, Gpelab, a matlab toolbox to solve gross–pitaevskii equations i: Computation of stationary solutions, Computer Physics Communications 185 (2014), no. 11, 2969–2991.
-  by same author, Robust and efficient preconditioned krylov spectral solvers for computing the ground states of fast rotating and strongly interacting bose–einstein condensates, Journal of Computational Physics 258 (2014), 509–523.
-  X. Antoine, A. Levitt, and Q. Tang, Efficient spectral computation of the stationary states of rotating bose–einstein condensates by preconditioned nonlinear conjugate gradient methods, Journal of Computational Physics 343 (2017), 92–109.
-  W. Bao, Mathematical models and numerical methods for Bose-Einstein condensation, Proceedings of the International Congress of Mathematicians—Seoul 2014. Vol. IV, Kyung Moon Sa, Seoul, 2014, pp. 971–996.
-  W. Bao and Y. Cai, Mathematical theory and numerical methods for bose-einstein condensation, Tech. Report 1212.5341, arxiv.org, 2012.
-  W. Bao, Y. Cai, and H. Wang, Efficient numerical methods for computing ground states and dynamics of dipolar bose–einstein condensates, Journal of Computational Physics 229 (2010), no. 20, 7874–7892.
-  W. Bao, I-L. Chern, and F.Y. Lim, Efficient and spectrally accurate numerical methods for computing ground and first excited states in bose–einstein condensates, Journal of Computational Physics 219 (2006), no. 2, 836–854.
-  W. Bao and Q. Du, Computing the ground state solution of bose–einstein condensates by a normalized gradient flow, SIAM Journal on Scientific Computing 25 (2004), no. 5, 1674–1697.
-  W. Bao, D. Jaksch, and P. Markowich, Numerical solution of the gross–pitaevskii equation for bose–einstein condensation, Journal of Computational Physics 187 (2003), no. 1, 318–342.
-  W. Bao, H. Li, and J. Shen, A generalized-laguerre–fourier–hermite pseudospectral method for computing the dynamics of rotating bose–einstein condensates, SIAM Journal on Scientific Computing 31 (2009), no. 5, 3685–3711.
-  W. Bao and J. Shen, A fourth-order time-splitting laguerre–hermite pseudospectral method for bose–einstein condensates, SIAM Journal on Scientific Computing 26 (2005), no. 6, 2010–2028.
-  W. Bao and W. Tang, Ground-state solution of bose–einstein condensate by directly minimizing the energy functional, Journal of Computational Physics 187 (2003), no. 1, 230–254.
-  C. Bernardi, J. Dakroub, G. Mansour, and T. Sayah, A posteriori analysis of iterative algorithms for a nonlinear problem, J. Sci. Comput. 65 (2015), no. 2, 672–697.
-  S.N. Bose, Plancks gesetz und lichtquantenhypothese, Zeitschrift für Physik 26 (1924), no. 1, 178–181.
-  M. Caliari, A. Ostermann, S. Rainer, and M. Thalhammer, A minimisation approach for computing the ground state of gross–pitaevskii systems, Journal of Computational Physics 228 (2009), no. 2, 349–360.
-  E. Cancès, R. Chakir, L. He, and Y. Maday, Two-grid methods for a class of nonlinear elliptic eigenvalue problems, IMA Journal of Numerical Analysis 38 (2017), no. 2, 605–645.
-  E. Cancès, R. Chakir, and Y. Maday, Numerical analysis of nonlinear eigenvalue problems, J. Sci. Comput. 45 (2010), no. 1-3, 90–117.
-  E. Cancès, R. Chakir, and Y. Maday, Numerical analysis of nonlinear eigenvalue problems, J. Sci. Comput. 45 (2010), 90–117.
E. Cancès, G. Dusson, Y. Maday, B. Stamm, and M. Vohralík,
A perturbation-method-based a posteriori estimator for the planewave discretization of nonlinear schrödinger equations, Comptes Rendus Mathematique 352 (2014), no. 11, 941–946.
-  C-S. Chien, H-T. Huang, B-W. Jeng, and Z-C. Li, Two-grid discretization schemes for nonlinear schrödinger equations, Journal of Computational and Applied Mathematics 214 (2008), no. 2, 549–571.
-  S. Congreve and T. P. Wihler, Iterative Galerkin discretizations for strongly monotone problems, Journal of Computational and Applied Mathematics 311 (2017), 457–472.
-  F. Dalfovo, S. Giorgini, L.P. Pitaevskii, and S. Stringari, Theory of bose-einstein condensation in trapped gases, Reviews of Modern Physics 71 (1999), no. 3, 463.
-  I. Danaila and P. Kazemi, A new sobolev gradient method for direct minimization of the gross–pitaevskii energy with rotation, SIAM Journal on Scientific Computing 32 (2010), no. 5, 2447–2467.
-  I. Danaila and B. Protas, Computation of ground states of the gross–pitaevskii functional via riemannian optimization, SIAM Journal on Scientific Computing 39 (2017), no. 6, B1102–B1129.
-  W. Dörfler, A convergent adaptive algorithm for Poisson’s equation, SINUM 33 (1996), 1106–1124.
-  A. Einstein, Quantentheorie des einatomigen idealen gases, SB Preuss. Akad. Wiss. phys.-math. Klasse (1924), 261–267.
-  L. El Alaoui, A. Ern, and M. Vohralík, Guaranteed and robust a posteriori error estimates and balancing discretization and linearization errors for monotone nonlinear problems, Comput. Methods Appl. Mech. Engrg. 200 (2011), no. 37-40, 2782–2795.
-  A. Ern and M. Vohralík, Adaptive inexact Newton methods with a posteriori stopping criteria for nonlinear diffusion PDEs, SIAM J. Sci. Comput. 35 (2013), no. 4, A1761–A1791.
-  E. M. Garau, P. Morin, and C. Zuppa, Convergence of an adaptive Kačanov FEM for quasi-linear problems, Appl. Numer. Math. 61 (2011), no. 4, 512–529.
-  X-G. Gong, L. Shen, D. Zhang, and A. Zhou, Finite element approximations for schrödinger equations with applications to electronic structure computations, Journal of Computational Mathematics (2008), 310–323.
-  P. Heid and T.P. Wihler, Adaptive iterative linearization Galerkin methods for nonlinear problems, Tech. Report 1808.04990, arxiv.org, 2018.
-  by same author, On the convergence of adaptive iterative linearized Galerkin methods, Tech. Report 1905.06682, arxiv.org, 2019.
-  P. Henning and D. Peterseim, Sobolev gradient flow for the gross-pitaevskii eigenvalue problem: global convergence and computational efficiency, Tech. Report 1812.00835, arxiv.org, 2018.
-  P. Houston and T. P. Wihler, Adaptive energy minimisation for -finite element methods, Comput. Math. Appl. 71 (2016), no. 4, 977 – 990.
-  by same author, An -adaptive newton-discontinuous-galerkin finite element approach for semilinear elliptic boundary value problems, Math. Comp. 87 (2018), no. 314, 2641–2674.
E. Jarlebring, S. Kvaal, and W. Michiels,
An inverse iteration method for eigenvalue problems with eigenvector nonlinearities, SIAM Journal on Scientific Computing 36 (2014), no. 4, A1978–A2001.
-  P. Kazemi and M. Eckart, Minimizing the gross-pitaevskii energy functional with the sobolev gradient - analytical and numerical results, International Journal of Computational Methods 7 (2010), no. 03, 453–475.
-  M. Lewin, Mean-field limit of bose systems: rigorous results, Tech. Report 1510.04407, arxiv.org, 2015.
-  E.H. Lieb, R. Seiringer, and J. Yngvason, Bosons in a trap: A rigorous derivation of the gross-pitaevskii energy functional, The Stability of Matter: From Atoms to Stars, Springer, 2001, pp. 685–697.
-  L. Lin and B. Stamm, A posteriori error estimates for discontinuous galerkin methods using non-polynomial basis functions. part ii: Eigenvalue problems, ESAIM: M2AN 51 (2017), no. 5, 1733–1753.
-  Y. Maday and C. Marcati, Regularity and discontinuous galerkin finite element approximation of linear elliptic eigenvalue problems with singular potentials, Tech. Report 1810.09010, arxiv.org, 2018.
-  N. Raza, S. Sial, S.S. Siddiqi, and T. Lookman, Energy minimization related to the nonlinear schrödinger equation, Journal of Computational Physics 228 (2009), no. 7, 2572–2577.
-  H. Xie and M. Xie, A multigrid method for ground state solution of bose-einstein condensates, Communications in Computational Physics 19 (2016), no. 3, 648–662.
-  R. Zeng and Y. Zhang, Efficiently computing vortex lattices in rapid rotating bose–einstein condensates, Computer Physics Communications 180 (2009), no. 6, 854–860.
-  A. Zhou, An analysis of finite-dimensional approximations for the ground state solution of Bose-Einstein condensates, Nonlinearity 17 (2004), no. 2, 541–550.