Gradient Flow Finite Element Discretizations with Energy-Based Adaptivity for the Gross-Pitaevskii Equation

TL;DR

This paper introduces an energy-based adaptive finite element method for solving the steady-state Gross-Pitaevskii equation, combining gradient flow iterations with mesh refinement to achieve high accuracy and optimal convergence.

Contribution

The paper develops a novel adaptive finite element approach based on energy minimization and gradient flow for the Gross-Pitaevskii equation, demonstrating improved accuracy and convergence.

Findings

Achieves highly accurate solutions for the Gross-Pitaevskii equation.

Demonstrates optimal convergence rates with respect to degrees of freedom.

Effective combination of gradient flow and adaptive mesh refinement.

Abstract

We present an effective adaptive procedure for the numerical approximation of the steady-state Gross-Pitaevskii equation. Our approach is solely based on energy minimization, and consists of a combination of gradient flow iterations and adaptive finite element mesh refinements. Numerical tests show that this strategy is able to provide highly accurate results, with optimal convergence rates with respect to the number of freedom.

Equations

$$\mathsf{E}(v):={\int }_{\Omega }\left(\frac{1}{2}|\nabla v{|}^2+V(\mathbf{x})|v{|}^2+\frac{\beta }{2}|v{|}^4\right)\mathsf{d}\mathbf{x},$$

$${\Vert v\Vert }_{{\mathrm{L}}^2(\Omega )}=1.$$

$$⟨{\mathsf{E}}^{\prime }(v),w⟩={\int }_{\Omega }\left(\nabla v\cdot \nabla w+2V(\mathbf{x})vw+2\beta |v{|}^{2}vw\right)\mathsf{d}\mathbf{x},v,w\in \mathbb{H},$$

$$\underset{v\in S_{\mathbb{H}}}{\mathrm{argmin}}\mathsf{E}(v),$$

$$v\in \mathbb{H}:⟨{\mathsf{E}}^{\prime }(v),w⟩=\lambda (v,w{)}_{{\mathrm{L}}^2(\Omega )}\forall w\in \mathbb{H}.$$

$$\lambda =⟨{\mathsf{E}}^{\prime }(v),v⟩=2\mathsf{E}(v)+\beta \Vert v{\Vert }_{{\mathrm{L}}^4(\Omega )}^4>0.$$

$$a_z(v,w):={\int }_{\Omega }\left(\nabla v\cdot \nabla w+2V(\mathbf{x})vw+2\beta |z{|}^{2}vw\right)\mathsf{d}\mathbf{x},v,w\in \mathbb{H}.$$

$$a_z({\mathsf{G}}_z(v),w)=(v,w{)}_{{\mathrm{L}}^2(\Omega )}\forall w\in \mathbb{H}.$$

$${\mathsf{P}}_z:\mathbb{H}\to {\mathbb{T}}_z,{\mathsf{P}}_z(v)=v-\frac{(z,v{)}_{{\mathrm{L}}^2(\Omega )}}{a_z({\mathsf{G}}_z(z),{\mathsf{G}}_z(z))}{\mathsf{G}}_z(z),$$

$${\mathbb{T}}_z:=\{w\in \mathbb{H}:(z,w{)}_{{\mathrm{L}}^2(\Omega )}=0\};$$

$$a_z(v-{\mathsf{P}}_z(v),w)=0\forall v\in \mathbb{H},\forall w\in {\mathbb{T}}_z,$$

$$\dot{u}(t)=-{\mathsf{P}}_{u(t)}(u(t)),t>0.$$

$$a_z({\nabla }_z\mathsf{E}(v),w)=⟨{\mathsf{E}}^{\prime }(v),w⟩\forall w\in \mathbb{H}.$$

$$a_{u(t)}(\dot{u}(t),w)$$

$$=a_{u(t)}(-u(t),{\mathsf{P}}_{u(t)}(w))$$

$$=⟨-{\mathsf{E}}^{\prime }(u(t)),{\mathsf{P}}_{u(t)}(w)⟩$$

$$=a_{u(t)}(-{\nabla }_{u(t)}\mathsf{E}(u(t)),{\mathsf{P}}_{u(t)}(w))$$

$$=a_{u(t)}({\mathsf{P}}_{u(t)}(-{\nabla }_{u(t)}\mathsf{E}(u(t))),w),$$

$$\dot{u}(t)\in {\mathbb{T}}_{u(t)}\forall t>0.$$

$$\frac{\mathsf{d}}{\mathsf{dt}}{\Vert u(t)\Vert }_{{\mathrm{L}}^2(\Omega )}^2=2(u(t),\dot{u}(t){)}_{{\mathrm{L}}^2(\Omega )}=0\forall t>0.$$

$$\frac{\mathsf{d}}{\mathsf{dt}}\mathsf{E}(u(t))=⟨{\mathsf{E}}^{\prime }(u(t)),\dot{u}(t)⟩=a_{u(t)}(u(t),\dot{u}(t))=-a_{u(t)}(-{\mathsf{P}}_{u(t)}(u(t)),\dot{u}(t)).$$

$$⟨{\mathsf{E}}^{\prime }(u(t)),\dot{u}(t)⟩=-a_{u(t)}(\dot{u}(t),\dot{u}(t))\le 0\forall t>0,$$

$$0\le {\int }_0^{\infty }{a}_{u(\tau )}(\dot{u}(\tau ),\dot{u}(\tau ))\mathsf{d}\tau =\mathsf{E}(u_0)-{\mathsf{E}}^{\star }<\infty ,$$

$${\lambda }^{\star }=\frac{\Vert u^{\star }{\Vert }_{{\mathrm{L}}^2(\Omega )}^2}{a_{u^{\star }}({\mathsf{G}}_{u^{\star }}(u^{\star }),{\mathsf{G}}_{u^{\star }}(u^{\star }))}.$$

$$⟨{\mathsf{E}}^{\prime }(u^{\star }),w⟩=a_{u^{\star }}(u^{\star },w)={\lambda }^{\star }{a}_{u^{\star }}({\mathsf{G}}_{u^{\star }}(u^{\star }),w)={\lambda }^{\star }(u^{\star },w{)}_{{\mathrm{L}}^2(\Omega )},$$

$$u^{n+1}$$

$${\hat{u}}^{n+1}$$

$$\mathbb{V}(\omega )=\left\{v\in \mathbb{H}:v{|}_{\kappa }\in {\mathbb{P}}_p(\kappa ),\kappa \in \omega ,v{|}_{\Omega \setminus \omega }=0\right\},$$

$$S_{\mathbb{V}(\omega )}=\left\{v\in \mathbb{V}(\omega ):\Vert v{\Vert }_{{\mathrm{L}}^2(\Omega )}=1\right\}$$

$$a_u({\mathsf{G}}^N(u),v)=(u,v{)}_{{\mathrm{L}}^2(\Omega )}\forall v\in {\mathbb{X}}_N;$$

$$u_N^{n+1}$$

$${\hat{u}}_N^{n+1}$$

$$a_{u_N^{\star }}(u_N^{\star },v)={\lambda }_N^{\star }(u_N^{\star },v{)}_{{\mathrm{L}}^2(\Omega )}\forall v\in {\mathbb{X}}_N.$$

Full text

The authors acknowledge the financial support of the Swiss National Science Foundation (SNF), Grant No. 200021 182524

1. Introduction

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 [17, 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 [43, 42] 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

[TABLE]

under the following normalization constraint for the single particle functions $v$:

[TABLE]

Here, $\Omega\subset\mathbb{R}^{d}$, $d=\{1,2,3\}$, is a bounded, connected, and open set with Lipschitz boundary, $V\in L^{\infty}(\Omega)$ is a potential function with $V\geq 0$ almost everywhere, and $\beta\geq 0$ is a constant. We note that the fourth-order term in $\mathsf{E}$ (causing the associated eigenvalue to be nonlinear if $\beta>0$) results from the interaction of particles. The global minimizer of $\mathsf{E}$ under the constraint (2) is called the (normalized) ground state of (1).

State-of-the-art of numerical methods for BEC

An overview of models and numerical schemes for the GPE is provided in [7]. In general terms, the minimizer of the Gross–Pitaevskii energy functional (1) can be found either

1. (i)

by solving the corresponding Euler-Lagrange formulation, i.e. the GPE, 2. (ii)

or by direct minimization (under the normalization constraint (2)).

In the context of (i), classical spatial discretization approaches such as finite element methods [15, 33, 47, 46], finite difference schemes [22, 46], Fourier methods [12, 20], or (pseudo-) spectral methods [14, 10, 13] 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 [21], also referred to as self-consistent field (SCF) iteration procedure. Alternatively, Newton’s method [18] or adaptations to the inverse power method [39] 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, 19, 33, 22].

Among the class of methods (ii) we point to the imaginary time method [12, 4, 3, 8] which, upon employing an imaginary time transformation, $t\to-\imath t$, 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 $\mathsf{E}$ from (1), the closely related gradient flow approaches [9, 11, 48, 36, 10, 25, 40, 46] can be derived. Further numerical solution methods in the scope of (ii) include a recently proposed preconditioned conjugate gradient method [5], direct energy minimization using symmetric properties [15] (which simplify the constraint minimization problem), or the combination of gradient flows and Riemannian optimization [26].

Contribution

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 [34, 35, 23, 1, 2, 38], whereby adaptive discretizations and iterative nonlinearity solvers are combined in an intertwined way; we also refer to the closely related works [30, 29, 16, 32].

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 [36] 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 [37], 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.

Outline

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 $\mathsf{E}$ from (1) is Fréchet differentiable on the Sobolev space $\mathbb{H}:=\mathrm{H}^{1}_{0}(\Omega)$ (defined to be the space of all functions with weak gradient in $\mathrm{L}^{2}(\Omega)^{d}$ and zero trace along the boundary $\partial\Omega$); indeed, a simple calculation reveals that

[TABLE]

where $\langle\cdot,\cdot\rangle$ signifies the dual product in $\mathbb{H}^{\star}\times\mathbb{H}$. The Euler-Lagrange formulation of the constrained minimization problem

[TABLE]

with $S_{\mathbb{H}}:=\{v\in\mathbb{H}:\|v\|_{{\rm L}^{2}(\Omega)}=1\}$ signifying the $\mathrm{L}^{2}(\Omega)$-unit sphere in $\mathbb{H}$, is given by

[TABLE]

Here, the scalar $\lambda$ takes the role of a Lagrange multiplier corresponding to the norm constraint (2). The nonlinear eigenvalue problem (4) is called the Gross–Pitaevskii equation (GPE). If $v\in S_{\mathbb{H}}$ is an eigenfunction of (4) associated with an eigenvalue $\lambda$, then we note that

[TABLE]

Given that $V\geq 0$ (almost everywhere in $\Omega$) and $\beta\geq 0$, the Gross–Pitaevskii eigenvalue problem (4) has a unique ($\mathrm{L}^{2}(\Omega)$-normalized) positive eigenfunction $u_{\mathrm{GS}}>0$ which is the ground state of the Bose–Einstein condensate (1), see [36, Lem. 5.4]; in particular, $u_{\mathrm{GS}}$ is an eigenfunction to the minimal (and simple) eigenvalue, signified by $\lambda_{\mathrm{GS}}$ of (4), see [20].

2.2. Continuous gradient flow

The ground state $u_{\mathrm{GS}}$ will be determined iteratively. To this end, we employ the projected gradient flow approach proposed in [36]. One of the key ideas is to introduce a weighted energy inner product on $\mathbb{H}\times\mathbb{H}$: For fixed $z\in\mathbb{H}$, we let

[TABLE]

Owing to the Riesz representation theorem, for any $v\in\mathbb{H}$, it exists a unique $\mathsf{G}_{z}(v)\in\mathbb{H}$ such that

[TABLE]

If $z\neq 0$, we notice that $a_{z}(\mathsf{G}_{z}(z),z)=\|z\|^{2}_{\mathrm{L}^{2}(\Omega)}>0$, and therefore $\mathsf{G}_{z}(z)\neq 0$. Hence, we may consider the linear mapping

[TABLE]

where

[TABLE]

for $z=0$, we have $\mathbb{T}_{0}=\mathbb{H}$, and $\mathsf{P}_{0}$ is the identity map. Using (7), it is fairly elementary to verify that

[TABLE]

i.e. $\mathsf{P}_{z}$ is the orthogonal projection onto the (tangential plane) $\mathbb{T}_{z}$ with respect to the $a_{z}$-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 $u:\,[0,\infty)\to\mathbb{H}$ which, for a given initial value $u(0)=u_{0}\in S_{\mathbb{H}}$, follows the dynamical system

[TABLE]

The existence of a solution $u$ has been discussed in [36, Sec. 3.1 & 3.2].

In order to motivate the gradient flow system (10), we briefly revisit some arguments from [36, Proof of Theorem 3.2]:

1. (i)

For any $v,z\in\mathbb{H}$, we define $\nabla_{z}\mathsf{E}(v)\in\mathbb{H}$ to be the Riesz representative of the steepest descend direction $\mathsf{E}^{\prime}(v)$ at $v$ with respect to the $a_{z}$-inner product, i.e.

[TABLE]

Then, using (10), employing (9) respectively (6), and recalling (3), for any $t>0$, we notice that

[TABLE]

for any $w\in\mathbb{H}$. This implies that the gradient flow from (10) follows the orthogonal projection of the steepest descend direction $-\nabla_{u(t)}\mathsf{E}(u(t))$, for $t>0$, onto the tangential plane $\mathbb{T}_{u(t)}$. 2. (ii)

Observe that

[TABLE]

Hence,

[TABLE]

In particular, it follows that $\left\|u(t)\right\|_{\mathrm{L}^{2}(\Omega)}=\left\|u(0)\right\|_{\mathrm{L}^{2}(\Omega)}=1$. This means that the gradient flow stays on the sphere $S_{\mathbb{H}}$ for any $t>0$, i.e., physically speaking, it is mass preserving; cf. [36, Lem. 3.3]. 3. (iii)

Using again (11) and invoking (9), for any $t>0$, it holds that

[TABLE]

Thus, implementing (10) yields

[TABLE]

i.e. the energy $\mathsf{E}(t)$ is monotone decreasing as $t\to\infty$; cf. [36, Lem. 3.3]. 4. (iv)

Since $\mathsf{E}(v)$ is nonnegative for any $v\in\mathbb{H}$, the monotonicity property from (iii) implies that there is $\mathsf{E}^{\star}\geq 0$ with $\lim_{t\to\infty}\mathsf{E}(u(t))=\mathsf{E}^{\star}<\infty$. Hence, applying (12), we obtain the identity

[TABLE]

which implies that $\int_{0}^{\infty}\|\nabla\dot{u}(\tau)\|^{2}_{\mathrm{L}^{2}(\Omega)}\,\mathsf{d}\tau$ is bounded. It follows that $u(t)$ has a limit $u^{\star}\in\mathbb{H}$ as $t\to\infty$ (with $\left\|u^{\star}\right\|_{{\rm L}^{2}(\Omega)}=1$, cf. (ii)), and $\mathsf{P}_{u^{\star}}(u^{\star})=0$ upon taking the limit in (10). Consequently, we deduce from (8) that $u^{\star}=\lambda^{\star}\mathsf{G}_{u^{\star}}(u^{\star})$, where

[TABLE]

Therefore, exploiting (7), for all $w\in\mathbb{H}$, we infer that

[TABLE]

i.e. $u^{\star}\not\equiv 0$ is an eigenfunction for the GPE (4) to the eigenvalue $\lambda^{\star}$. We note that $u^{\star}$ can be any eigenfunction satisfying (13), respectively (4); in particular, it is not necessarily the ground state.

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 $u^{0}\in S_{\mathbb{H}}$ this yields a sequence of functions $\{u^{n}\}_{n\geq 0}\subset S_{\mathbb{H}}$ which, for $n\geq 0$, is defined by

[TABLE]

Here, $\{\tau_{n}\}_{n\geq 0}$ is a sequence of positive (discrete) time steps that is assumed uniformly bounded from above and below with bounds $\tau_{\max}$ and $\tau_{\min}$, respectively, such that $0<\tau_{\min}\leq\tau_{n}\leq\tau_{\max}<\infty$. The scheme (14) is called discrete gradient flow iteration (GFI).

Remark 2.1**.**

Provided that the maximal time step $\tau_{\max}$ is sufficiently small, it can be seen that the GFI scheme (14) yields guranteed energy reduction. More precisely, there exists $0<\tau_{\max}=\mathcal{O}(\min\{\beta^{-1},\mathsf{E}(u^{0})^{-\nicefrac{{1}}{{2}}}\})$ such that for all $\tau_{n}\leq\tau_{\max}<2$ it holds that $\mathsf{E}(u^{n+1})\leq\mathsf{E}(\widehat{u}^{n+1})\leq\mathsf{E}(u^{n})$, see [36, Lemma 4.7]. Moreover, if $\tau_{n}\leq 1$ for all $n\geq 0$ then, for any starting value $u^{0}\in S_{\mathbb{H}}$ with $u^{0}\geq 0$, the (full) sequence $\{u^{n}\}_{n\geq 0}$ generated by the GFI (14) satisfies $u^{n}\geq 0$ for all $n\geq 0$ and converges strongly in $\mathbb{H}$ to the unique positive ground state $u_{\mathrm{GS}}$, see [36, Theorem 5.1].

3. Adaptive gradient flow finite element discretization

We now focus on the adaptive spatial discretization of the gradient flow iteration scheme (14).

3.1. Finite element discretization

Consider a sequence of conforming and shape-regular partitions $\{\mathcal{T}_{N}\}_{N\in\mathbb{N}}$ of the domain $\Omega$ into simplicial elements $\mathcal{T}_{N}=\{\kappa\}_{\kappa\in\mathcal{T}_{N}}$ (i.e. triangles for $d=2$ and tetrahedra for $d=3$). Moreover, for a (fixed) polynomial degree $p\in\mathbb{N}$ and any subset $\omega\subset\mathcal{T}_{N}$, we introduce the finite element space

[TABLE]

with $\mathbb{P}_{p}(\kappa)$ signifying the (local) space of all polynomials of maximal total degree $p$ on $\kappa$, $\kappa\in\mathcal{T}_{N}$. Furthermore, similarly as before, we denote by

[TABLE]

the $\mathrm{L}^{2}(\Omega)$-unit sphere in $\mathbb{V}(\omega)$. In the sequel, we apply the notations $\mathbb{X}_{N}:=\mathbb{V}(\mathcal{T}_{N})$ and $S_{N}:=S_{\mathbb{V}(\mathcal{T}_{N})}$. We further denote by $\mathsf{E}_{N}:=\mathsf{E}|_{\mathbb{X}_{N}}$ the restriction of the energy functional $\mathsf{E}$ from (1) to the Galerkin space $\mathbb{X}_{N}$. Then, due to the compactness of $S_{N}$, it exists a minimizer $u_{N}\in\mathbb{X}_{N}$ of $\mathsf{E}_{N}$, i.e. $\mathsf{E}(u_{N})=\min_{v\in S_{N}}\mathsf{E}(v)$, with $(u_{N},1)_{{\rm L}^{2}(\Omega)}\geq 0$. It is not known, however, if $u_{N}$ is unique (up to the sign), see, e.g., [20]. Furthermore, if $\{\mathbb{X}_{N}\}_{N\in\mathbb{N}}$ is a dense family of finite element subspaces of $\mathbb{H}$, then any sequence of minimizers $u_{N}\in S_{N}$ of $\mathsf{E}_{N}$ with $(u_{N},1)_{{\rm L}^{2}(\Omega)}\geq 0$ converges in $\mathrm{H}^{1}(\Omega)$ to the ground state $u_{\mathrm{GS}}$ of $\mathsf{E}$; we refer to [20, Theorem 1] or [50, Theorem 3.1].

3.2. Discrete GFI

Let us define the (space) discrete version of the gradient flow iteration (14) on a finite element subspace $\mathbb{X}_{N}\subset\mathbb{H}$. For $u\in\mathbb{X}_{N}$ we denote by $\mathsf{G}^{N}(u)\in\mathbb{X}_{N}$ the unique solution of

[TABLE]

cf. (7). For given $u\in\mathbb{X}_{N}$, note that the computation of $\mathsf{G}^{N}(u)$ is a standard linear source problem; it can be solved by any linear solver at the disposal of the user. Then, for $n\geq 0$, the space discrete GFI is given by

[TABLE]

with a sequence of discrete time steps $\{\tau_{N}^{n}\}_{n\geq 0}$ as in (14).

Remark 3.1**.**

Consider a fixed mesh $\mathcal{T}_{N}$ and associated approximation space $\mathbb{X}_{N}$. Let $\{u_{N}^{n}\}_{n\geq 0}\subset\mathbb{X}_{N}$ be the sequence generated by the discrete GFI (16) with some initial value $u_{N}^{0}\in S_{N}$. If $\tau_{N}^{n}\leq\tau_{\max}$, with $\tau_{\max}$ as in Remark 2.1, for all $n\geq 0$, then the corresponding energies are strictly monotone decreasing, and there exists a limit energy $\mathsf{E}_{N}^{\star}=\lim_{n\to\infty}\mathsf{E}(u_{N}^{n})$. Furthermore, up to subsequences, we have $u_{N}^{n}\to u_{N}^{\star}$ strongly in $\mathbb{H}$, where $u_{N}^{\star}\in S_{N}$, with $\mathsf{E}(u_{N}^{\star})=\mathsf{E}_{N}^{\star}$, is a discrete eigenfunction of the corresponding GPE, i.e. there is $\lambda_{N}^{\star}$ so that

[TABLE]

We refer to [36, Corollary 4.11] for details.

Remark 3.2**.**

In practical computations, in order to guarantee a positive energy decay in each iteration step, we propose the time step strategy within (16) given by

[TABLE]

where, for $0<s\leq 1$, we write $u_{N}^{n+1}(s)$ to denote the output of the discrete GFI (16) based on the time step $\tau_{N}^{n}=s$ and on the previous approximation $u_{N}^{n}$. We observed in several examples that for the choice $\tau_{N}^{0}=1$, i.e. using $m=0$ above, no time correction was needed; we also refer to [36, Remark 4.8] for a discussion of the fixed time step $\tau=1$. For that reason, and for the sake of keeping the computational cost minimal, we fix the time step $\tau=1$ 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 $\kappa\in\mathcal{T}_{N}$ we consider the open patch $\omega_{\kappa}$ comprising of $\kappa$ and its immediate face-wise neighbours. Moreover, given $\kappa\in\mathcal{T}_{N}$, we define the modified patch $\widetilde{\omega}_{\kappa}$ by uniformly (red) refining the element $\kappa$ into a (fixed) number of subelements; here, we assume that the introduction of any hanging nodes in $\omega_{\kappa}$ is removed by doing (e.g. green) refinements, see Figure 1. We remark that the notions of red and green refinements refer to standard element subdivision techniques in automatic mesh adaptation; further details can be found, e.g., in [41, §4.10.2.2] (or [6], where red refinement is termed bisection-type mesh refinement).

We consider basis functions $\{\xi^{1}_{\kappa},\ldots,\xi^{m_{\kappa}}_{\kappa}\}$ of the locally supported space $\mathbb{V}(\widetilde{\omega}_{\kappa})$. Furthermore, for any given $v\in\mathbb{V}(\mathcal{T}_{N})$, we introduce the extended space

[TABLE]

Suppose we have found an accurate approximation $u_{N}^{n}\in\mathbb{X}_{N}$ of the discrete GPE (17), for some $n\geq 0$. Then, by performing one local discrete GFI-step in $\widehat{\mathbb{V}}(\widetilde{\omega}_{\kappa};u_{N}^{n})\subset\mathbb{H}$ we obtain a new local approximation, denoted by $\widetilde{u}_{N,\kappa}^{n}\in\widehat{\mathbb{V}}(\widetilde{\omega}_{\kappa};u_{N}^{n})$, with $\left\|\widetilde{u}_{N,\kappa}^{n}\right\|_{\mathrm{L}^{2}(\Omega)}=1$. We emphasize that $\widehat{\mathbb{V}}(\widetilde{\omega}_{\kappa};u_{N}^{n})$ has a small dimension, and hence the discrete GFI (16) based on $\widehat{\mathbb{V}}(\widetilde{\omega}_{\kappa};u_{N}^{n})$ entails hardly any computational cost (for instance, for dimension $d=2$ and polynomial degree $p=1$, the dimension of the locally refined space $\widehat{\mathbb{V}}(\widetilde{\omega}_{\kappa};u_{N}^{n})$, cf. Figure 1, is typically 3 or 4).

By modus operandi of the discrete GFI (16), the above construction leads in general, see Remark 3.2, to the (local) energy decay

[TABLE]

for all $\kappa\in\mathcal{T}_{N}$. The value $\Delta\mathsf{E}_{N}^{n}(\kappa)$ indicates the potential energy reduction due to a refinement of the element $\kappa$. This observation motivates the energy-based adaptive mesh refinement procedure outlined in Algorithm 1.

3.4. Adaptive strategy

From a practical viewpoint, once the discrete GFI approximation from (16) is close to a solution of (17), on a given finite element space, we expect that any further GFI steps will no longer contribute an essential decay to the energy in (19). In this case, in order to further reduce the energy, we need to enrich the finite element space appropriately. More specifically, for $N\geq 1$, suppose that we have performed a reasonable number $n\geq 1$ (possibly depending on $N$) of GFI-iterations (16) in $\mathbb{X}_{N-1}$. Consider now a (hierarchically) refined mesh $\mathcal{T}_{N}$ of $\mathcal{T}_{N-1}$, for example, obtained by Algorithm 1. Then we may embed the final guess $u^{n}_{N-1}\in\mathbb{X}_{N-1}$ on the previous space into the enriched finite element space $\mathbb{X}_{N}$ in order to obtain an initial guess on the refined mesh $\mathcal{T}_{N}$:

[TABLE]

For each GFI-iteration $n$ we monitor two quantities. Firstly, we introduce the increment on each iteration given by

[TABLE]

Secondly, we compare $\mathrm{inc}_{N}^{n}$ to the energy loss as compared to the previous mesh refinement, i.e.

[TABLE]

We stop the iteration for $n\geq 1$ as soon as $\mathrm{inc}_{N}^{n}$ becomes small compared to $\Delta\mathsf{E}_{N}^{n}$, i.e. once there is no notable benefit (relatively speaking) in performing any more discrete GFI steps on the current space $\mathbb{X}_{N}$. Specifically, for $n\geq 1$, this is expressed by the bound

[TABLE]

for some parameter $0<\gamma<1$. We implement this procedure in Algorithm 2.

3.5. Computational complexity

We comment on the computational cost of one loop occurring in Algorithm 2. This is essentially comprised of a number of GFI steps on a given Galerkin space $\mathbb{X}_{N}$, and on one adaptive (local) mesh refinement (using Algorithm 1).

• •

The cost of one GF iteration is dominated by the computation of the Riesz-representative $\mathsf{G}^{N}(u_{N}^{n})$ from (15), and, thereby, depends primarily on the dimension $\textrm{dim}(\mathbb{X}_{N})$ of the finite element space $\mathbb{X}_{N}$. Specifically, the solution of (15) (for given $u^{n}_{N}$) amounts to a computational work that scales like $\mathcal{O}(\textrm{dim}(\mathbb{X}_{N})^{\alpha})$, with a parameter $\alpha\geq 1$ depending on the linear solver employed.

• •

The local finite element space $\widehat{\mathbb{V}}(\widetilde{\omega}_{\kappa};u_{N}^{n})$ contains $m_{\kappa}$ local and one global basis function, namely $u_{N}^{n}$, which is the same for all elements $\kappa$. The mesh-refinement procedure contains one GF iteration as described in (16), however, on the local space $\widehat{\mathbb{V}}(\widetilde{\omega}_{\kappa};u_{N}^{n})$. This, in turn, requires the solution of a linear system involving $m_{\kappa}$ local and one global degrees of freedom. In particular, only one entry of the corresponding matrix requires a global integration, which can be computed element-by-element (and, therefore, in parallel); all other matrix entries are represented by local integrals. Finally, note that the entry requiring global integration is the same for all local GF iterations and all elements $\kappa$, and, thus, needs to be computed only once within each step of the loop in Algorithm 2.

A similar observation holds for the computation of the local energy decays. In fact, for $u_{N}^{n}\in\mathbb{X}_{N}$, denoting, as in Algorithm 1, the locally improved approximation by $\widetilde{u}_{N,\kappa}^{n}$, we have the linear combination

[TABLE]

for suitable $\mu,\mu_{i}\in\mathbb{R}$, $i\in\{1,\dotsc,m_{\kappa}\}$. Then, it is elementary to verify that

[TABLE]

where $\mathsf{E}_{\widetilde{\omega}_{\kappa}}(\cdot)$ represents the contribution from the elements comprising the patch $\widetilde{\omega}_{\kappa}$ to the total energy $\mathsf{E}(\cdot)$. We emphasize that the second and third terms on the right-hand side only require the computation of a local integral. Moreover, for the term $\mathsf{E}(\mu u_{N}^{n})$ we observe that $\mathsf{E}(\mu u_{N}^{n})=:\mu^{2}\mathsf{E}^{n}_{N,1}+\mu^{4}\mathsf{E}^{n}_{N,2},$ with

[TABLE]

Thus, for each patch, the energy $\mathsf{E}(\mu u_{N}^{n})$ is a combination of the above two global integrals, which need to be computed once only, and can be split into elementwise contributions.

Altogether, the computational cost for the local GF iterations and corresponding local energy decays therefore scales linearly with the number of elements in the mesh, i.e. with $\textrm{dim}(\mathbb{X}_{N})$, and can be performed fully in parallel. The same observation applies to the element marking and local mesh refinement procedure, which compares to the evaluation of standard residual based a posteriori error estimators.

In summary, we find that the computational work for one loop in Algorithm 2 scales with $\mathcal{O}(\textrm{dim}(\mathbb{X}_{N})^{\alpha})$, where $\alpha\geq 1$ depends on the global finite element solver. In particular, our scheme exhibits a similar complexity as standard adaptive finite element discretization procedures for linear problems (provided that the number of GFI steps remain reasonably modest).

4. Numerical Experiments

We apply Algorithm 2 for some numerical computations in two space dimensions, i.e. $d=2$, with Cartesian coordinates denoted by $\bm{x}=(x,y)\in\mathbb{R}^{2}$. In all examples, we choose the initial guess $u_{0}^{0}\in S_{0}$ such that $u_{0}^{0}(\bm{x})=c$ for any node $\bm{x}$ in the interior of the corresponding (coarse and uniform) initial mesh $\mathcal{T}_{0}$, where $c>0$ is the appropriate constant to fulfil the norm constraint (2); we remark that numerical experiments based on other positive initial guesses (not presented here) have resulted in analogous convergence rates of the minimal energy approximation, thereby indicating a certain robustness with respect to the starting values. Moreover, we set $\theta=0.5$, as well as $\gamma=0.1$ and $\epsilon=10^{-8}$ 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 $\mathbb{X}_{N}$) exceeds $10^{6}$.

4.1. Laplace EVP on $L$-shaped domain

We begin by testing our algorithm for the Laplace eigenvalue problem, which is to find $u\in\mathbb{H}=\mathrm{H}^{1}_{0}(\Omega)$ and $\lambda>0$ such that

[TABLE]

here, $\Omega=(0,2)^{2}\setminus[1,2]\times[0,1]$ is an $L$-shaped domain. This problem is of interest since the eigenfunction to the lowest eigenvalue $\lambda_{1}$, i.e. the ground state, has a singularity at the re-entrant corner point $(1,1)$. From [31] it is known that $9.6397238\leq\lambda_{1}\leq 9.6397239$. Taking the lower bound (divided by 2, cf. (5) for $\beta=0$) as reference value, Figure 2 demonstrates the optimal convergence rate for the minimal energy approximation.

4.2. Linear GPE with smooth potential $V$

We consider the case where $\beta=0$ in (1), and $V(x,y)=\nicefrac{{1}}{{2}}(x^{2}+9y^{2})$ is a smooth function:

[TABLE]

Note that the associated eigenvalue problem (4) is linear for this example. It is known that, for $\Omega=\mathbb{R}^{2}$, the energy of the ground state is given by $\mathsf{E}(u_{\mathrm{GS}})=2$, see, e.g., [7]. We note that the mass of $u_{\mathrm{GS}}$ is essentially concentrated in a vicinity of the origin $\bm{0}$ due to the global minimum of $V$ at $\bm{0}$. Therefore, the restriction to the bounded domain $\Omega:=(-10,10)^{2}\subset\mathbb{R}^{2}$, which we use in our computations, has almost no effect on the minimal value of the ground state energy. Figure 3 (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 $u_{\mathrm{GS}}$ is concentrated, see Figure 3 (right). These results underline that the proposed adaptive gradient flow procedure effectively detects the local behaviour of the model.

4.3. Linear GPE with singular potential $V$.

We perform another experiment with $\beta=0$ in (1), i.e. the associated eigenvalue problem (4) is again linear. In contrast to the previous example, however, we consider a potential $V(\bm{x})=(2|\bm{x}|)^{-1}$ which features a severe point singularity at the origin $\bm{0}=(0,0)$. Specifically, the energy functional is given by

[TABLE]

with $\Omega=\left(\nicefrac{{-1}}{{2}},\nicefrac{{1}}{{2}}\right)^{2}$. This experiment has been conducted already in [45], where the authors obtained an approximated minimal eigenvalue $\lambda_{\mathrm{GS}}=25.934923921168$ for the corresponding GPE (4); we compare our results to this reference value (divided by 2 on account of (5)). We can see from Figure 4 that the proposed Algorithm 2 achieves a sequence of energy approximations for the ground state which decays at an (almost) optimal rate.

4.4. Linear GPE with potential wells

We test a final example with $\beta=0$. The potential $V$ is given by the sum of four Gaussian bells, see Figure 5 (left). This experiment is borrowed from [44, Experiment 4.2], however, with a (constant) shift such that $V\geq 0$ in the underlying domain $\Omega=(0,2\pi)^{2}$. As we can see from Figure 5 (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 [44]. Moreover, the energy-based adaptive mesh refinement has properly resolved the two local hills featured in the ground state $u_{\mathrm{GS}}$, see Figure 6.

The authors from [44] have computed an approximated minimal eigenvalue $\lambda_{\mathrm{GS}}=16.6879$, 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 7 (left). We remark, for this example, that the performance of Algorithm 2 crucially depends on the choice of the parameter $\gamma$ in Line 6. Indeed, if $\gamma=0.5$ (instead of $\gamma=0.1$) is selected, then the numerical results exhibit a considerably less favourable asymptotic convergence regime, see Figure 7 (right). An explanation can be inferred from Figure 8: We see that the choice $\gamma=0.1$ leads to a significantly higher number of gradient flow steps on each Galerkin space $\mathbb{X}_{N}$, which seems to be essential for the effective numerical solution of this problem.

4.5. Nonlinear GPE with harmonic confinement potential

We now consider a nonlinear Bose-Einstein condensate, i.e. $\beta=1000\gg 0$ in (1), and use the smooth potential $V(\bm{x})=\nicefrac{{1}}{{2}}|\bm{x}|^{2}$:

[TABLE]

the domain is given by $\Omega=(-6,6)^{2}$. This experiment was conducted previously in [36], where an approximation $\mathsf{E}(u_{\mathrm{GS}})\approx 11.9860647$ for the energy of the ground state has been documented. 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 $\mathsf{E}(u_{\mathrm{GS}})\approx\underline{11.98605}121\dotsc$ based on an adaptively refined mesh with $\mathcal{O}(10^{7})$ degrees of freedom; here, the underlined digits are stable (i.e. the computations indicate that they do not change any more as the iterations continue). Using this as reference value, we obtain optimal convergence for the approximation of the energy of the ground state, see Figure 9.

In order to study the dependence of the computational work on the strength of the nonlinearity, we have performed the same experiment for different values of $\beta$. As can be seen from Table 4.5, our algorithm does not seem to be sensitive with respect this parameter. Indeed, after an initial phase, the number of GF iterations on a given mesh is independent of $\beta$. Moreover, for $\gamma=0.5$, this holds true already from the start, i.e. even for the initial mesh (not displayed here). We remark that analogous results have been observed recently in [49] for a multigrid method.

[TABLE]

\captionof

tableExperiment 4.5. Number of GF iterations on the first 15 finite element spaces for different values of $\beta$.

4.6. Nonlinear GPE with optical lattice potential

As before, we choose $\Omega=(-6,6)^{2}$ and $\beta=1000$, with an oscillating potential function $V$, see Figure 11 (left) for its contour plot. More precisely, the energy functional is given by

[TABLE]

This experiment was also considered in [36] with an asserted approximation $\mathsf{E}(u_{\mathrm{GS}})\approx 30.40965$ of the ground state energy. Based on a sufficiently fine uniform initial mesh (with $128$ elements), our algorithm yields the approximation $\mathsf{E}(u_{\mathrm{GS}})\approx\underline{30.387}533\dotsc$ for $\mathcal{O}(10^{7})$ degrees of freedom. In Figure 10 (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.

As in Example 4.4, we have also run this experiment for both values $\gamma=0.1$ and $\gamma=0.5$. In contrast to the previous test, we observe no considerable difference in the performance of the corresponding computations. Indeed, except for the first two meshes, the number of gradient flow steps on each Galerkin space is the same for both cases, see Figure 12.

4.7. Nonlinear energy functional with a nonsymmetric potential $V$

Finally, we revisit the test problem [11, Example 4.3.II]:

[TABLE]

with the symmetric domain $\Omega=(-8,8)^{2}$. In [11] an approximated value of $\mathsf{E}(u_{\mathrm{GS}})\approx 5.8507$ for the energy of the ground state $u_{\mathrm{GS}}$ has been obtained. Like in the previous experiments, for an adaptively refined mesh with $\mathcal{O}(10^{7})$ degrees of freedom, a smaller value of $\mathsf{E}(u_{\mathrm{GS}})\approx\underline{5.85058}738\dotsc$ is obtained. Applying this value as reference ground state energy, we observe optimal convergence in Figure 13.

Furthermore, following the previous Experiment 4.5, we have run this test problem for different values of $\beta$, with similar results, see Table 4.7.

[TABLE]

\captionof

tableExperiment 4.7. Number of GF iterations on the first 15 finite element spaces for different values of $\beta$.

5. Conclusions

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. Moreover, in our experiments, the effect of the parameter $\beta$ (steering the strength of the nonlinearity) on the computational work seems to be negligible.