Sobolev gradient flow for the Gross-Pitaevskii eigenvalue problem: global convergence and computational efficiency

TL;DR

This paper introduces a new Sobolev gradient flow for the Gross-Pitaevskii eigenvalue problem, demonstrating exponential convergence to ground states and an efficient numerical method that outperforms existing approaches.

Contribution

It proposes a novel normalized Sobolev gradient flow with a time-dependent inner product, providing global convergence and an efficient discretization for the nonlinear eigenvalue problem.

Findings

Flow converges exponentially fast to eigenfunctions.

Numerical method reduces energy at each step and converges globally.

Method is competitive with established discretizations.

Abstract

We propose a new normalized Sobolev gradient flow for the Gross-Pitaevskii eigenvalue problem based on an energy inner product that depends on time through the density of the flow itself. The gradient flow is well-defined and converges to an eigenfunction. For ground states we can quantify the convergence speed as exponentially fast where the rate depends on spectral gaps of a linearized operator. The forward Euler time discretization of the flow yields a numerical method which generalizes the inverse iteration for the nonlinear eigenvalue problem. For sufficiently small time steps, the method reduces the energy in every step and converges globally in $H^1$ to an eigenfunction. In particular, for any nonnegative starting value, the ground state is obtained. A series of numerical experiments demonstrates the computational efficiency of the method and its competitiveness with established discretizations arising from other gradient flows for this problem.

Figures

Tables

Table 1. Table 1 : Approximate ground state energies E GS subscript 𝐸 GS E_{\mbox{\rm\tiny GS}} and corresponding ground state eigenvalues λ GS subscript 𝜆 GS \lambda_{\mbox{\rm\tiny GS}} for Model Problem 1 with h = 12 ⋅ 2 − 8 ℎ ⋅ 12 superscript 2 8 h=12\cdot 2^{-8} and different values for β 𝛽 \beta .

$\beta$	$E_{\text{GS}}$	${\lambda }_{\text{GS}}$
10	0.79620688	2.06380
100	1.97298868	5.75977
1000	5.99303235	17.9771

Table 2. Table 2 : Model Problem 1: computation of ground states. The table shows the number of iterations obtained for the various methods for τ = 0.5 𝜏 0.5 \tau=0.5 and τ = 1 𝜏 1 \tau=1 . The entry “ ∞ \infty ” means that the iteration did not converge. The spatial mesh size was selected as h = 12 ⋅ 2 − 6 ℎ ⋅ 12 superscript 2 6 h=12\cdot 2^{-6} . The entries in brackets show the iteration count for higher spatial resolution h = 12 ⋅ 2 − 8 ℎ ⋅ 12 superscript 2 8 h=12\cdot 2^{-8} .

$\tau$	$\beta$	GF$L^2$	GF$H^1$	GF$a_0$	GF$a_z$
1.0	10	9 (9)	$\infty$ ($\infty$)	7 (7)	6 (7)
0.5	10	11 (11)	$\infty$ ($\infty$)	14 (14)	14 (14)
1.0	100	11 (12)	$\infty$ ($\infty$)	$\infty$ ($\infty$)	9 (9)
0.5	100	13 (13)	$\infty$ ($\infty$)	$\infty$ ($\infty$)	18 (18)
1.0	1000	15 (15)	$\infty$ ($\infty$)	$\infty$ ($\infty$)	11 (11)
0.5	1000	16 (16)	$\infty$ ($\infty$)	$\infty$ ($\infty$)	22 (22)

Table 3. Table 3 : Model Problem 1: computation of ground states. The table shows the minimum number of iterations N 𝑁 N that the methods required to reach the error tolerance. Alongside N 𝑁 N we list one possible step size τ 𝜏 \tau for which this number is reached. The spatial mesh size is fixed with h = 12 ⋅ 2 − 6 ℎ ⋅ 12 superscript 2 6 h=12\cdot 2^{-6} .

	$\tau$	$N$	$\tau$	$N$	$\tau$	$N$
	$\beta =10$		$\beta =100$		$\beta =1000$
GF$L^2$	$5.0$	$7$	$5.0$	$10$	$5.0$	$15$
GF$H^1$	$0.25$	$41$	$0.31$	$25$	$0.12$	$48$
GF$a_0$	$1.0$	$7$	$0.34$	$24$	$0.1$	$109$
GF$a_z$	$1.2$	$5$	$1.1$	$8$	$1.1$	$10$

Equations

$$-△z^{\ast }+V{z}^{\ast }+\beta |z^{\ast }{|}^2{z}^{\ast }=\lambda z^{\ast }.$$

$$E(v):=\frac{1}{2}{\int }_{{\mathbb{R}}^d}|\nabla v{|}^2+V|v{|}^2+\frac{\beta }{2}|v{|}^{4}dx,$$

$$z^{\prime }(t)=-P_{z,X}({\nabla }_{X}E(z(t))),$$

$$E(v):=\frac{1}{2}{\int }_{\Omega }|\nabla v{|}^2+V|v{|}^2+\frac{\beta }{2}|v{|}^{4}dx.$$

$$⟨E^{\prime }(v),w⟩={\int }_{\Omega }\nabla v\cdot \nabla w+Vvw+\beta |v{|}^{2}vwdx.$$

$$⟨E^{\prime }(z^{\ast }),v⟩={\lambda }^{\ast }(z^{\ast },v{)}_{L^2(\Omega )}\text{mbox}forallv\in H_0^1(\Omega ).$$

$${\lambda }^{\ast }=2E(z^{\ast })+\frac{\beta }{2}\Vert z^{\ast }{\Vert }_{L^4(\Omega )}^4.$$

$$E(z_{\text{mbox}\mathrm{GS}})=\inf \{E(z)|z\in H_0^1(\Omega )\text{mbox}\mathrm{and}\Vert z{\Vert }_{L^2(\Omega )}=1\}.$$

$$({\nabla }_{X}E(z),v{)}_X=⟨E^{\prime }(z),v⟩\text{mbox}forallv\in X.$$

$$T_{z,X}:=\{v\in X|(v,z{)}_{L^2(\Omega )}=0\}.$$

$$\frac{1}{2}\frac{\text{mbox}d}{\text{mbox}dt}\Vert z(t){\Vert }_{L^2(\Omega )}^2=(z^{\prime }(t),z(t){)}_{L^2(\Omega )}=(P(z(t)),z(t){)}_{L^2(\Omega )}=0,$$

$$(P_{z,X}(v),\psi {)}_X=(v,\psi {)}_X\text{mbox}forall\psi \in T_{z,X},$$

$$P_{z,X}(v)=v-\frac{(z,v{)}_{L^2(\Omega )}}{(z,R_X(z){)}_{L^2(\Omega )}}{R}_X(z).$$

$$z^{\prime }(t)=-(P_{z(t),X}\circ {\nabla }_{X}E)(z(t))\text{mbox}fort\ge 0\text{mbox}andz(0)=z_0.$$

$$P_{z,L^2}(v)=v-\frac{(z,v{)}_{L^2(\Omega )}}{(z,z{)}_{L^2(\Omega )}}z.$$

$${\nabla }_{L^2}E(z)=-△z+Vz+\beta |z{|}^{2}z.$$

$$z^{\prime }(t)=-A_{z}z+⟨A_{z}z,z⟩\frac{z}{\Vert z{\Vert }_{L^2(\Omega )}^2}\text{mbox}fort\ge 0$$

$$z^{n+1}={\alpha }_n(1+\tau A_{z^n}{)}^{-1}{z}^n\text{mbox}where{\alpha }_n:=\Vert (1+\tau A_{z^n}{)}^{-1}{z}^n{\Vert }_{L^2(\Omega )}^{-1}.$$

$$(\nabla R_{H^1}(z),\nabla v{)}_{L^2(\Omega )}=(z,v{)}_{L^2(\Omega )}\text{mbox}forallv\in H_0^1(\Omega ).$$

$$z^{\prime }=-{\nabla }_{H^1}E(z)+\frac{({\nabla }_{H^1}E(z),z{)}_{L^2(\Omega )}}{(R_{H^1}(z),z{)}_{L^2(\Omega )}}{R}_{H^1}(z)\text{mbox}fort\ge 0\text{mbox}andz(0)=z_0.$$

$${\hat{z}}^{n+1}=z^n-\tau {\nabla }_{H^1}E(z^n)+\tau \frac{({\nabla }_{H^1}E(z^n),z^n{)}_{L^2(\Omega )}}{(R_{H^1}(z^n),z^n{)}_{L^2(\Omega )}}{R}_{H^1}(z^n)$$

$$a_0(v,w):=(\nabla v,\nabla w{)}_{L^2(\Omega )}+(Vv,w{)}_{L^2(\Omega )}.$$

$$a_0(R_{a_0}(z),v)=(z,v{)}_{L^2(\Omega )}\text{mbox}forallv\in H_0^1(\Omega ).$$

$$z^{\prime }$$

$${\hat{z}}^{n+1}=(1-\tau )z^n-\tau R_{a_0}(\beta |z^n{|}^2{z}^n)+\tau \frac{(z^n,z^n+R_{a_0}(\beta |z^n{|}^2{z}^n){)}_{L^2(\Omega )}}{(z^n,R_{a_0}{z}^n{)}_{L^2(\Omega )}}{R}_{a_0}{z}^n,$$

$$a_z(v,w):={\int }_{\Omega }\nabla v\cdot \nabla w+Vvw+\beta |z{|}^{2}vwdx$$

$$z^{\prime }(t)=-P_{z(t)}{\nabla }_{a_{z(t)}}E(z(t))=-P_{z(t)}(z(t)).$$

$$P_z(h)=h-\frac{(z,h{)}_{L^2(\Omega )}}{(z,{\mathcal{G}}_{z}z{)}_{L^2(\Omega )}}{\mathcal{G}}_{z}z,$$

$$a_z({\mathcal{G}}_z(f),v)=(f,v{)}_{L^2(\Omega )}$$

$$z^{\prime }(t)=-z(t)+{\gamma }_{z(t)}{\mathcal{G}}_{z(t)}z(t),\text{mbox}where{\gamma }_z:=\frac{(z,z{)}_{L^2(\Omega )}}{a_z({\mathcal{G}}_{z}z,{\mathcal{G}}_{z}z)}>0.$$

Full text

1 Introduction

A Bose-Einstein condensate (BECs) is an extreme state of matter formed by a dilute gas of bosons at ultra-cold temperatures, very close to absolute zero [16, 25, 29, 42]. In a BEC, individual particles (i.e. their wave packages) overlap, lose their identity, and form one single “super atom”. BECs allow to study macroscopic quantum phenomena such as superfluity (i.e. the frictionless flow of a fluid) on an observable scale. This is why BECs are a very relevant research area of modern quantum physics [4, 1, 32, 38, 40]. For a general mathematical description of BECs and corresponding analytical results we refer exemplarily to [2, 10, 18, 39, 42].

In this paper we consider stationary states of a BEC modeled by the Gross-Pitaevskii eigenvalue problem (GPE) in real-valued variables. In non-dimensional form, the GPE seeks $L^{2}$-normalized eigenfunctions $z^{\ast}\in H^{1}(\mathbb{R}^{d})$ and corresponding eigenvalues $\lambda^{\ast}\in\mathbb{R}$ such that

[TABLE]

In the context of Bose-Einstein condensates, a solution $z^{\ast}$ represents a stationary quantum state of the condensate, $|z^{\ast}|^{2}$ is the corresponding density and $\lambda^{\ast}$ the so-called chemical potential. The function $V$ represents an external confining potential and the parameter $\beta$ depends on physical properties of the particles that form the BEC. Its sign determines the type of particle interactions. In this paper, we shall only consider the defocusing GPE, which covers the regime $\beta\geq 0$, resembling repulsive particle interactions. The normalization constraint $\int_{\mathbb{R}^{d}}|z^{\ast}|^{2}\operatorname{\,d}x=1$ is such that the total mass of the condensate equals the number of constituting particles (with probability $1$).

The numerical solution of the stationary GPE has been studied extensively in recent years; see e.g. [3, 8, 11, 12, 13, 14, 15, 19, 20, 21, 22, 23, 24, 26, 27, 28, 31, 33, 35, 36, 37, 43, 44, 45] and the references therein. Typically, the problem is rephrased in terms of the energy functional

[TABLE]

where one is interested in finding the critical points of $E$ under the normalization constraint $\|v\|_{L^{2}(\mathbb{R}^{d})}=1$. The unique global minimizer (the state of minimal energy) is called the ground state, whereas all other critical points are called excited states. The identification of critical points of $E$ can be accomplished by the construction of appropriate gradient flows of the form

[TABLE]

where $\nabla_{X}E$ is the Sobolev gradient of the energy functional $E$ with respect to some inner product $(\cdot,\cdot)_{X}$ and where $P_{z,X}$ is the projection onto the tangent space associated with the normalization constraint. Depending on the choice of $(\cdot,\cdot)_{X}$ and the numerical time integration of the arising gradient flow, several numerical methods arise (cf. [26]). Presumably, the most popular method in the context of the GPE is the Discrete Normalized Gradient Flow (DNGF) [12] which is based on the choice $(\cdot,\cdot)_{X}=(\cdot,\cdot)_{L^{2}}$ of the $L^{2}$-inner product and a backward Euler-type time discretization with explicit treatment of the nonlinear term. Other approaches combine a forward Euler discretization with the choice $(\cdot,\cdot)_{X}=(\cdot,\cdot)_{H^{1}}$ [37, 44] or the choice $(\cdot,\cdot)_{X}=(\nabla\cdot,\nabla\cdot)_{L^{2}}+(V\cdot,\cdot)_{L^{2}}$ [26]. These examples and their discrete version are briefly discussed in Section 3. For further variants, we refer to [3, 8, 14, 27, 33, 43].

Although the aforementioned schemes for the GPE are empirically successful, their numerical analysis lacks a proof of global convergence to a critical point of $E$ and any quantification of convergence rates. There is not even a proof of monotonic energy dissipation of the iteration in analogy to the continuous gradient flow (1). The only result that comes close is for DNGF (based on the $L^{2}$-gradient) [12]. In the absence of any spatial discretization, the reduction of a modified energy is shown which deviates from the exact energy by a term of the form $\tfrac{\beta}{4}\|v\|_{L^{4}}^{4}$. Since this result exploits elliptic regularity theory, its generalization to a fully discrete setting involving e.g. a finite element discretization is not straight forward.

In Section 4 of this paper, we present a new choice for the Sobolev gradient, where the inner product $(\cdot,\cdot)_{X}$ is not fixed, but evolves with time. It is selected in such a way that the Sobolev gradient equals the identity, thus leading to an optimal preconditioning of the flow. We show that the arising continuous gradient flow of the form (1) is well-posed. Thanks to the optimal preconditioning, the problem can be discretized by the forward Euler scheme (cf. Section 4). The time-discrete method reduces the (correct) energy monotonically and converges globally in $H^{1}$ to a critical point of $E$ for sufficiently small time steps. These unique results remain valid even after Galerkin discretization in space. Furthermore, in Section 5 we prove that, for any non-negative initial value $z_{0}$, the method must necessarily converge to a strictly positive eigenfunction of the GPE. Since there exist no positive excited state, the method is guaranteed to converge to the ground state whenever $z_{0}\geq 0$.

Exponential convergence of the new discrete Gradient flow with respect to the number of iterations (i.e. reduction of the error by a fixed factor in each step) remains open but is observed in numerical experiments. It is worth mentioning that, for a particular choice of the time step, the method recovers the inverse iteration for the nonlinear eigenvalue problem. Moreover, for this very time step, the method is equivalent to DNGF which indicates its competitiveness with the established approaches for the GPE. In some scenarios we even observe superior performance (see Section 6). This is particularly true when the time step is chosen adaptively by some standard line search strategy which appears to be cost neutral.

2 Model problem and established gradient flows

We shall introduce the precise setup of the model problem of this paper and briefly recall the projected $L^{2}$- and $H^{1}$-Sobolev gradient flows at hand. Note that all functions and functionals considered in this paper are real-valued.

2.1 Gross-Pitaevskii eigenvalue problem

Since confinement potentials $V$ cause a localization of stationary states, it is common to consider the GPE on a bounded domain $\Omega\subset\mathbb{R}^{d}$, for $d=1,2,3$, together with a homogenous Dirichlet boundary condition. In addition to the boundedness, we shall also assume that $\Omega$ is either a convex Lipschitz domain or a domain with a smooth boundary. The latter assumption is natural in this context and prevents singular behavior of stationary states at the artificial boundary. We also assume that the nonlinearity is defocusing, i.e., $\beta\geq 0$, and that the potential is bounded almost everywhere, i.e. $V\in L^{\infty}(\Omega)$. Without loss of generality, we assume that $V\geq 1$ a.e., as a constant shift of $V$ would not affect the eigenfunctions but only shift the spectrum accordingly. Note that this assumption implies that all eigenvalues satisfy $\lambda^{\ast}>1$. We note that we only use $V\geq 1$ instead of $V\geq 0$ to avoid a repeated usage of the Poincaré inequality in our estimates.

We define the non-negative energy for a function $v\in H^{1}_{0}(\Omega)$ by

[TABLE]

The energy functional is strictly convex and Fréchet differentiable, where the first Fréchet derivative is given by

[TABLE]

Here, $\langle\cdot,\cdot\rangle$ denotes the dual pairing between $H^{-1}(\Omega)$ and $H^{1}_{0}(\Omega)$. The Gross-Pitaevskii eigenvalue problem (GPE) seeks the critical points $z^{\ast}\in H^{1}_{0}(\Omega)$ of the energy functional $E$ subject to the constraint $\|z^{\ast}\|_{L^{2}(\Omega)}=1$. A function $z^{\ast}$ is a critical point if there is a $\lambda^{\ast}\in\mathbb{R}$ such that

[TABLE]

For an $L^{2}$-normalized eigenfunction $z^{\ast}$, the energy is related to the corresponding eigenvalue $\lambda^{\ast}$ through the equality

[TABLE]

Classical Ljusternik-Schnirelman theory (cf. [46]) for even, positive, convex functionals guarantees that problem (2) has infinitely many eigenvalues $0<\lambda_{1}^{\ast}<\lambda_{2}^{\ast}\leq\lambda_{3}^{\ast}\leq\cdots<\infty$. Of particular interest is the ground state of $E$ (the global minimizer) with ground state eigenvalue $\lambda_{\mbox{\rm\tiny GS}}=\lambda_{1}^{\ast}$. The following result can be e.g. found in [21].

Proposition 2.1**.**

Under the general assumptions of this paper there exists a ground state $z_{\mbox{\rm\tiny GS}}\in H^{1}_{0}(\Omega)$ with $\|z_{\mbox{\rm\tiny GS}}\|_{L^{2}(\Omega)}=1$ such that

[TABLE]

The (normalized) ground state is unique up to its sign, it is Hölder-continuous on $\overline{\Omega}$ and it satisfies $|z_{\mbox{\rm\tiny GS}}|>0$ in $\Omega$. Furthermore, the Lagrange multiplier $\lambda_{\mbox{\rm\tiny GS}}$ given by (3), is the smallest eigenvalue of the GPE (2) with corresponding eigenfunction $z_{\mbox{\rm\tiny GS}}$. This ground state eigenvalue $\lambda_{\mbox{\rm\tiny GS}}$ is simple.

We stress the nontrivial observation that a normalized eigenfunction to the smallest eigenvalue of the GPE is always a global minimizer of the energy functional $E$. We are not aware of any result that ensures that the ordering of the eigenvalues $\lambda_{n}^{\ast}$ by size still corresponds with the ordering of the energies $E(z_{n}^{\ast})$ by size for $n>1$.

Other than the ground state, excited states are not unique (up to sign) in general. E.g., on a circular domain with an isotropic quadratic potential, the eigenvalues that correspond to excited states can even have an infinite multiplicity due to rotational invariance of $E$.

2.2 Projected Sobolev gradient flows

We shall briefly recall the basic concept of projected gradient flows. For a detailed introduction to the topic in the context of the Gross-Pitaevskii equation, we refer to [37].

We consider the energy functional $E$ along with a Hilbert space $H^{1}_{0}(\Omega)\subset X\subset L^{2}(\Omega)$ with inner product $(\cdot,\cdot)_{X}$ as the energy dissipation mechanism. Various choices for $X$ are possible and lead to different gradient flows. With this, let $\nabla_{X}E(z)$ denote the Riesz-representative of $E^{\prime}(z)$ in the space $X$, i.e., $\nabla_{X}E(z)\in X$ satisfies

[TABLE]

The operator $\nabla_{X}E:H^{1}_{0}(\Omega)\rightarrow X$ is called the Sobolev gradient of $E$ with respect to $X$. For the sake of mass conservation along the flow, we define the tangent space of the constraint $\|z\|_{L^{2}(\Omega)}^{2}=1$ in $X$ by

[TABLE]

Note that $T_{z,X}$ is the null space of the Fréchet derivative of the functional $v\mapsto\|v\|_{L^{2}(\Omega)}^{2}$ on $X$ evaluated at $z$. If $z^{\prime}=P(z)$ for some operator $P:H^{1}_{0}(\Omega)\rightarrow T_{z,X}$, then we have $z^{\prime}\in X\subset L^{2}(\Omega)$ and, hence,

[TABLE]

i.e., we have mass conservation with $\|z(t)\|_{L^{2}(\Omega)}^{2}=\|z(0)\|_{L^{2}(\Omega)}^{2}$ for all $t\geq 0$. This motivates to seek the best approximation of the Sobolev gradient $\nabla_{X}E(z)$ in the tangent space $T_{z,X}$. The $X$-orthogonal projection $P_{z,X}(v)\in T_{z,X}$ of $v$ onto $T_{z,X}$ is given by

[TABLE]

and be expressed in terms of the Riesz-representative $R_{X}(z)\in X$ of $z$ in $X$ by

[TABLE]

Given some sufficiently smooth initial value $z_{0}$, the projected Sobolev gradient flow is then characterized by

[TABLE]

We shall discuss three choices of spaces $(X,(\cdot,\cdot)_{X})$ along with suitable time discretizations in Sections 2.2.1–2.2.3 below.

2.2.1 Projected $L^{2}$-gradient flow

The most popular choice $X=L^{2}(\Omega)$ leads to the ordinary $L^{2}$-gradient flow. In this case, $R_{L^{2}}(z)=z$ and the projection reads

[TABLE]

The $L^{2}$-gradient is given by the Gross-Pitaevskii differential operator. In particular, for any $z\in H^{2}(\Omega)$, we have

[TABLE]

With $A_{z}v:=-\triangle v+V\hskip 2.0ptv+\beta\hskip 2.0pt|z|^{2}v$, the projected $L^{2}$-gradient flow is given by

[TABLE]

and some initial value $z(0)=z_{0}\in H^{1}_{0}(\Omega)\cap H^{2}(\Omega)$. This is the normalized gradient flow of [12, Section 2.3]. For focusing nonlinearities, Faou and Jézéquel [31] proved exponential $L^{2}$-convergence of the flow to an eigenfunction if the starting value was selected sufficiently close. This is to our best knowledge the only convergence result for the projected $L^{2}$-gradient flow in the context of nonlinear Schrödinger equations. Applying a certain first order splitting method together with a semi-implicit backward Euler discretization with step size $\tau>0$, the DNGF approach is obtained [12]. For the sake of consistent notation we will refer to it as GF$L^{2}$.

Definition 2.2** (Method: GF$L^{2}$ (known as DNGF)).**

Let $z^{0}\in H^{1}_{0}(\Omega)$ be given. Then the GF$L^{2}$ iteration for $n\geq 0$ reads

[TABLE]

By construction, the continuous flow (5) is mass-conservative and energy-dissipative. However, on the time-discrete level, energy dissipation is only established for a modified energy (cf. [12, Lemma 2.10]) that can be seen as an eigenvalue functional. There is no proof of global convergence in $H^{1}(\Omega)$ to a critical point of $E$, neither for the GF$L^{2}$ iteration nor the continuous flow. However, local convergence in $H^{1}(\Omega)$ was established in [31] for focusing nonlinearities, i.e. for $\beta=-1$.

Promising computational improvements of GF$L^{2}$ by using sophisticated preconditioners were obtained in [8].

2.2.2 Projected $H^{1}$-Sobolev gradient flow

In the second example, we consider the choice $X=H^{1}_{0}(\Omega)$ equipped with the standard inner product $(\nabla\cdot,\nabla\cdot)_{L^{2}(\Omega)}$ (cf. [37]). Then the Ritz-projection $R_{H^{1}}:L^{2}(\Omega)\rightarrow H^{1}_{0}(\Omega)$ is characterized by

[TABLE]

The Sobolev gradient $\nabla_{H^{1}}E(z)\in H^{1}_{0}(\Omega)$ is defined according to (4) and the continuous flow reads

[TABLE]

[37, Thms. 5 and 6] reports well-posedness and exponential convergence of the flow to a critical point of $E$ in $H^{1}$. The discretization of the continuous flow (7) using the forward Euler method leads to the GF$H^{1}$ approach.

Definition 2.3** (Method: GF$H^{1}$).**

Let $z^{0}\in H^{1}_{0}(\Omega)$ be given. Then the GF$H^{1}$ iteration for $n\geq 0$ reads

[TABLE]

with the normalization $z^{n+1}=\hat{z}^{n+1}/\|\hat{z}^{n+1}\|_{L^{2}(\Omega)}$ after each iteration.

To the best of our knowledge, there is neither a proof of energy dissipation of the GF$H^{1}$ iteration (8) nor convergence in $H^{1}(\Omega)$ to a critical point of $E$. In numerical experiments, the desired properties are observed for sufficiently small $\tau$ (cf. [37, 26] or the numerical experiments in Section 6.1).

2.2.3 Projected $a_{0}$-Sobolev gradient flow

In the final example we choose $X=H^{1}_{0}(\Omega)$ again, but equip it with an inner product that incorporates the potential $V$. We set $(\cdot,\cdot)_{X}=a_{0}(\cdot,\cdot)$, where

[TABLE]

This choice was proposed in [26] in a more general setup that involves angular momentum rotation. Define the Ritz-projection $R_{a_{0}}(z)\in H^{1}_{0}(\Omega)$ by

[TABLE]

Then, the continuous projected gradient flow reads

[TABLE]

completed by the initial condition $z(0)=z_{0}\in H^{1}_{0}(\Omega)$. Well-posedness of this gradient flow follows from [26, Theorem 3.2].

A forward Euler discretization leads to the following method

Definition 2.4** (Method: GF$a_{0}$).**

Let $z^{0}\in H^{1}_{0}(\Omega)$ be given. Then the GF$a_{0}$ iteration for $n\geq 0$ reads

[TABLE]

together with the normalization $z^{n+1}=\hat{z}^{n+1}/\|\hat{z}^{n+1}\|_{L^{2}(\Omega)}$.

Proofs of energy reduction or the convergence of $z^{n}$ to a stationary point of $E$ are not available in the literature.

3 Continuous Projected $a_{z}$-Soblev Gradient Flow

In this section we propose and analyze a new Sobolev gradient flow in $X=H^{1}_{0}(\Omega)$ based on an inner product that changes with the flow itself. For any $z\in H^{1}_{0}(\Omega)$, we define the weighted energy inner product $a_{z}(\cdot,\cdot)$ by

[TABLE]

for $v,w\in H^{1}_{0}(\Omega)$. Since $\langle E^{\prime}(z),v\rangle=a_{z}(z,v)$ for any $z\in H^{1}_{0}(\Omega)$, the Sobolev gradient of $E(z)$ with respect to $a_{z}(\cdot,\cdot)$ is the identity, i.e., $\nabla_{a_{z}}E(z)=z$. The gradient flow of $E$ with respect to $\nabla_{a_{z}}$ projected into the tangent space associated with mass constraint $\int_{\Omega}|z|^{2}\operatorname{\,d}x=1$ is thus characterized by

[TABLE]

The projection $P_{z}$ can be written as

[TABLE]

where $\mathcal{G}_{z}$ is just the Green’s operator (or Ritz-projection) associated with the time-dependent inner product $a_{z}(\cdot,\cdot)$, i.e., for any $f\in L^{2}(\Omega)$, $\mathcal{G}_{z}(f)\in H^{1}_{0}(\Omega)$ satisfies

[TABLE]

for all $v\in H^{1}_{0}(\Omega)$. Altogether, this yields the following projected gradient flow problem.

Definition 3.1** (Projected $a_{z}$-Sobolev Gradient Flow).**

Given $z_{0}\in H^{1}_{0}(\Omega)$ with $\|z_{0}\|_{L^{2}(\Omega)}=1$, find a differentiable function $z\in C^{1}([0,\infty);H^{1}_{0}(\Omega))$ with $z(0)=z_{0}$ such that, for all $t>0$,

[TABLE]

The subsequent theorem states the well-posedness of this flow and all its important properties.

Theorem 3.2**.**

For any initial value $z_{0}\in H^{1}_{0}(\Omega)$ with $\|z_{0}\|_{L^{2}(\Omega)}=1$, there exists a unique global solution $z$ to the Sobolev gradient flow problem stated in Definition 3.1. The flow is mass-conservative, i.e., $\|z(t)\|_{L^{2}(\Omega)}=1$ for all $t\in\mathbb{R}_{\geq 0}$, and energy-dissipative, i.e. $E(\hskip 1.0ptz(t)\hskip 1.0pt)\leq E(\hskip 1.0ptz(s)\hskip 1.0pt)$ for all $0\leq s\leq t<\infty$. Moreover, $z$ converges globally in $H^{1}(\Omega)$ to an eigenfunction $z^{\ast}$ with eigenvalue $\lambda^{\ast}=\gamma_{z^{\ast}}$ of the Gross-Pitaevskii equation (2).

If $z^{\ast}=z_{\mbox{\rm\tiny GS}}$ is the unique positive ground state eigenfunction from Proposition 2.1 (which can be obtained by selecting $z_{0}$ as a nonnegative functions) then the convergence rate is asymptotically exponential in the following sense. For all $0<\delta\leq 1$, there exists some $c_{\delta}>0$ and a finite time $0<t_{\delta}<\infty$, such that for all $t\geq t_{\delta}$

[TABLE]

Here, $\lambda_{\mbox{\rm\tiny GS}}>0$ is the ground state eigenvalue and $\mu_{2}>\lambda_{\mbox{\rm\tiny GS}}$ is the second eigenvalue of the linearized eigenvalue problem seeking $w_{i}\in H^{1}_{0}(\Omega)$ with $\|w_{i}\|_{L^{2}(\Omega)}=1$ and $\mu_{i}\in\mathbb{R}$ such that

[TABLE]

The remainder of this section is devoted to the proof of the theorem.

3.1 Energy decay, mass conservation and local well-posedness

This subsection shows that a well-defined flow $z$ is energy-diminishing and mass-conserving, as expected.

Lemma 3.3** (Mass conservation and energy reduction).**

We consider the weighted Sobolev gradient flow of Definition 3.1. If $z(t)$ is well-defined on an interval $[0,T)$ for some $T>0$ then, for all $0\leq t\leq t^{\prime}<T$,

[TABLE]

Proof.

Let $t\in[0,T)$ be arbitrary but fixed. Noting $a_{z}(\mathcal{G}_{z}z,\mathcal{G}_{z}z)=(z,\mathcal{G}_{z}z)_{L^{2}(\Omega)}$ and testing in the $L^{2}$-variational formulation of (12) with $z$ yields

[TABLE]

This implies conservation of mass. By definition, $\mathcal{G}_{z}z\in H^{1}_{0}(\Omega)$ and, hence,

[TABLE]

We can therefore use $z^{\prime}$ as a test function in the energy-inner product to obtain

[TABLE]

This implies

[TABLE]

The combination of the previous equalities readily yields

[TABLE]

which shows that energy is reduced along the flow. ∎

To prove local existence of $z(t)$ in a neighborhood of $z_{0}$ for some maximal time $T>0$, we define the bounded and closed set $H_{0,M}^{1}(\Omega)$ for a given $M>0$ by

[TABLE]

We need to show that the operator that describes the right-hand side of the flow problem (12) is Lipschitz-continuous on $H_{0,M}^{1}(\Omega)$. For the ease of notation we set $\mathcal{K}z:=\mathcal{G}_{z}z$.

Lemma 3.4**.**

The operator $\mathcal{K}:H^{1}_{0}(\Omega)\rightarrow H^{1}_{0}(\Omega)$ is Lipschitz-continuous and bounded on $H_{0,M}^{1}(\Omega)$. In particular, there exists a constant $L_{M}>0$ that depends on $M$, $z_{0}$, $\beta$, $\Omega$ and $d$, such that

[TABLE]

Proof.

Let $v,w\in H_{0,M}^{1}(\Omega)$ and set $M_{0}:=M+\|z_{0}\|_{H^{1}(\Omega)}$. Since

[TABLE]

we conclude that

[TABLE]

where $C$ only depends on the Poincaré-Friedrichs constant. With this, we have

[TABLE]

Using the Hölder inequality and embedding estimates, the second term on the right-hand side can be bounded by

[TABLE]

Using (13) and $v,w\in H_{0,M}^{1}(\Omega)$ we conclude the existence of some $L_{M}=\mathcal{O}(1+\beta M_{0}^{2})$ such that the Lipschitz-continuity holds true. ∎

The Lipschitz-continuity of Lemma 3.4 and the trivial observation that $\mathcal{K}v=0$ if and only if $v=0$ imply that there exists a sufficiently small neighbourhood $H_{0,M}^{1}(\Omega)$ of $z_{0}$ and a constant $c_{M}>0$ such that

[TABLE]

In such a neighborhood we have that the normalization factor $\gamma_{v}$, as a function in $v$, is also Lipschitz-continuous.

Lemma 3.5**.**

For any sufficiently small $M>0$ the functional $\gamma_{\cdot}:H^{1}_{0}(\Omega)\rightarrow\mathbb{R}$ is Lipschitz-continuous and bounded on $H_{0,M}^{1}(\Omega)$. In particular, there is $\tilde{L}_{M}>0$ (depending only on $M$, $z_{0}$, $\beta$, $V$, $\Omega$ and $d$) such that

[TABLE]

Proof.

The error in the difference of $\gamma_{v}$ and $\gamma_{w}$ can be expressed as

[TABLE]

and

[TABLE]

Using the splitting

[TABLE]

and the norm inequality $\|a\|^{2}-\|b\|^{2}\leq\|a-b\|(\|a\|+\|b\|)$ (which follows from the inverse triangle inequality) we see that there exists a constant $\tilde{C}_{M}>0$ such that for all $v,w\in H_{0,M}^{1}(\Omega)$ it holds

[TABLE]

The Lipschitz-continuity of $\mathcal{K}$ on $H_{0,M}^{1}(\Omega)$ as shown in Lemma 3.4 finishes the proof. ∎

The combination of Lemmas 3.4 and 3.5 shows that the right-hand side $g(z):=-z+\gamma_{z}\mathcal{G}_{z}z$ of the gradient flow problem of Definition 4.1 is Lipschitz-continuous in the neighborhood $H_{0,M}^{1}(\Omega)$ of $z_{0}$. Thus, the classical Picard-Lindelöf Theorem for Hilbert spaces implies local existence and uniqueness for some time $T>0$.

Lemma 3.6** (Local well-posedness).**

For any $z_{0}\in H^{1}_{0}(\Omega)$ with $\|z_{0}\|_{L^{2}(\Omega)}=1$, there exists a maximum time $T>0$ such that there is a unique solution $z$ to (12) on the time interval $[0,T)$.

3.2 Global well-posedness

Starting from the local existence of $z(t)$ guaranteed by Lemma 3.6, Lemma 3.3 allows us to pass to a global existence result. Note that as soon as such a global existence result is established, Lemma 3.3 implies mass conservation and energy-reduction for all times $t\in[0,\infty)$.

Proof of Theorem 3.2 - Global well-posedness.

Recall that Lemma 3.6 guarantees the existence of a unique solution on a time interval $[0,T)$ and assume that $T$ is finite and maximal in the sense that the problem is no longer well-posed for $t\geq T$. The energy reduction in Lemma 3.3 guarantees that $E_{T}:=\lim_{t\rightarrow T}E(z(t))$ exists. Next, let $t_{1},t_{2}\in[0,T)$ be arbitrary with $t_{1}\leq t_{2}$. Then, using the estimate $\|\int_{t_{1}}^{t_{2}}v(t)\hskip 2.0ptdt\|_{H^{1}(\Omega)}\leq\int_{t_{1}}^{t_{2}}\|v(t)\|_{H^{1}(\Omega)}\hskip 2.0ptdt$ (cf. [30, Appendix E.5]) and the construction of $z$ we see that

[TABLE]

This implies boundedness of $z(t)$ in $H^{1}_{0}(\Omega)$. Hence, we have the existence of a sequence $\{t^{n}\}_{n\in\mathbb{N}}$ with $t^{n}\rightarrow T$ and a function $z_{T}\in H^{1}_{0}(\Omega)$ so that $z(t^{n})\rightharpoonup z_{T}$ weakly in $H^{1}_{0}(\Omega)$. This implies

[TABLE]

Consequently

[TABLE]

Since the Hilbert space $H^{1}_{0}(\Omega)$ is uniformly convex, the weak convergence together with $\lim_{n\rightarrow{\infty}}\|z(t^{n})\|_{H^{1}(\Omega)}\leq\|z_{T}\|_{H^{1}(\Omega)}$ guarantee that $z(t^{n})\rightarrow z_{T}$ strongly in $H^{1}_{0}(\Omega)$. The continuity of $z$ in $t$ implies independence of this limit on the choice of the sequence $t^{n}$. Consequently, we have $z(t)\rightarrow z(T):=z_{T}$ for $t\rightarrow T$, strongly in $H^{1}_{0}(\Omega)$. Since we assumed $T<\infty$, we could use $z(T)$ as a new starting value to guarantee existence of $z$ on an extended interval $[0,T+\delta)$ for some $\delta>0$. This contradicts the assumed maximality of $T$ and, hence, shows that the problem admits a unique solution for all times. ∎

3.3 Global convergence and exponential decay to the ground state

With the previous results, we can now prove global $H^{1}$-convergence of $z(t)$ to a critical point of $E$ that fulfills the normalization constraint. In a first step, we need to make an identification of the limit.

Proof of Theorem 3.2 - Limit is eigenfunction of the GPE.

It remains to identify the limit as an eigenfunction of the GPE problem. Since the (non-negative) energy is decreasing along the flow there exists some limit $E_{\infty}:=\lim_{t\rightarrow{\infty}}E(z(t))$. With $a_{z(t)}(z^{\prime}(t),z^{\prime}(t))=-\frac{\mbox{\scriptsize d}}{\mbox{\scriptsize d}t}E(\hskip 1.0ptz(t)\hskip 1.0pt)$ we can conclude that

[TABLE]

This implies that $\int_{0}^{\infty}\|z^{\prime}\|^{2}_{H^{1}(\Omega)}\operatorname{\,d}t$ is finite and there exists $z^{\ast}\in H^{1}_{0}(\Omega)$ such that $z(t)\rightarrow z^{\ast}$ strongly in $H^{1}(\Omega)$. Obviously, the limit fulfills $z^{\ast}=\gamma_{z^{\ast}}\mathcal{G}_{z^{\ast}}z^{\ast}$ and, hence,

[TABLE]

i.e, $(z^{\ast},\gamma_{z^{\ast}})$ is an eigensolution of (2). ∎

We are now ready to prove the exponential convergence to the ground state.

Proof of Theorem 3.2 - Exponential convergence to ground state.

Assume that the strong $H^{1}$-limit of the flow $z(t)$ coincides with the unique positive ground state, i.e. $z^{\ast}=z_{\mbox{\rm\tiny GS}}$, where $z_{\mbox{\rm\tiny GS}}>0$ is characterized as in Proposition 2.1. From the first part of the proof of Theorem 3.2 we also know that $\gamma_{z}$ converges to the ground state eigenvalue $\lambda^{\ast}=\lambda_{\mbox{\rm\tiny GS}}>0$ for $t\rightarrow\infty$.

The proof of exponential convergence is based on a Grönwall-type argument. For that, we define the function

[TABLE]

and want to show that $f^{\prime}\leq c\hskip 2.0ptf$ for some positive constant $c$. Since $2\int_{0}^{\infty}f(t)\operatorname{\,d}t=E(z_{0})-E_{\infty}$ is finite, we know that $f(t)\rightarrow 0$ for $t\rightarrow\infty$ and hence $f(t)^{3/2}\leq f(t)$ for all sufficiently large times. Using this fact, we can conclude that for any $\varepsilon_{0}>0$, there exists a finite time $t(\varepsilon_{0})\geq 0$ such that for all $t\geq t(\varepsilon_{0})$ it holds

[TABLE]

Next, recall the projection $P_{z}$ from (11). With $P_{z}(z)=z-\gamma_{z}\mathcal{G}_{z}z=-z^{\prime}$, we can rewrite $f=\tfrac{1}{2}a_{z}(P_{z}(z),P_{z}(z))$ and, hence,

[TABLE]

where $\frac{\mbox{d}}{\mbox{d}t}P_{z}(z)$ is given by

[TABLE]

and where $a_{z}^{\prime}\langle v,(w_{1},w_{2})\rangle$ is the Fréchet derivative of $a_{z}(\cdot,\cdot)$ wrt. $z$, which can be computed as

[TABLE]

Consequently, we have with (15)

[TABLE]

where we used the Sobolev embedding $H^{1}_{0}(\Omega)\hookrightarrow L^{6}(\Omega)$ (for $d\leq 3$) in the last step. Next, we investigate the middle term. We use

[TABLE]

to see that

[TABLE]

The latter term can be written as

[TABLE]

where the first two terms are of higher order (due to the strong $H^{1}$-convergence to $z_{\mbox{\rm\tiny GS}}$) and the last term is strictly positive. Hence, for any $\varepsilon_{1}>0$ and sufficiently large times $t\geq t(\varepsilon_{1})$, we have the crude estimate

[TABLE]

The crucial estimate is now for $\|z^{\prime}\|_{L^{2}(\Omega)}^{2}$. First, we note that

[TABLE]

for $t\rightarrow\infty$, which shows that $z^{\prime}(t)$ converges strongly in $H^{1}$ to zero. Furthermore, we have

[TABLE]

where $C_{\Omega}$ is the Poincaré-Friedrichs constant on $\Omega$. We want to derive a sharper estimate for

[TABLE]

Assume that $(t_{n})_{n\in\mathbb{N}}$ with $t_{n}\rightarrow\infty$ is a corresponding minimal sequence so that $C_{\inf}$ is reached. In this case we have

[TABLE]

Note that $z^{\prime}(t_{n})$ is also bounded in $H^{1}(\Omega)$ and hence we can assume without loss of generality that there exists a weak limit $\hat{z}\in H^{1}_{0}(\Omega)$ with $\|\hat{z}\|_{L^{2}(\Omega)}=1$ such that for $n\rightarrow\infty$

[TABLE]

Together with the strong $H^{1}$-convergence of $z(t)$ to the ground state, this implies for $n\rightarrow\infty$

[TABLE]

Hence, the function $\hat{z}$ is orthogonal to the ground state $z_{\mbox{\rm\tiny GS}}$ both with respect to the $L^{2}$- and the $a_{z_{\mbox{\rm\tiny GS}}}(\cdot,\cdot)$-inner product. We conclude (with the lower semi-continuity of weakly converging sequences) that

[TABLE]

where $\mbox{span}\{z_{\mbox{\rm\tiny GS}}\}^{\perp}$ is the $a_{z_{\mbox{\rm\tiny GS}}}(\cdot,\cdot)$-orthogonal complement of the first eigenspace. Hence, with the Courant-Fischer Theorem we have

[TABLE]

where $\mu_{2}>\mu_{1}:=\lambda_{\mbox{\rm\tiny GS}}$ is the second eigenvalue of linear eigenvalue problem: find $w_{i}\in H^{1}_{0}(\Omega)$ with $\|w_{i}\|_{L^{2}(\Omega)}=1$ and $\mu_{i}\in\mathbb{R}$

[TABLE]

Here we exploited that $\mu_{1}=\lambda_{\mbox{\rm\tiny GS}}$ is the smallest eigenvalue of the linearized problem and that it is also simple (cf. [21, Lemma 2]). We can summarize that

[TABLE]

and hence, for all $\varepsilon_{2}>0$ there is a sufficiently large time $t\geq t(\varepsilon_{2})$ such that

[TABLE]

Here we used that $\gamma_{z}\rightarrow\lambda_{\mbox{\rm\tiny GS}}$. Combining (3.3), (14), (18) and (21) we have

[TABLE]

Selecting $\varepsilon_{0}$, $\varepsilon_{1}$ and $\varepsilon_{3}$ sufficiently small and the corresponding times sufficiently large, we see that for any $\delta>0$ there exists a finite time $t_{\delta}$ such that for all $t\geq t_{\delta}$

[TABLE]

By Grönwall’s lemma we obtain

[TABLE]

Hence, for every $\delta>0$ there exists a constant $c_{\delta}$ and a finite time $t_{\delta}$, such that

[TABLE]

Finally, we obtain for $0<\delta<1$ that

[TABLE]

This finishes the proof. ∎

4 Discrete Projected $a_{z}$-Sobolev Gradient Flow

In this section we propose and analyze a forward Euler discretization of the projected $a_{z}$-Sobolev gradient flow from Definition 3.1. For this purpose, let $\{\tau_{n}\}_{n\in\mathbb{N}}$ be a sequence of positive time steps that is bounded from above and below by

[TABLE]

The time steps can be seen as parameters that should be selected sufficiently large for the sake of computational efficiency. In the following, we use the simplifying notation and write

[TABLE]

With this, we consider the following forward Euler discretization of the continuous $a_{z}$-gradient flow.

Definition 4.1** (Method: GF$a_{z}$).**

Let $z^{0}\in H^{1}_{0}(\Omega)$ be given with $\|z^{0}\|_{L^{2}(\Omega)}=1$. Then for $n\geq 0$ the GF$a_{z}$-iteration $z^{n+1}\in H^{1}_{0}(\Omega)$ is defined as

[TABLE]

Since $z^{n}\in H^{1}_{0}(\Omega)$ and $\mathcal{G}^{n}z^{n}\in H^{1}_{0}(\Omega)$, the iterates are well-defined.

Remark 4.2** (Nonlinear inverse iteration).**

For the particular choice $\tau_{n}=1$ the iteration can be rewritten as

[TABLE]

which is the simplest form of the nonlinear inverse iteration (inverse power method). In this sense, GF$a_{z}$ is a generalized inverse iteration.

We emphasize that a (near) optimal $\tau_{n}$ can be cheaply computed by (nearly) minimizing the energy $E(z^{n+1})$ as a function of $\tau_{n}$ along the given search direction as described in the following remark.

Remark 4.3** (Adaptive GF$a_{z}$).**

The proposed method GF$a_{z}$ can be easily combined with an adaptive step size control to compute optimal values for $\tau_{n}$ in each step. This involves the minimization of the function

[TABLE]

w.r.t. $\tau_{n}\in(0,2)$. This can be done efficiently. Let us define

[TABLE]

and

[TABLE]

and also

[TABLE]

The terms $\alpha_{i}$, $\beta_{i}$ and $\zeta_{i}$ have to be precomputed only once per time step (with a single grid walk). With these terms and the function

[TABLE]

we can see that $f(\tau_{n})$ is given by

[TABLE]

This quantity can be evaluated cheaply once that $\alpha_{i}$, $\beta_{i}$ and $\zeta_{i}$ were precomputed. The minimization of $f(\tau_{n})$ on $(0,2)$ using e.g. golden section search leads to the (approximate) minimum $f(\tau_{n}^{\ast})$. The energy of $z^{n+1}$ is then given by $f(\tau_{n}^{\ast})$. Note that even without adaptivity, the quantity $f(\tau_{n})$ has to be computed, which is of the same order of complexity as the preprocessing step in the adaptive version. Hence, the computational overhead for using adaptivity is negligible. In particular, no additional linear system needs to be solved when using adaptivity with GF$a_{z}$. In contrast we observe that GF$L^{2}$ can typically not be efficiently combined with adaptivity, however, this also not as crucial as for the other methods since any sufficiently large choice for $\tau_{n}=\tau$ yields automatically a nearly optimal number of iterations for GF$L^{2}$ in terms of $\tau$.

The remaining parts of this section are devoted to the numerical analysis of this scheme.

4.1 Intermediate mass growth

While conservation of mass is guaranteed by normalization in each step of the iteration, it is worth studying the change of mass that is associated with the map $z^{n}\mapsto\hat{z}^{n+1}$. It will turn out that mass cannot be diminished under this operation. Multiplying equation (22) with $z^{n}$ and integrating over $\Omega$ yields

[TABLE]

This implies

[TABLE]

Hence $1=\|z^{n}\|_{L^{2}(\Omega)}^{2}\leq\|z^{n}\|_{L^{2}(\Omega)}^{2}+\|\hat{z}^{n+1}-z^{n}\|_{L^{2}(\Omega)}^{2}=\|\hat{z}^{n+1}\|_{L^{2}(\Omega)}^{2}$. We summarize it:

Lemma 4.4** (Intermediate mass growth).**

For all $n$ it holds $1=\|z^{n}\|_{L^{2}(\Omega)}\leq\|\hat{z}^{n+1}\|_{L^{2}(\Omega)}$. Furthermore, the normalization error can be expressed as

[TABLE]

The previous lemma implies that the normalization of $\hat{z}^{n}$ to unit mass necessarily decreases the energy. Moreover, if the mass is not increased, i.e. $\|\hat{z}^{n+1}\|_{L^{2}(\Omega)}=\|z^{n}\|_{L^{2}(\Omega)}$, then this implies that $z^{n}=\hat{z}^{n+1}=z^{n+1}$.

4.2 Energy dissipation

The proof of energy reduction is established in several steps. First, using the result from the previous subsection, applying the energy inner product $a_{z^{n}}(\cdot,\cdot)$ to (22) and testing with $\hat{z}^{n+1}-z^{n}$ yields

[TABLE]

This leads to a preliminary lower bound for the energy difference.

Lemma 4.5** (Sharp lower bounds for the energy difference).**

If $\tau_{n}\leq 2/5$ then

[TABLE]

For $\tau_{n}\geq 2$, either $E(\hat{z}^{n+1})>E(z^{n+1})$ or $z^{n}$ is already a critical point.

Proof.

Set $\tau:=\tau_{n}$. We get

[TABLE]

which implies

[TABLE]

Observe that

[TABLE]

and

[TABLE]

Combining everything yields

[TABLE]

For $\tau\geq 2$, the right hand side is negative which implies a growth of energy. It only remains to find a lower bound for the first term. Here we estimate

[TABLE]

which yields

[TABLE]

where we used $\tau\leq 2/5$. This finishes the proof. ∎

Remark 4.6** (Adaptive time steps).**

Observe that if the time steps $\tau_{n}$ are chosen adaptively, then asymptotically any choice $\tau_{n}<2$ is admissible. This however requires that the previous time steps (with typically smaller step size) were such that the iterates $z^{n}$ are in a sufficiently small neighbourhood of a critical point. This is because in the convergent regime, the first term in (25) (which is of order four) is eventually negligible, compared to the dominant second term which is only of second order. We will not exploit this observation, but believe that it is worth mentioning.

With Lemma 4.5 we can now prove energy reduction for sufficiently small time steps.

Lemma 4.7** (Energy reduction).**

There exists $0<\tau_{\mbox{\rm\tiny max}}<2$ (which depends on $\beta$ and $E(z^{0})$) such that for all $\tau_{n}\leq\tau_{\mbox{\rm\tiny max}}$

[TABLE]

Remark 4.8** (Energy reduction for $\tau_{n}=1$).**

Numerically, we could observe the coupling between $\tau_{\mbox{\rm\tiny max}}$ and the energy of $z^{0}$ at several occasions, i.e. if $E(z_{0})$ was large then the step size $\tau_{n}$ had to be reduced to obtain reduction of the energy. However, we never observed that $\tau_{\mbox{\rm\tiny max}}$ dropped below one. In this connection we shall note that, for $\tau_{n}=1$, the GF$a_{z}$ can be seen as a GF$L^{2}$ realization applied to the Schrödinger operator whose spectrum was shifted by $+1$. Consider the GPE with the modified potential $V_{\mbox{\tiny\rm mod}}:=V-1$. Applying GF$L^{2}$ to this modified problem gives the same iterations as applying the GF$a_{z}$ iterations to the GPE with original potential (for the particular choice $\tau_{n}=1$). Hence, both methods produce the same approximations $z^{n}$. Using the results obtained in [12, Lemma 2.10] for GF$L^{2}$ with $\tau_{n}=1$ we can hence argue that the GF$a_{z}$ iterates are guaranteed to reduce a functional of the form $\tilde{E}(v)=E(v)+\tfrac{\beta}{4}\|v\|_{L^{4}(\Omega)}^{4}$. Since $\lambda=2E(u)+\tfrac{\beta}{2}\|u\|_{L^{4}(\Omega)}^{4}=2\tilde{E}(v)$, this can be seen as minimizing an “eigenvalue functional” instead of the original energy functional.

Proof of Lemma 4.7 (by induction).

Step $n=0$: With (4.2) we have

[TABLE]

If $\tau_{0}^{2}\leq(4E(z^{0}))^{-1}$, we have $a_{0}(\hat{z}^{1}-z^{0},\hat{z}^{1}-z^{0})\leq 1$ and hence

[TABLE]

Together with (25), we conclude

[TABLE]

Hence, there exists $\tau_{\mbox{\rm\tiny max}}\lesssim\min\{\beta^{-1},E(z^{0})^{-1/2}\}$ such that for all $\tau_{0}\leq\tau_{\mbox{\rm\tiny max}}$ we have

[TABLE]

Step $n\mapsto n+1$: Let $E(z^{n})\leq E(z^{0})$ and $\tau_{n}\leq\tau_{\mbox{\rm\tiny max}}$ with $\tau_{\mbox{\rm\tiny max}}$ as for $n=0$. Using (4.2) and $E(z^{n})\leq E(z^{0})$ we have

[TABLE]

Analogously as for $n=0$, we have

[TABLE]

and, hence,

[TABLE]

Note that $c_{\tau}\rightarrow\infty$ for $\tau\rightarrow 0$. ∎

4.3 Global convergence

We have the following main result on the global convergence of the discrete gradient flow.

Theorem 4.9**.**

We consider the GF$a_{z}$-approach stated in Definition 4.1. Assume that the time steps fulfill $\tau_{n}\leq\tau_{\mbox{\rm\tiny max}}$ as in Lemma 4.7 and that they are non-degenerate in the sense that $\tau_{n}\geq\tau_{\mbox{\rm\tiny min}}>0$. Then there exists a limit energy $E^{\ast}:=\lim_{n\rightarrow\infty}E(z^{n})$. Furthermore, there exists a subsequence $\{z^{n_{i}}\}_{i\in\mathbb{N}}$ of $\{z^{n}\}_{n\in\mathbb{N}}$, such that $z^{n_{i}}\rightarrow z^{\ast}$ strongly in $H^{1}_{0}(\Omega)$ to some limit $z^{\ast}\in H^{1}_{0}(\Omega)$ with $\|z^{\ast}\|_{L^{2}(\Omega)}=1$ and $E(z^{\ast})=E^{\ast}$. With

[TABLE]

we have that $z^{\ast}$ is an eigenfunction to the Gross-Pitaevskii equation and fulfills

[TABLE]

Any other convergent subsequence of $\{z^{n}\}_{n\in\mathbb{N}}$ will also converge strongly in $H^{1}(\Omega)$ to an $L^{2}$-normalized eigenfunction of the GPE with energy level $E^{\ast}$. However, the corresponding eigenvalue might be different from $\lambda^{\ast}$ above.

Remark 4.10**.**

Note that $\mathcal{G}^{\ast}(z^{\ast})=\frac{1}{\lambda^{\ast}}z^{\ast}$.

Recall that the limit energy in Theorem 4.9 depends crucially on $z^{0}$, but potentially it can also depend on the choice of the sequence $\{\tau_{n}\}_{n\in\mathbb{N}}$. If there exists only one eigenfunction (up to normalization and multiplication with $-1$) for the energy level $E^{\ast}$, then we have convergence of the full sequence in Theorem 4.9.

Proof of Theorem 4.9.

As $E(z^{n})$ is a monotonically decreasing sequence the limit $E^{\ast}:=\lim_{n\rightarrow\infty}E(z^{n})$ exists. This means that $\{z^{n}\}_{n\in\mathbb{N}}$ is a bounded sequence in $H^{1}_{0}(\Omega)$ from which we can extract a subsequence, for brevity still denoted by $\{z^{n}\}_{n\in\mathbb{N}}$, that converges weakly in $H^{1}_{0}(\Omega)$ to some limit function $z^{\ast}\in H^{1}_{0}(\Omega)$ with $\|z^{\ast}\|_{L^{2}(\Omega)}=1$. For space dimension $d\leq 3$ the Rellich-Kondrachov theorem guarantees that $z^{n}$ converges to $z^{\ast}$, strongly in $L^{p}(\Omega)$ for $p<6$. First, we note that (26) implies that $\|\hat{z}^{n+1}-z^{n}\|_{L^{2}(\Omega)}\rightarrow 0$ and consequently

[TABLE]

Furthermore, it holds for any $v\in H^{1}_{0}(\Omega)$

[TABLE]

Using the aforementioned Rellich-Kondrachov embedding we have that $|z^{n}|^{2}\rightarrow|z^{\ast}|^{2}$ strongly in $L^{2}(\Omega)$ and hence

[TABLE]

The later equation implies that $\mathcal{G}^{n}z^{n}$ converges weakly in $H^{1}_{0}(\Omega)$ (and strongly in $L^{2}(\Omega)$) to $\mathcal{G}^{\ast}z^{\ast}$. Combining the strong $L^{2}$-convergence of $\mathcal{G}^{n}z^{n}$ and $z^{n}$ we obtain

[TABLE]

Combining all the results we can use

[TABLE]

and pass to the limit for any $v\in H^{1}_{0}(\Omega)$ in

[TABLE]

Thus, $a_{z^{\ast}}(z^{\ast},v)=\gamma^{\ast}(z^{\ast},v)_{L^{2}(\Omega)}$ for all $v\in H^{1}_{0}(\Omega)$. To verify the convergence of the energy, i.e. $E^{\ast}=E(z^{\ast})$, observe that

[TABLE]

Using this expression, we have

[TABLE]

For all terms on the right-hand side we verified (strong) convergence. Consequently we have

[TABLE]

The strong convergence of $z^{n}$ in $H^{1}(\Omega)$ follows readily from the previous result as it implies $\lim_{n\rightarrow\infty}\|z^{n}\|_{H^{1}(\Omega)}=\|z^{\ast}\|_{H^{1}(\Omega)}$. ∎

It is easily seen that all proofs in this section remain valid, if we replace the space $H^{1}_{0}(\Omega)$ in the GF$a_{z}$-approach by a finite dimensional subspace, e.g. in a spatial finite element discretization.

Corollary 4.11** (Convergence of the fully discrete GF$a_{z}$).**

Let $V_{h}\subset H^{1}_{0}(\Omega)$ be a finite dimensional subspace and let $\mathcal{G}_{z}^{h}(z_{h})\in V_{h}$ solve $a_{z}(\mathcal{G}_{z}^{h}(z_{h}),v_{h})=(z_{h},v_{h})_{L^{2}(\Omega)}$ for all $v_{h}\in V_{h}$. For $z_{h}^{0}\in V_{h}$ with $\|z_{h}^{0}\|_{L^{2}(\Omega)}=1$ we consider the GF$a_{z}$ iteration

[TABLE]

If $\tau\leq\tau_{\mbox{\rm\tiny max}}$ then the energy is strictly reduced and there exists a limit energy $E^{\ast}_{h}:=\lim_{n\rightarrow\infty}E(z^{n}_{h})$. Furthermore, up to subsequences, we have $z^{n}_{h}\rightarrow z^{\ast}_{h}$ strongly in $H^{1}_{0}(\Omega)$ where $z^{\ast}_{h}\in V_{h}$ with $\|z^{\ast}_{h}\|_{L^{2}(\Omega)}=1$ and $E(z^{\ast}_{h})=E^{\ast}_{h}$ is a discrete eigenfunction of the GPE, i.e. there is $\lambda_{h}^{\ast}$ so that

[TABLE]

The convergence of approximate eigenpairs $(\lambda^{\ast}_{h},z^{\ast}_{h})$ to the true ones has been investigated and analyzed in [21, 35, 20].

5 Global convergence to the ground state

Theorem 4.9 shows uniqueness of the limit of the discrete flow $(z^{n})$ under uniqueness of the eigenfunction (up to normalization) on the energy level $E^{\ast}$. The latter assumption can be relaxed in the particular case of eigenstates that are strictly positive in the interior of $\Omega$. As we have already seen, there exists at least one such state, which is the ground state of the energy functional $E$. In this section we will prove that it is also the only one. This is crucial for the following main result.

Theorem 5.1**.**

Consider the GF$a_{z}$-approach. Let the assumptions of Theorem 4.9 hold and $\tau_{n}\leq 1$ for all $n$. Then for any starting value $z^{0}\in H^{1}_{0}(\Omega)$ with $\|z^{0}\|_{L^{2}(\Omega)}=1$ and $z^{0}\geq 0$ the (full) sequence $(z^{n})$ converges strongly in $H^{1}(\Omega)$ to the positive ground state $z_{\mbox{\rm\tiny GS}}$ (which is unique according to Proposition 2.1).

Remark 5.2**.**

Theorem 5.1 guarantees global convergence to the ground state, provided that the starting value $z^{0}$ is not changing its sign. Additionally, starting from a non-negative $z^{0}$, the gradient flow does not converge to an excited state. A sign-changing starting value is compulsory for the computation of excited states.

Before we can prove Theorem 5.1, a few auxiliary results are required. The first result relates the positive eigenfunctions in the spectrum of the GPE to the ground states of a linear operator obtained by freezing the density. The lemma can be proved analogously to a similar result obtained in [21, Lemma 2].

Lemma 5.3**.**

Let $z^{\ast}\in H^{1}_{0}(\Omega)$ with $\|z^{\ast}\|_{L^{2}(\Omega)}=1$ be an eigenstate of the GPE with eigenvalue $\lambda^{\ast}>0$, i.e.

[TABLE]

If $z^{\ast}\geq 0$ in $\Omega$, then $z^{\ast}$ can be characterized as the as the unique positive ($L^{2}$-normalized) ground state to the linear operator $\mathcal{G}_{z^{\ast}}^{-1}$ (see Remark 4.10) and $z^{\ast}$ must be even strictly positive in the interior of $\Omega$.

Next, we prove that the positive eigenstate is unique and hence always the ground state.

Lemma 5.4** (Uniqueness of positive eigenstates).**

There is a unique positive eigenfunction to the GPE (2), which is the ground state.

Proof.

We recall the Picone identity (cf. [41, 17]), which implies that for two functions $u,v\in H^{1}(\Omega)$ with $u\geq 0$ and $v>0$ in $\Omega$ it holds

[TABLE]

From Lemma 5.3 we know that any nonnegative eigenfunction must be even strictly positive. Let us therefore assume we have two positive $L^{2}$-normalized eigenfunctions $z_{\mbox{\rm\tiny GS}},z_{\mbox{\rm\tiny ES}}\in H^{1}_{0}(\Omega)$ to the Gross-Pitaevskii equation, where $z_{\mbox{\rm\tiny GS}}$ is the unique ground state with minimal energy $E_{\mbox{\rm\tiny GS}}$ and eigenvalue $\lambda_{\mbox{\rm\tiny GS}}$ and $z_{\mbox{\rm\tiny ES}}$ is an excited state with energy $E_{\mbox{\rm\tiny ES}}>E_{\mbox{\rm\tiny GS}}$. Using $\|z_{\mbox{\rm\tiny GS}}\|_{L^{2}(\Omega)}=1$ it holds

[TABLE]

We conclude that

[TABLE]

This is a contradiction to the assumption that $z_{\mbox{\rm\tiny ES}}$ was an excited state with $E_{\mbox{\rm\tiny ES}}>E_{\mbox{\rm\tiny GS}}$. Hence, we have $z_{\mbox{\rm\tiny ES}}=z_{\mbox{\rm\tiny GS}}$ which is unique. ∎

The next (fairly obvious) result shows that positivity is preserved by the iteration.

Lemma 5.5**.**

Let $v,z\in H^{1}_{0}(\Omega)$. Then

[TABLE]

In particular, if $z^{n}\geq 0$ and $\tau_{n}\leq 1$ then $z^{n+1}\geq 0$.

Proof.

We can characterize $\mathcal{G}_{z}(v)$ as the unique minimizer of

[TABLE]

among all $w\in H^{1}_{0}(\Omega)$. However, since it holds $F(|\mathcal{G}_{z}(v)|)\leq F(\mathcal{G}_{z}(v))$ we conclude by uniqueness $|\mathcal{G}_{z}(v)|=\mathcal{G}_{z}(v)$, which guarantees that $\mathcal{G}_{z}(v)$ cannot become negative. The positivity of $z^{n+1}$ follows immediately with $\hat{z}^{n+1}=(1-\tau_{n})z^{n}+\tau_{n}\gamma^{n}\mathcal{G}^{n}z^{n}$, where $\gamma^{n},z^{n},\mathcal{G}^{n}z^{n}\geq 0$ and $\tau_{n}\leq 1$. ∎

We are now ready to prove the main result of this section.

Proof of Theorem 5.1.

Let $z^{\ast}\in H^{1}_{0}(\Omega)$ be any of the limits of a subsequence of $z^{n}$ whose existence is guaranteed by Theorem 4.9 with $E^{\ast}=E(z^{\ast})$. Then we have that for any $z\in H^{1}_{0}(\Omega)$ with $\|z\|_{L^{2}(\Omega)}=\|z^{\ast}\|_{L^{2}(\Omega)}=1$ it holds

[TABLE]

Here $\lambda^{\ast}$ is the eigenvalue to the eigenfunction $z^{\ast}$. Since $z^{\ast}$ is the strong $H^{1}$-limit of a sequence of positive functions $z^{n_{i}}$, pointwise convergence almost everywhere ensures that $z^{\ast}\geq 0$. Hence, we can apply Lemma 5.3 that guarantees $z^{\ast}>0$ and that $\lambda^{\ast}>0$ is the ground state eigenvalue of the linear operator $\mathcal{G}_{z^{\ast}}^{-1}$. Hence, it holds $\langle\mathcal{G}_{z^{\ast}}^{-1}v,v\rangle\geq\lambda^{\ast}(v,v)_{L^{2}(\Omega)}$ for any $v\in H^{1}_{0}(\Omega)$ or respectively

[TABLE]

Using this finding in (LABEL:energy-identity-z-ast) implies

[TABLE]

where the global convergence of the energies is ensured by Theorem 4.9. Since $z^{n},z^{\ast}\geq 0$, we conclude convergence of the whole sequence $z^{n}$ to $z^{\ast}$. That means that all strong $H^{1}$-limits of subsequences in Theorem 4.9 must coincide. Lemma 5.4, the uniqueness of the nonnegative eigenstates, finishes the proof. ∎

Remark 5.6**.**

Elliptic regularity theory provides $H^{2}$- and $L^{\infty}$-bounds for $\mathcal{G}_{z}(v)$ which are of the form

[TABLE]

This implies that in the energy diminishing regime, the iterates $z^{n}$ remain pointwise uniformly bounded, with a bound that depends on $\beta$, $V$ and $E(z^{0})$.

6 Numerical experiments

This section concerns the numerical performance of the proposed projected $a_{z}$-Sobolev gradient flow GF$a_{z}$ defined in (22). For a better assessment, we compare with established gradient flows, the GF$L^{2}$ iteration (or DNGF) from (6), the $H^{1}$-Sobolev gradient version GF$H^{1}$ from (8) and the $a_{0}$-Sobolev gradient version GF$a_{0}$ from (10) that incorporates the potential $V$. For the sake of a fair comparison of all methods, we use the (otherwise impractical) stopping criterion that the relative error with respect to some highly accurate (accuracy order $10^{-8}$) reference energy falls below the tolerance $\mbox{\rm\small TOL}=10^{-5}$. For the sake of simplicity, we measure performance in terms of number of iterations required to match this stopping criterion. This is a reasonable complexity indicator because the computational cost per iteration is essentially the same for all methods if a uniform step size $\tau$ is used. While GF$L^{2}$ and GF$a_{z}$ require the assembly of a new stiffness matrix from the previous density $|z^{n}|^{2}$ and one linear solve, GF$a_{0}$ and GF$H_{1}$ require two solves but the system matrices are invariant and do not need to be re-assembled. Our practical experience is that GF$L^{2}$ and GF$a_{z}$ iterations are slightly faster than the other two but this will not be taken into account in the following comparison.

As a general model, we solve the following Gross-Pitaevskii eigenvalue problem: find $z^{\ast}\in H^{1}_{0}(\Omega)$ with $\|z^{\ast}\|_{L^{2}(\Omega)}=1$ and

[TABLE]

for all $v\in H^{1}_{0}(\Omega)$ and in a bounded domain $\Omega$ of $\mathbb{R}^{2}$. Note that the kinetic part, i.e. $(\nabla z^{\ast},\nabla v)_{L^{2}(\Omega)}$, has an additional scaling factor $1/2$ compared to previously considered problem (2). The particular choices of $\Omega$, $V$ and $\beta$ are specified separately in the various experiments. All problems are discretized using a $P1$-Lagrange finite element method on a uniform grid of width $h$ specified below. Although adaptivity (as explained in Remark 4.3) can be used to improve the performance of GF$a_{z}$ (and also GF$H_{1}$, GF$a_{0}$), our comparisons focus on equidistant steps $\tau$.

Remark 6.1**.**

We stress that our comparison only aims at comparing the basic versions of the gradient flow methods and that each of these methods can be improved significantly with various techniques and hence the overall picture might change in this case. Here, we refer for example to the improvements of GF$L^{2}$ by using preconditioners and conjugated gradients as suggested in [8] or the improvements of GF$a_{0}$ by using Riemannian conjugate gradients as proposed in [27]. Such improvement can boost the performance dramatically compared to the basic versions of the gradient flow methods (cf. the numerical experiments in [8, 27]). Furthermore, adaptive mesh refinement strategies can improve the efficiency even further [34]. Another strategy, which can be particularly beneficial for excited states, is to use a different linearization technique that is based on the derivative of a scaling-invariant version of the Gross-Pitaevskii operator and which reacts more favorably to spectral shifts [5, 36].

6.1 Model problem 1 - Ground states for a harmonic potential

In the first model problem, we consider (30) for a harmonic trapping potential with trapping frequencies $1/2$, i.e.

[TABLE]

The repulsion parameter $\beta$ is selected with three different values $\beta=10,100,1000$. Computing the corresponding Thomas-Fermi radii of the problem we restrict the computations to a square domain of the size $\Omega=(-6,6)^{2}$. The initial value $z^{0}$ is selected as the Thomas-Fermi density computed according to [9] using the exact ground state for $\beta=0$. Since this is a nonnegative initial value, we expect all numerical approximations to converge to the unique positive ground state of $E$ (if $\tau$ is in the convergent regime). The ground state energies and eigenvalues for different values of $\beta$ are listed in Table 1.

Throughout our numerical experiments we observed that the stability regions for GF$H^{1}$ and GF$a_{0}$ are notably smaller than the ones for GF$L^{2}$ and GF$a_{z}$. Furthermore, the size of the spatial mesh size $h$ has essentially no influence on the convergence and number of steps required to fall below the tolerance. Both of these findings become visible in the results depicted in Table 2 where we compare the different methods for the ad-hoc parameter choices $\tau=\tau_{n}=0.5$ and $\tau=\tau_{n}=1$ and for the mesh sizes $h=12\cdot 2^{-6}$ and $h=12\cdot 2^{-8}$. With the default choice $\tau=1$, GF$L^{2}$ and GF$a_{z}$ perform equally well. In general we observe that GF$a_{z}$ is more sensitive with respect to the step size parameter $\tau$.

Since the tables show only the results for two exemplary choices of $\tau$, it is more interesting to investigate what, for a fixed setup, is the minimum number of iterations that the methods require to reach the tolerance. We keep the step size $\tau$ constant. Corresponding results are depicted in Table 3 for the three different values of $\beta$. We observe that the bigger $\beta$, the more iterations are required, though the growth is only moderate. We see that GF$a_{z}$ requires the fewest iterations, closely followed by GF$L^{2}$. Both GF$H^{1}$ and GF$a_{0}$ perform decently, though they are considerably behind the other two approaches. We made the same observation in various experiments and assume that this is linked to the smaller stability domain of the GF$H^{1}$ and GF$a_{0}$, enforcing smaller values for $\tau$ and hence smaller updates in modulus. Optimal values for $\tau$ can be found by solving a minimization problem for $\tau$ in each time step (cf. [26, Section 4]). If this is not done, the GF$H^{1}$ approach can be tough to use, because a stable constant time step is rather small.

In Lemma 4.5 we observed the expected divergence (energy blow-up) for the GF$a_{z}$-approach for time steps $\tau\geq 2$. This bound seems to be pretty sharp according to further numerical experiments not presented here. All $H^{1}$ gradient flows share such a time step restriction. Only GF$L^{2}$ is unconditionally stable for all $\tau<\infty$.

The remaining experiments focus on a comparison between GF$a_{z}$ and GF$L^{2}$.

6.2 Model problem 2 - Ground state in a lattice potential

In the second model problem, again based on (30), we investigate how the GF$a_{z}$ and GF$L^{2}$ methods perform when using a more complicated potential $V$ which consists of a harmonic part and an additional optical lattice. The potential is visualized in Figure 1 (left) and reads

[TABLE]

Furthermore, we use again $\Omega=(-6,6)^{2}$ and $\beta=1000$ and fix the mesh size $h=12\cdot 2^{-8}$. To compute the ground state of the corresponding energy functional, we start the different iterations with a Thomas-Fermi density that was computed according to [9] (for the case of general potentials, which includes the lattice part in our case). The final ground state density is depicted in Figure 1 (right), where we identified the ground state energy with approximately $E_{\mbox{\rm\tiny GS}}=15.204825$ and the corresponding ground state eigenvalue with $\lambda_{\mbox{\rm\tiny GS}}=36.708$.

In Table 4(a) we see how the number of GF$a_{z}$ iterations vary depending on the selected step size $\tau$. The method is unstable for $\tau\leq 1.7$. For smaller time steps, the number of iterations decreases uniformly from 23 iterations for $\tau=0.8$ to $12$ iterations for $\tau=1.5$. Even though not contained in the table, the number of iterations for GF$a_{z}$ increases dramatically for $\tau\leq 0.7$ and the method is no longer competitive in this regime. In practice we always recommend the usage of adaptivity (cf. Remark 4.3) to find a good value for $\tau$. In this case only $11$ iterations were needed to achieve the error tolerance. The GF$L^{2}$ is less sensitive to the choice of time step. However, the minimal number of time steps $26$ is considerably higher to what is achieved by GF$a_{z}$. With the right choice of the step size, GF$a_{z}$ performs up to twice as fast. Using adaptivity the appropriate time step regime is easily reached.

Our general conclusion is that for simple test problems the GF$L^{2}$ and GF$a_{z}$ perform basically evenly. On the other hand, the GF$a_{z}$ can have visible advantages for more challenging test cases involving poor choices for the starting value $z^{0}$ or more complicated potentials.

Remark 6.2** (Negative potentials, shift and invert).**

Shifting the potential $V$ by $-20$ leads to a negative potential but does not affect the eigenfunctions. All energies and corresponding eigenvalues are simply shifted by $-20$ as well. The ground state energy level then reads $E_{\mbox{\rm\tiny GS}}=5.204825$ and the corresponding eigenvalue $\lambda_{\mbox{\rm\tiny GS}}=16.708$. Still, the negative potential causes problems for numerical simulation. We observed strong energy oscillations for the GF$L^{2}$ if the step size was not selected sufficiently small ($\tau<0.7$ in our tests). Such oscillations cannot happen if $V\geq 0$. Therefore it is reasonable to first shift $V$ so that it becomes positive, apply the methods to compute e.g. the ground state and afterwards shift the energy and the eigenvalue back to the original setup. This is equivalent to using a suitable shift parameter in a conventional inverse iteration method.

As with linear eigenvalue problems, such a shift may as well be used to speed up convergence by increasing the relative sizes of spectral gaps.

6.3 Model Problem 3 - Anderson Localization

Our final numerical experiments is devoted to the phenomenon of Anderson localization [7], which describes the exponential localization of waves in a disordered medium. In the context of the Gross-Pitaevskii eigenvalue problem this Anderson effect is reflected by strongly localized peaks in the ground state eigenfunction, provided that the potential $V$ is sufficiently disordered.

We consider (30) and let $\Omega=(-6,6)^{2}$ and $\beta=10$. The potential $V$ is a random disorder potential that is obtained by dividing $\Omega$ into $400\times 400$ square cells with edge length $\varepsilon=0.03$. In each cell independently, the potential takes either the value $V(x)=1$ or $V(x)=\varepsilon^{-2}$ with equal probability. The scaling is selected according to the theoretical findings in [6]. The particular (deterministic) realization of $V$ used in our experiment is depicted in Figure 2, together with the corresponding ground state $z_{\mbox{\rm\tiny GS}}$. We can clearly see the expected Anderson localization, as $z_{\mbox{\rm\tiny GS}}$ consists of few exponentially fast decaying peaks and is essentially zero elsewhere. With a highly accurate reference computation we obtained the ground state energy with $E_{\mbox{\rm\tiny GS}}=4.84223025$ and the ground state eigenvalue with $\lambda_{\mbox{\rm\tiny GS}}=10.826242$. The uniform mesh in our computations has the mesh size $h=12\cdot 2^{-8}$ which is fine enough to resolve the variations of the potential. The initial value was again selected as a suitable Thomas-Fermi approximation.

In Table 5 we can see the number of iterations for GF$a_{z}$ and GF$L^{2}$. Again, we observe a similar performance of both methods, where GF$a_{z}$ shows stronger variations in the number of iterations. However, comparing the peak performance of the approaches, we see that GF$L^{2}$ is around 27% slower than GF$a_{z}$. It is interesting to note that we observed convergence of GF$a_{z}$ until very close to the theoretical upper limit of $\tau=2$. Combining GF$a_{z}$ with an adaptive step size control as described at the beginning of this section, the number of iterations dropped even further from $56$ to $52$. In general we can conclude that both GF$L^{2}$ and GF$a_{z}$ are well-suited for an efficient computation of Anderson localized ground states, where the GF$a_{z}$ with adaptivity shows clearly the best performance.

Acknowledgements. The authors thank Robert Altmann for the fruitful discussions and valuable comments on some of the proofs. Furthermore, we thank the anonymous reviewers for their very insightful comments that greatly improved the contents of this paper.