Computing Ground States of Spin-1 Bose-Einstein Condensates by the Normalized Gradient Flow

# Computing Ground States of Spin-1 Bose-Einstein Condensates by the Normalized Gradient Flow

Weizhu Bao Department of Mathematics and Center for Computational Science and Engineering, National University of Singapore, Singapore 117543 (bao@math.nus.edu.sg    Fong Yin Lim . Department of Mathematics and Center for Computational Science and Engineering, National University of Singapore, Singapore 117543 (fongyin.lim@nus.edu.sg)
###### Abstract

In this paper, we propose an efficient and accurate numerical method for computing the ground state of spin-1 Bose-Einstein condensates (BEC) by using the normalized gradient flow or imaginary time method. The key idea is to find a third projection or normalization condition based on the relation between the chemical potentials so that the three projection parameters used in the projection step of the normalized gradient flow are uniquely determined by this condition as well as the other two physical conditions given by the conservation of total mass and total magnetization. This allows us to successfully extend the most popular and powerful normalized gradient flow or imaginary time method for computing the ground state of single component BEC to compute the ground state of spin-1 BEC. An efficient and accurate discretization scheme, the backward-forward Euler sine-pseudospectral method (BFSP), is proposed to discretize the normalized gradient flow. Extensive numerical results on ground states of spin-1 BEC with ferromagnetic/antiferromagnetic interaction and harmonic/optical lattice potential in one/three dimensions are reported to demonstrate the efficiency of our new numerical method.

Key words. Spin-1 Bose-Einstein condensate, coupled Gross-Pitaevskii equations, ground state, normalized gradient flow, backward-forward Euler sine-pseudospectral method

AMS subject classifications. 35Q55, 65T99, 65Z05, 65N12, 65N35, 81-08

## 1 Introduction

Research in low temperature dilute atomic quantum gases remains active for more than ten years after the experimental realizations of Bose-Einstein condensation (BEC) in alkali atomic gases in 1995 [2, 12, 17]. Extensive theoretical and experimental studies have been carried out to investigate various novel phenomena of the condensates. In earlier BEC experiments, the atoms were confined in magnetic trap [2, 12, 17], in which the spin degrees of freedom is frozen. The particles are described by a scalar model and the wavefunction of the particles is governed by the Gross-Pitaevskii equation (GPE) within the mean-field approximation [16, 11, 23]. In recent years, experimental achievement of spin-1 and spin-2 condensates [10, 18, 21, 26, 28] offers new regimes to study various quantum phenomena that are generally absent in a single component condensate. The spinor condensate is achieved experimentally when an optical trap, instead of a magnetic trap, is used to provide equal confinement for all hyperfine states.

The theoretical studies of spinor condensate have been carried out in several papers since the achievement of it in experiments [19, 20, 22, 27]. In contrast to single component condensate, a spin- () condensate is described by a generalized coupled GPEs which consists of equations, each governing one of the hyperfine states within the mean-field approximation. For a spin-1 condensate, at temperature much lower than the critical temperature , the three-components wavefunction are well described by the following coupled GPEs [27, 29, 30, 31],

 iℏ∂tψ1(x,t) = [−ℏ22m∇2+V(x)+(c0+c2)(|ψ1|2+|ψ0|2)+(c0−c2)|ψ−1|2]ψ1 \hb@xt@.01(1.1) +c2¯ψ−1ψ20, iℏ∂tψ0(x,t) = [−ℏ22m∇2+V(x)+(c0+c2)(|ψ1|2+|ψ−1|2)+c0|ψ0|2]ψ0 \hb@xt@.01(1.2) +2c2ψ−1¯ψ0ψ1, iℏ∂tψ−1(x,t) = [−ℏ22m∇2+V(x)+(c0+c2)(|ψ−1|2+|ψ0|2)+(c0−c2)|ψ1|2]ψ−1 \hb@xt@.01(1.3) +c2ψ20¯ψ1.

Here, is the Cartesian coordinate vector, is time, is the Planck constant, is the atomic mass, and is the external trapping potential. When a harmonic trap potential is considered,

 V(x)=m2(ω2xx2+ω2yy2+ω2zz2), \hb@xt@.01(1.4)

with , and being the trap frequencies in the -, - and -direction, respectively. and Re denote the conjugate and real part of the function , respectively. There are two atomic collision terms, and , expressed in terms of the -wave scattering lengths, and , for scattering channel of total hyperfine spin 0 (anti-parallel spin collision) and spin 2 (parallel spin collision), respectively. The usual mean-field interaction, , is positive for repulsive interaction and negative for attractive interaction. The spin-exchange interaction, , is positive for antiferromagnetic interaction and negative for ferromagnetic interaction. The wave function is normalized according to

 ∥Ψ∥2:=∫R3|Ψ(x,t)|2dx=∫R31∑l=−1|ψl(x,t)|2dx:=1∑l=−1∥ψl∥2=N, \hb@xt@.01(1.5)

where is the total number of particles in the condensate.

By introducing the dimensionless variables: with , with , (), we get the dimensionless coupled GPEs from (LABEL:eq:CGPE_1)-(LABEL:eq:CGPE_3) as [30, 32, 9]:

 i∂tψ1(x,t) = [−12∇2+V(x)+(βn+βs)(|ψ1|2+|ψ0|2)+(βn−βs)|ψ−1|2]ψ1 \hb@xt@.01(1.6) +βs¯ψ−1ψ20, i∂tψ0(x,t) = [−12∇2+V(x)+(βn+βs)(|ψ1|2+|ψ−1|2)+βn|ψ0|2]ψ0 \hb@xt@.01(1.7) +2βsψ−1¯ψ0ψ1, i∂tψ−1(x,t) = [−12∇2+V(x)+(βn+βs)(|ψ−1|2+|ψ0|2)+(βn−βs)|ψ1|2]ψ−1 \hb@xt@.01(1.8) +βsψ20¯ψ1;

where , and with , and . Similar as those in single component BEC [24, 8, 3, 6], in a disk-shaped condensation, i.e. and (, and with ), the 3D coupled GPEs (LABEL:GPEs80)-(LABEL:GPEs82) can be reduced to a 2D coupled GPEs; and in a cigar-shaped condensation, i.e. and (, and with ), the 3D coupled GPEs (LABEL:GPEs80)-(LABEL:GPEs82) can be reduced to a 1D coupled GPEs. Thus here we consider the dimensionless coupled GPEs in -dimensions ():

 i∂tψ1(x,t) = [−12∇2+V(x)+(βn+βs)(|ψ1|2+|ψ0|2)+(βn−βs)|ψ−1|2]ψ1 \hb@xt@.01(1.9) +βs¯ψ−1ψ20, i∂tψ0(x,t) = [−12∇2+V(x)+(βn+βs)(|ψ1|2+|ψ−1|2)+βn|ψ0|2]ψ0 \hb@xt@.01(1.10) +2βsψ−1¯ψ0ψ1, i∂tψ−1(x,t) = [−12∇2+V(x)+(βn+βs)(|ψ−1|2+|ψ0|2)+(βn−βs)|ψ1|2]ψ−1 \hb@xt@.01(1.11) +βsψ20¯ψ1.

In the equations above, is a real-valued potential whose shape is determined by the type of system under investigation, and correspond to the dimensionless mean-field (spin-independent) and spin-exchange interaction, respectively. Three important invariants of (LABEL:eq:dCGPE_1)-(LABEL:eq:dCGPE_3) are the mass (or normalization) of the wave function

 N(Ψ(⋅,t)):=∥Ψ(⋅,t)∥2:=∫Rd1∑l=−1|ψl(x,t)|2dx≡N(Ψ(⋅,0))=1,t≥0, \hb@xt@.01(1.12)

the magnetization (with )

 M(Ψ(⋅,t)):=∫Rd[|ψ1(x,t)|2−|ψ−1(x,t)|2]dx≡M(Ψ(⋅,0))=M \hb@xt@.01(1.13)

and the energy per particle

 E(Ψ(⋅,t)) = ∫Rd{1∑l=−1(12|∇ψl|2+V(x)|ψl|2)+(βn−βs)|ψ1|2|ψ−1|2 \hb@xt@.01(1.14) +βn2|ψ0|4+βn+βs2[|ψ1|4+|ψ−1|4+2|ψ0|2(|ψ1|2+|ψ−1|2)] +βs(¯ψ−1ψ20¯ψ1+ψ−1¯ψ20ψ1)}dx≡E(Ψ(⋅,0)),t≥0.

A fundamental problem in studying BEC is to find the condensate stationary states , in particular the ground state which is the lowest energy stationary state. In other words, the ground state, , is obtained from the minimization of the energy functional subject to the conservation of total mass and magnetization:

Find such that

 Eg:=E(Φg)=minΦ∈SE(Φ), \hb@xt@.01(1.15)

where the nonconvex set is defined as

 S={Φ=(ϕ1,ϕ0,ϕ−1)T | ∥Φ∥=1, ∫Rd[|ϕ1(x)|2−|ϕ−1(x)|2]=M, E(Φ)<∞}. \hb@xt@.01(1.16)

This is a nonconvex minimization problem. When and and , the existence of a minimizer of the nonconvex minimization problem (LABEL:eq:minimize) follows from the standard theory [25]. For understanding the uniqueness question note that for all with . Thus additional constraints have to be introduced to show the uniqueness.

As derived in [9], by defining the Lagrangian

 \@fontswitchL(Φ,μ,λ):=E(Φ)−μ(∥ϕ1∥2+∥ϕ0∥2+∥ϕ−1∥2−1)−λ(∥ϕ1∥2−∥ϕ−1∥2−M), \hb@xt@.01(1.17)

we get the Euler-Lagrange equations associated to the minimization problem (LABEL:eq:minimize):

 (μ+λ)ϕ1(x) = [−12∇2+V(x)+(βn+βs)(|ϕ1|2+|ϕ0|2)+(βn−βs)|ϕ−1|2]ϕ1 \hb@xt@.01(1.18) +βs¯ϕ−1ϕ20:=H1ϕ1,
 μϕ0(x) = [−12∇2+V(x)+(βn+βs)(|ϕ1|2+|ϕ−1|2)+βn|ϕ0|2]ϕ0 \hb@xt@.01(1.19) +2βsϕ−1¯ϕ0ϕ1:=H0ϕ0,
 (μ−λ)ϕ−1(x) = [−12∇2+V(x)+(βn+βs)(|ϕ−1|2+|ϕ0|2)+(βn−βs)|ϕ1|2]ϕ−1 \hb@xt@.01(1.20) +βsϕ20¯ϕ1:=H−1ϕ−1.

Here and are the Lagrange multipliers (or chemical potentials) of the coupled GPEs (LABEL:eq:dCGPE_1)-(LABEL:eq:dCGPE_3). In addition, (LABEL:GPEs30)-(LABEL:GPEs32) is also a nonlinear eigenvalue problem with two constraints

 ∥Φ∥2:=∫Rd|Φ(x)|2dx=∫Rd1∑l=−1|ϕl(x)|2dx:=1∑l=−1∥ϕl∥2=1, \hb@xt@.01(1.21) ∥ϕ1∥2−∥ϕ−1∥2:=∫Rd[|ϕ1(x)|2−|ϕ−1(x)|2]dx=M. \hb@xt@.01(1.22)

In fact, the nonlinear eigenvalue problem (LABEL:GPEs30)-(LABEL:GPEs32) can also be obtained from the coupled GPEs (LABEL:eq:dCGPE_1)-(LABEL:eq:dCGPE_3) by plugging () with

 μ1=μ+λ,μ0=μ,μ−1=μ−λ⟺μ1+μ−1=2μ0. \hb@xt@.01(1.23)

Thus it is also called as time-independent coupled GPEs. In physics literatures, any eigenfucntion of the nonlinear eigenvalue problem (LABEL:GPEs30)-(LABEL:GPEs32) under constraints (LABEL:con11) and (LABEL:con22), whose energy is larger than the energy of the ground state is called as an excited state of the coupled GPEs (LABEL:eq:dCGPE_1)-(LABEL:eq:dCGPE_3).

A widely used numerical method for computing the ground state of a single component condensate is the imaginary time method followed by an appropriate discretization scheme [15, 5, 3] to evolve the resulted gradient flow equation under normalization of the wavefunction, which is mathematically justified by using the normalized gradient flow [5, 3]. However, it is not obvious that this most popular and powerful normalized gradient flow (or imaginary time method) could be directly extended to compute the ground state of spin-1 BEC. The reason is that we only have two normalization conditions (i.e. the two constraints: conservation of total mass and magnetization) which are insufficient to determine the three projection constants for the three components of the wavefunction used in the normalization step. In physics literatures, the imaginary time method is still applied to compute the ground state of spin-1 BEC through the introduction of a random variable to choose the three projection parameters in the projection step [30, 32]. Of course, this is not a determinate and efficient way to compute the ground state of spin-1 BEC due to the choice of the random variable. Recently, Bao and Wang [9] have proposed a continuous normalized gradient flow (CNGF) for computing the ground state of spin-1 BEC. The CNGF is discretized by Crank-Nicolson finite difference method with a proper and very special way to deal with the nonlinear terms and thus the discretization scheme can be proved to be mass and magnetization conservative and energy diminishing in the discretized level [9]. However, at each time step, a fully nonlinear system must be solved which is a little tedious from computational point of view since the CNGF is an integral-differential equations (see details in (LABEL:cngf1)-(LABEL:FPhi)) which involves implicitly the Lagrange multipliers in the normalized gradient flow evolution [9]. The aim of this paper is to introduce a third normalization condition based on the relation between the chemical potentials of spin-1 BEC, in addition to the two existing normalization conditions given by the conservation of the total mass and magnetization. Thus we can completely determine the three projection constants used in the normalization step for the normalized gradient flow. This allows us to develop the most popular and powerful normalized gradient flow or imaginary time method to compute the ground state of spin-1 BEC.

The paper is organized as follows. In section LABEL:s2, the normalized gradient flow is constructed by introducing the third projection or normalization condition for computing the ground state of spin-1 BEC. In section LABEL:s3, the backward-forward Euler sine-pseudospectral method (BESP) is presented to discretize the normalized gradient flow. In section LABEL:s4, ground states of spin-1 BEC are reported with ferromagnetic/antiferromagnetic interaction and harmonic/optical lattice potential in one/three dimensions, respectively. Finally, some conclusions are drawn in section LABEL:s5.

## 2 The normalized gradient flow

In this section, we will construct the normalized gradient flow for computing the ground state of spin-1 BEC by introducing the third normalization condition.

Various algorithms for computing the minimizer of the nonconvex minimization problem (LABEL:eq:minimize) have been studied in literature. For instance, a continuous normalized gradient flow (CNGF) and its discretization that preserve the total mass and magnetization conservation and energy diminishing properties were presented in [9]. Perhaps one of the more popular and efficient techniques for dealing with the normalization constraints in (LABEL:eq:S) is through the following construction: choose a time step and denote time steps as for  . To adapt an efficient algorithm for the solution of the usual gradient flow to the minimization problem under constraints, it is natural to consider the following splitting (or projection) scheme, which was widely used in the physics literature for computing the ground state of BECs:

 ∂tϕ1(x,t) = [12∇2−V(x)−(βn+βs)(|ϕ1|2+|ϕ0|2)−(βn−βs)|ϕ−1|2]ϕ1 \hb@xt@.01(2.1) −βs¯ϕ−1ϕ20, ∂tϕ0(x,t) = [12∇2−V(x)−(βn+βs)(|ϕ1|2+|ϕ−1|2)−βn|ϕ0|2]ϕ0 \hb@xt@.01(2.2) −2βsϕ−1¯ϕ0ϕ1,x∈Rd,tn−1≤t

followed by a projection step as

 ϕ1(x,tn):=ϕ1(x,t+n)=σn1ϕ1(x,t−n), \hb@xt@.01(2.4) ϕ0(x,tn):=ϕ0(x,t+n)=σn0ϕ0(x,t−n),x∈Rd,n≥1, \hb@xt@.01(2.5) ϕ−1(x,tn):=ϕ−1(x,t+n)=σn−1ϕ−1(x,t−n); \hb@xt@.01(2.6)

where () and () are projection constants and they are chosen such that

 ∥Φ(⋅,tn)∥2=1∑l=−1∥ϕl(⋅,tn)∥2=1,∥ϕ1(⋅,tn)∥2−∥ϕ−1(⋅,tn)∥2=M. \hb@xt@.01(2.7)

In fact, the gradient flow (LABEL:eq:GFDN_1)-(LABEL:eq:GFDN_3) can be viewed as applying the steepest decent method to the energy functional in (LABEL:energy) without constraints, and (LABEL:eq:projection_1c)-(LABEL:eq:projection_3c) project the solution back to the unit sphere in order to satisfy the constraints in (LABEL:eq:S). In addition, (LABEL:eq:GFDN_1)-(LABEL:eq:GFDN_3) can also be obtained from the coupled GPEs (LABEL:eq:dCGPE_1)-(LABEL:eq:dCGPE_3) by the change of variable , that is why the algorithm is usually called as the imaginary time method in the physics literatures [15, 5, 3].

Plugging (LABEL:eq:projection_1c)-(LABEL:eq:projection_3c) into (LABEL:cong1), we obtain

 1∑l=−1(σnl)2∥ϕl(⋅,t−n)∥2=1, \hb@xt@.01(2.8) (σn1)2∥ϕ1(⋅,t−n)∥2−(σn−1)2∥ϕ−1(⋅,t−n)∥2=M. \hb@xt@.01(2.9)

There are three unknowns and only two equations in the above nonlinear system, so the solution is undetermined! In order to determine the projection constants (), we need to find an additional equation. Based on the relation between the chemical potentials in (LABEL:chem8) and the continuous normalized gradient flow proposed in [9] for computing the ground state of spin-1 BEC, see details in Appendix A, we propose to use the following equation as the third normalization condition

 σn1σn−1=(σn0)2. \hb@xt@.01(2.10)

Solving the nonlinear system (LABEL:cong2), (LABEL:cong3) and (LABEL:cong5), see details in Appendix B, we get explicitly the projection constants as

 σn0=√1−M2[∥ϕ0(⋅,t−n)∥2+√4(1−M2)∥ϕ1(⋅,t−n)∥2∥ϕ−1(⋅,t−n)∥2+M2∥ϕ0(⋅,t−n)∥4]1/2, \hb@xt@.01(2.11)
 σn1=√1+M−(σn0)2∥ϕ0(⋅,t−n)∥2√2 ∥ϕ1(⋅,t−n)∥,σn−1=√1−M−(σn0)2∥ϕ0(⋅,t−n)∥2√2 ∥ϕ−1(⋅,t−n)∥. \hb@xt@.01(2.12)

From the numerical point of view, the gradient flow (LABEL:eq:GFDN_1)-(LABEL:eq:GFDN_3) can be solved via traditional techniques, and the normalization of the gradient flow is simply achieved by a projection at the end of each time step.

## 3 Backward-forward Euler sine-pseudospectral method

In this section, we will present the backward-forward Euler sine-pseudospectral method (BESP) to discretize the normalized gradient flow (LABEL:eq:GFDN_1)-(LABEL:eq:GFDN_3), (LABEL:eq:projection_1c)-(LABEL:eq:projection_3c) and (LABEL:eq:constant1)-(LABEL:eq:constant2).

Due to the trapping potential given by (LABEL:poten), the solution decays to zero exponentially fast when . Thus in practical computation, we truncate the problem into a bounded computational domain (chosen as an interval in 1D, a rectangle in 2D, and a box in 3D, with , , , , and sufficiently large) with homogeneous Dirichlet boundary conditions.

For simplicity of notation we introduce the method for the case of one spatial dimension () defined over the interval with homogeneous Dirichlet boundary conditions. Generalization to higher dimension are straightforward for tensor product grids, and the results remain valid without modifications. For , we choose the spatial mesh size with for an even positive integer, and let the grid points be

 xl:=a+jh,j=0,1,⋯,M.

Let be the approximation of , , and be the solution vector with component . In the discretization, we use sine-pseudospectral method for spatial derivatives and backward/forward Euler scheme for linear/nonlinear terms in time discretization. The gradient flow (LABEL:eq:GFDN_1)-(LABEL:eq:GFDN_3) is discretized, for and , as

 ϕ∗1,j−ϕn−11,jΔt=12Dsxxϕ∗1|x=xj−α1ϕ∗1,j+Gn−11,j, \hb@xt@.01(3.1) ϕ∗0,j−ϕn−10,jΔt=12Dsxxϕ∗0|x=xj−α0ϕ∗0,j+Gn−10,j,, \hb@xt@.01(3.2) ϕ∗−1,j−ϕn−1−1,jΔt=12Dsxxϕ∗−1|x=xj−α−1ϕ∗−1,j+Gn−1−1,j; \hb@xt@.01(3.3)

where

 Gn−11,j = [α1−V(xj)−(βn+βs)(|ϕn−11,j|2+|ϕn−10,j|2)−(βn−βs)|ϕn−1−1,j|2]ϕn−11,j \hb@xt@.01(3.4) Gn−10,j = [α0−V(xj)−(βn+βs)(|ϕn−11,j|2+|ϕn−1−1,j|2)−βn|ϕn−10,j|2]ϕn−10,j \hb@xt@.01(3.5) −2βsϕn−1−1,j¯ϕn−10,jϕn−11,j, Gn−1−1,j = [α−1−V(xj)−(βn+βs)(|ϕn−1−1,j|2+|ϕn−10,j|2)−(βn−βs)|ϕn−11,j|2]ϕn−1−1,j \hb@xt@.01(3.6) −βs(ϕn−10,j)2¯ϕn−11,j.

Here, , a pseudospectral differential operator approximation of , is defined as

 DsxxU|x=xj=−M−1∑m=1μ2m(^U)msin(μm(xj−a)),j=1,2,⋯,M−1,

where (), the sine transform coefficients of the vector satisfying , are defined as

 μm=πmb−a,(^U)m=2MM−1∑j=1Ujsin(μm(xj−a)),m=1,2,⋯,M−1;

and () are the stabilization parameters which are chosen in the ‘optimal’ form (such that the time step can be chosen as large as possible) as [4]

 α1=12(bmax1+bmin1),α0=12(bmax0+bmin0),α−1=12(bmax−1+bmin−1); \hb@xt@.01(3.7)

with

 bmax1=max1≤j≤M−1[V(xj)+(βn+βs)(|ϕn−11,j|2+|ϕn−10,j|2)+(βn−βs)|ϕn−1−1,j|2], bmin1=min1≤j≤M−1[V(xj)+(βn+βs)(|ϕn−11,j|2+|ϕn−10,j|2)+(βn−βs)|ϕn−1−1,j|2], bmax0=max1≤j≤M−1[V(xj)+(βn+βs)(|ϕn−11,j|2+|ϕn−1−1,j|2)+βn|ϕn−10,j|2], bmin0=min1≤j≤M−1[V(xj)+(βn+βs)(|ϕn−11,j|2+|ϕn−1−1,j|2)+βn|ϕn−10,j|2],

The homogeneous Dirichlet boundary conditions are discretized as

 ϕ∗1,0=ϕ∗1,M=ϕ∗0,0=ϕ∗0,M=ϕ∗−1,0=ϕ∗−1,M=0. \hb@xt@.01(3.8)

The projection step (LABEL:eq:projection_1c)-(LABEL:eq:projection_1c) is discretized, for and , as

 ϕn1,j=σn1ϕ∗1,j,ϕn0,j=σn0ϕ∗0,j,ϕn−1,j=σn−1ϕ∗−1,j, \hb@xt@.01(3.9)

where

 σn0=√1−M2[∥ϕ∗0∥2+√4(1−M2)∥ϕ∗1∥2∥ϕ∗−1∥2+M2∥ϕ∗0∥4]1/2, \hb@xt@.01(3.10)
 σn1=√1+M−α20∥ϕ∗0∥2√2 ∥ϕ∗1∥,σn−1=√1−M−α20∥ϕ∗0∥2√2 ∥ϕ∗−1∥; \hb@xt@.01(3.11)

with

 ∥ϕ∗1∥2=hM−1∑j=1|ϕ∗1,j|2,∥ϕ∗0∥2=hM−1∑j=1