Contents lists available atSciVerse ScienceDirect
Journal of Computational Physics
www.elsevier.com/locate/jcp
Efficient numerical methods for computing ground states
of spin-1 Bose–Einstein condensates based
on their characterizations
Weizhu Bao
a, I-Liang Chern
b,
c, Yanzhi Zhang
d,
∗
aDepartment of Mathematics and Center for Computational Science and Engineering, National University of Singapore, Singapore 119076, Singapore
bDepartment of Applied Mathematics and Center of Mathematical Modeling and Scientific Computing, National Chiao Tung University, Hsinchu 30010, Taiwan
cDepartment of Mathematics, National Taiwan University, Taipei 10617, Taiwan
dDepartment of Mathematics and Statistics, Missouri University of Science and Technology, Rolla, MO 65409-0020, USA
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 18 May 2012
Received in revised form 24 April 2013 Accepted 27 June 2013
Available online 10 July 2013
Keywords:
Spin-1 Bose–Einstein condensate Ground state
Ferromagnetic Antiferromagnetic Single-mode approximation
Gradient flow with discrete normalization
In this paper, we propose efficient numerical methods for computing ground states of spin-1 Bose–Einstein condensates (BECs) with/without the Ioffe–Pritchard magnetic field B(x). When B(x)=0, a numerical method is introduced to compute the ground states and it is also applied to study properties of ground states. Numerical results suggest that the densities of mF = ±1 components in ground states are identical for any nonzero B(x).
In particular, if B(x)≡B=0 is a constant, the ground states satisfy the single-mode approximation. When B(x)≡0, efficient and simpler numerical methods are presented to solve the ground states of spin-1 BECs based on their ferromagnetic/antiferromagnetic characterizations. Numerical simulations show that our methods are more efficient than those in the literature. In addition, some conjectures are made from our numerical observations.
©2013 Elsevier Inc. All rights reserved.
1. Introduction
Since its first realization in 1995, Bose–Einstein condensation (BEC) has become an important tool to investigate behav-iors of quantum many-body system. In earlier BEC experiments, atoms were confined in a magnetic trap, where their spin degree of freedom was frozen[1,12]. Recently, the development of optical trapping techniques has enabled to confine atoms independently of their spin orientation and thus result in so-called spinor condensates. The spinor BEC has revealed nu-merous exciting new phenomena which are not possessed by single-component (spin-frozen) condensates. It has provided a unique possibility of exploring fundamental concepts of quantum mechanics in a remarkably controllable and tunable environment[32,33,31].
In the mean-field approximation, a spin-F (F
∈ N
) condensate can be described by coupled Gross–Pitaevskii equations (CGPEs) consisting of 2F+
1 equations, and each of them governs one of the 2F+
1 Zeeman states (mF= −
F,
−
F+
1,
. . . ,
F−
1,
F ). For example, a spin-1 BEC is described by the following dimensionless CGPEs[14,27,8,6,7,32]:*
Corresponding author.E-mail addresses:[email protected](W. Bao),[email protected](I-L. Chern),[email protected](Y. Zhang). URL:http://www.math.nus.edu.sg/~bao/(W. Bao).
0021-9991/$ – see front matter ©2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2013.06.036
i
∂ψ
1(
x,
t)
∂
t=
H+ β
s|ψ
1|
2+ |ψ
0|
2− |ψ
−1|
2ψ1
+ β
sψ
02ψ
−∗1+
Bψ0,
i∂ψ0(
x,
t)
∂
t=
H+ β
s|ψ
1|
2+ |ψ
−1|
2ψ0
+
2β
sψ1ψ
0∗ψ
−1+
B∗ψ1
+
Bψ
−1, i∂ψ
−1(x,
t)
∂
t=
H+ β
s|ψ
−1|
2+ |ψ
0|
2− |ψ
1|
2ψ
−1+ β
sψ
02ψ
1∗+
B∗ψ
0,
(1.1)where x
∈ R
d (for d=
1,
2 or 3) is the Cartesian coordinate vector, t0 is time and
ψ
j
(
x,
t)
is the complex-valued wave function of the j-th ( j=
1,
0,
−
1) component. The operator H is defined byH
= −
1 2∇
2+
V d(
x)
+ β
n 1 j=−1|ψ
j|
2,
(1.2)where Vd
(
x)
represents a d-dimensional external trapping potential and it is determined by the type of system under investigation. For instance, if a three-dimensional (3D) harmonic potential is considered, it takes the form V3(
x)
=
12(
γ
x2x2+
γ
2yy2
+
γ
z2z2)
withγ
x,γ
y andγ
z being the dimensionless trapping frequencies in x-, y- and z-directions, respectively. The constantsβ
n andβ
s describe the spin-independent interaction and spin-exchange interaction, respectively, and they are proportional to the total number of atoms N. Ifβ
n>
0 (resp.<
0), the spin-independent interaction is repulsive (resp. attractive), while whenβ
s>
0 (resp.<
0), the spin-exchange interaction is antiferromagnetic (resp. ferromagnetic). The dimensionless function B(
x)
∈ C
represents the external Ioffe–Pritchard magnetic field[13,16,18]. In addition, f∗represents the complex conjugate of a function f .There are two important invariants of(1.1): the normalization of the wave functions, i.e.,
Ψ (
·,
t)
2=
1 j=−1ψ
j(
·,
t)
2:=
1 j=−1 Rdψ
j(
x,
t)
2dx
≡
Ψ (
·,
0)
2=
1,
t0,
(1.3) withΨ (
x,
t)
= (ψ
1(
x,
t), ψ
0(
x,
t), ψ
−1(
x,
t))
T, and the energyE
Ψ (
·,
t)
:=
E0Ψ (
·,
t)
+
2 ReRd B
ψ
1∗ψ0
+ ψ
0∗ψ
−1 dx≡
EΨ (
·,
0)
,
t0,
(1.4)with Re
(
f)
denoting the real part of a function f andE0
Ψ (
·,
t)
:=
Rd1 j=−1 1 2
|∇ψ
j|
2+
V d(
x)
|ψ
j|
2+
β
n 2|ψ
1|
2+ |ψ
0|
2+ |ψ
−1|
2 2+
β
s 2|ψ
1|
2− |ψ
−1|
2 2+ β
s|ψ
0|
2|ψ
1|
2+ |ψ
−1|
2+
2β
sReψ
1∗ψ
02ψ
−∗1dx
,
(1.5)for t
0, i.e., E0(Ψ (
·,
t))
represents the energy when B(
x)
≡
0. In the case of B(
x)
≡
0, the energy E0(Ψ (
·,
t))
is alsoconserved, i.e., E0
Ψ (
·,
t)
≡
E0Ψ (
·,
0)
,
t0,
when B(
x)
≡
0.
(1.6)Furthermore, the total magnetization is conserved when B
(
x)
≡
0, i.e., MΨ (
·,
t)
:=
Rd
ψ1(
x,
t)
2
−
ψ
−1(x,
t)
2dx
≡
MΨ (
·,
0)
=
M,
t0,
(1.7)with
−
1M1.Ground state, a stationary state with the lowest energy, plays an important role in understanding the properties of BECs. There have been many experimental and mathematical studies on ground states of spin-1 condensates. In [30], the phase diagram of the ground states of spin-1 BECs was first reported in the Thomas–Fermi regimes. The phenomena of broken axisymmetry phase were observed in [26] for spin-1 ferromagnetic condensates. Recently, Matuszewski et al. compared the phase separation of the ground states in the ferromagnetic and antiferromagnetic systems [22,23]. Cao et al. proved the existence of the ground states in one-dimensional condensates[9]. On the other hand, some numerical methods have been proposed in recent literature to compute ground states of spin-1 BECs in the absence of the external Ioffe–Pritchard magnetic field (i.e., B
(
x)
≡
0 in(1.1)). For instance, Bao and Wang proposed in[6]a continuous normalized gradient flow (CNGF) and constructed a Crank–Nicolson finite difference scheme to discretize it. In [7,19], Bao and Lim introduced a gradient flow with discrete normalization (GFDN) and they used the sine pseudo-spectral method to discretize it. It hasbeen pointed out in[7]that GFDN type method is computationally more efficient than CNGF in[6]. Chen et al. proposed a pseudo-arclength continuation method to compute the ground states of spin-1 BECs[10]. To the best of our knowledge, all these methods focus on computing the ground states of spin-1 BECs where B
(
x)
=
0 and so far there are still no numerical reports about the ground states when B(
x)
=
0. In addition, to obtain the ground states in the cases of B(
x)
=
0, all the available methods solve three-component CGPE type equations, which makes the simulations very costly. We notice that when B(
x)
=
0, the ground states of spin-1 BECs can be simplified to a single-mode (resp. two-component) reduction for ferromagnetic (resp. antiferromagnetic) systems[17,34,20]. Thus, the ground states in this case can be effectively found by solving the reduced single or two-component systems instead of the original three-component one.In this paper, we propose (i) a numerical method for computing ground states of spin-1 BECs when B
(
x)
=
0; (ii) effi-cient and simpler methods for the case when B(
x)
≡
0, by taking into account their ferromagnetic and antiferromagnetic characterizations. It is organized as follows. In Section2, we propose a numerical method to compute the ground states when the Ioffe–Pritchard magnetic field B(
x)
=
0. While when B(
x)
≡
0, the ground states of the three-component system (1.1)are characterized by those of the corresponding reduced systems, i.e., the single-mode and two-component reductions for the ferromagnetic and antiferromagnetic condensates, respectively. The reductions of the ground states for ferromagnetic and antiferromagnetic spin-1 BECs are discussed in Section3, followed by their numerical discretizations. Numerical results of ground states as well as comparison between different methods are presented in Section4. In Section5, we draw some conclusions and conjectures based on our numerical observations.2. Numerical methods for ground states with nonzero B
Some numerical methods have been recently proposed in the literature[6,7,19,10]to compute ground states of spin-1 BECs in the absence of the Ioffe–Pritchard magnetic field B
(
x)
. However, there is still no numerical report on the ground states when B(
x)
=
0. In this section, we propose a numerical method for computing ground states of spin-1 BECs with nonzero B(
x)
.When B
(
x)
=
0, the ground stateΦ
g(
x)
= (φ
1,g(
x), φ
0,g(
x), φ
−1,g(
x))
T is defined by minimizing the energy functional E in(1.4)subject to the normalization of wave functions in(1.3), i.e.,Find (
Φ
g∈
S), such thatEg
:=
E(Φ
g)
=
minΦ∈SE
(Φ),
(2.1)where the set S is defined by
S
:=
Φ
= (φ
1, φ
0, φ
−1)
TΦ
2=
1,
E(Φ) <
∞
.
It is easy to see that the ground state
Φ
g defined in(2.1)satisfies the following Euler–Lagrange equationsμφ1(
x)
=
H+ β
s|φ
1|
2+ |φ
0|
2− |φ
−1|
2φ1
+ β
sφ
02φ
−∗1+
Bφ0,
μφ0(
x)
=
H+ β
s|φ
1|
2+ |φ
−1|
2φ0
+
2β
sφ1φ
∗0φ
−1+
B∗φ1
+
Bφ
−1, x∈ R
d,
μφ
−1(x)
=
H+ β
s|φ
−1|
2+ |φ
0|
2− |φ
1|
2φ
−1+ β
sφ
20φ
1∗+
B∗φ0,
(2.2)with the constraint of normalization
Φ(
·)
2=
Rd
φ1(
x)
2
+
φ0(
x)
2
+
φ
−1(x)
2dx
=
1,
(2.3)where the operator H is defined in(1.2). The eigenvalue
μ
is the Lagrange multiplier (or called chemical potential) corre-sponding to the constraint in(2.3), which can be computed from its eigenfunctionΦ
byμ
:=
μ(Φ)
=
E(Φ)
+
Rdβ
n 2|φ
1|
2+ |φ
0|
2+ |φ
−1|
2 2+
β
s 2|φ
1|
2− |φ
−1|
2 2+ β
s|φ
0|
2|φ
1|
2+ |φ
−1|
2+
2β
sReφ
1∗φ
02φ
−∗1dx.
In fact, the Euler–Lagrange equations in(2.2) can also be obtained from the time-dependent GPEs in(1.1)by substituting the ansatz
ψ
j(
x,
t)
=
e−iμtφ
j(
x),
j=
1,
0,
−
1.
(2.4) The eigenfunctionsΦ
of the constrained nonlinear eigenvalue problem (2.2)–(2.3) are usually called stationary states of spin-1 BECs. Among all stationary states, the eigenfunctions with minimum energy is called the ground state and those with larger energies are usually called excited states.Various algorithms have been proposed in the literature to find minimizers of the energy functional under constraints. While the imaginary time method (i.e., replacing t with
−
iτ in(1.1)) is one of the most popular approaches in studying the ground states of BECs. It can be mathematically justified by the normalized gradient flow [11,3]. In this paper, we will develop our numerical methods for computing ground states of spin-1 BECs based on the gradient flow with discrete normalization; see more information in [3]. Choose a time stept
>
0 and define the time sequence as tn=
nt for n
=
0,
1, . . .
. Then in each time interval[
tn,
tn+1]
, the gradient flow with discrete normalization (GFDN) is given by∂φ1(
x,
t)
∂
t= −
H+ β
s|φ
1|
2+ |φ
0|
2− |φ
−1|
2φ
1− β
sφ
20φ
∗−1−
Bφ
0,
∂φ0(
x,
t)
∂
t= −
H+ β
s|φ
1|
2+ |φ
−1|
2φ0
−
2β
sφ1φ
0∗φ
−1−
B∗φ1
−
Bφ
−1,∂φ
−1(x,
t)
∂
t= −
H+ β
s|φ
−1|
2+ |φ
0|
2− |φ
1|
2φ
−1− β
sφ
02φ
1∗−
B∗φ0,
(2.5)followed by a projection step as
φ
j x,
tn++1:=
φ
j(
x,
t − n+1)
Φ(·,
tn−+1)
,
j=
1,
0,
−
1,
(2.6)where
φ
j(
x,
tn±+1)
=
limt→tn±+1φ
j(
x,
t)
( j=
1,
0,
−
1). The gradient flow in (2.5)–(2.6) can be viewed as first applying the steepest descent method to the energy functional in(1.4) without constraint and then projecting the solution back to the unit sphere to satisfy the normalization constraint in(2.3).In order to solve the GFDN numerically, we discretize(2.5)by using the sine pseudo-spectral method for spatial deriva-tives and the backward/forward Euler scheme for linear/nonlinear terms of temporal derivaderiva-tives[4,35]. In the following, we will give a detailed description of our numerical method. Notice that because of the confinement of the external trapping potentials, the wave function
Φ
in(2.5) decays to zero exponentially fast when|
x| → ∞
. Thus, in practical computations we can truncate the problem into a bounded computational domainΩ
with homogeneous Dirichlet boundary conditions, i.e.,φ
j(
x,
t)
|
∂Ω=
0,
t0,
j=
1,
0,
−
1.
(2.7)For simplicity of notations, we will only present the scheme in one-dimensional (1D) cases with a bounded computa-tional domain
Ω
= [
a,
b]
. Generalizations to higher dimensions are straightforward for tensor product grids. For an even integer K>
0, define the spatial mesh sizex
= (
b−
a)/
K>
0 and grid points xk=
a+
kx for 0kK . Let
φ
nj,k be the numerical approximation ofφ
j(
xk,
tn)
andΦ
nj be a vector consisting ofφ
nj,k for the j-th component. DenoteΦ
n a vector with sub-vectorsΦ
nj for j=
1,
0,
−
1. Then over each time interval[
tn,
tn+1]
, we discretize(2.5)asφ
(j1,k)− φ
nj,kt
=
1 2D s xxΦ
(1) jx=xk
−
α
n jφ
(j1,k)− φ
nj,k+
Pnj,k,
1kK−
1,
(2.8)for j
=
1,
0,
−
1, where Pnj,k( j=
1,
0,
−
1) are defined by Pn1,k:= −
V1
(
xk)
+ β
n 1 j=−1φ
nj,k2
+ β
sφ
n1,k2
+
φ
0n,k2
−
φ
−n1,k2
φ
n1,k− β
sφ
n0,k2φ
−n1,k∗−
B(
xk)φ
n0,k,
Pn0,k:= −
V1(xk
)
+ β
n 1 j=−1φ
nj,k2
+ β
sφ
n1,k2
+
φ
−n1,k2
φ
n0,k−
2β
sφ
n1,kφ
0n,k∗φ
−n1,k−
B∗(
xk)φ
n1,k−
B(
xk)φ
n−1,k,
Pn−1,k:= −
V1(xk
)
+ β
n 1 j=−1φ
nj,k2
+ β
sφ
n−1,k2
+
φ
0n,k2
−
φ
n1,k2
φ
−n1,k− β
sφ
0n,k2φ
n1,k∗−
B∗(
xk)φ
0n,k.
Dsxxis a sine pseudo-spectral differential operator approximating
∂
xx, which is defined byDsxxU
x=x k
=
K−1 l=1−
μ
2lUl sinμ
l(
xk−
a)
,
1kK−
1,
(2.9)where
Ul denotes the l-th coefficient of the discrete sine transform of the vector U= (
U1,
U2, . . . ,
UK−1)
T, i.e., Ul=
2 K K−1 k=1 Uksin
μ
l(
xk−
a)
andμ
l=
lπ
b−
a,
1lK−
1.
The constant
α
nj
0 ( j=
1,
0,
−
1) is a stabilization parameter, which is chosen in the “optimal” form (such as the time step can be chosen as large as possible) (see, e.g.,[35,7,19]).The projection step in(2.6)is discretized as
φ
nj,+k1=
φ
(1) j,kΦ
(1)with
Φ
(1)=
x 1 j=−1 K
−1 k=1
φ
(j1,k)2
,
(2.10)for
−
1j1 and 1kK−
1. In addition, the initial condition is discretized asφ
0j,k= φ
j(
xk,
0),
0kK,
j=
1,
0,
−
1,
(2.11) and the boundary conditions areφ
nj,0= φ
nj,K=
0,
n=
0,
1, . . . ,
j=
1,
0,
−
1.
(2.12)The discrete system(2.8),(2.11) and(2.12) can be efficiently solved by the discrete sine transform. In fact, taking the discrete sine transform at both sides of(2.8), we obtain
φ
(j1,l)−
φ
nj,lt
= −
1 2μ
2 lφ
(1) j,l−
α
njφ
(1) j,l−
φ
nj,l+
Pnj,l,
l=
1,
2, . . . ,
K−
1,
j=
1,
0,
−
1,
(2.13)which immediately gives that
φ
(j1,l)=
(
1+
α
n jt
)
φ
nj,l+
tP n j,l 1+ (
α
nj+
μ
2l/
2)
t,
j=
1,
0,
−
1 (2.14) for l=
1,
2, . . . ,
K−
1 and n=
0,
1, . . .
. Since the discrete sine transform is used, the memory required to solve the above system is O(
K)
and computational cost per time step is O(
K ln(
K))
. The simulation is stopped by requiring thatmax
−1j11maxkK−1
|φ
n+1 j,k− φ
nj,k|
t
<
ε,
(2.15)where
ε
is a chosen small tolerance. The resulting solutionΦ
:=
limn→∞Φ
n+1 is the ground state of the spin-1 BECs.3. Numerical methods for ground states with B
=
0In Section 2, we present a numerical method to compute ground states of spin-1 BECs when B
(
x)
=
0, while in this section we will study the ground states when B(
x)
≡
0. In the former case, ground states are minimizers of the energy E in (1.4)under the constraints of normalization in(1.3). However, when B(
x)
≡
0, they also need to satisfy the conservation of magnetization defined in(1.7). In detail, the ground stateΦ
0,g(
x)
when B(
x)
≡
0 is defined byFind (
Φ
0,g∈
S0) such thatE0,g
:=
E0(Φ
0,g)
=
min Φ∈S0E0
(Φ),
(3.1)where E0is the energy functional defined in(1.5)and S0 is a nonconvex set defined as
S0
:=
Φ
= (φ
1, φ0, φ−1)TΦ(
·)
2=
1,
φ1(
·)
2−
φ
−1(·)
2=
M,
E0(Φ) <∞
,
with
−
1M1 a given fixed total magnetization. Note that whenβ
n>
0,|β
s| β
nand lim|x|→∞Vd(
x)
= ∞
, the existence of a minimizer of the nonconvex minimization problem(3.1)follows from the standard theory[9,21]. In fact, E0(
α
· Φ
0,g)
=
E0(Φ
0,g)
for all constant vectorα
= (
eiθ1,
eiθ0,
eiθ−1)
T withθ
1+ θ
−1−
2θ
0= ±
mπ for any integer m.Similarly, the ground state
Φ
0,g(
x)
defined in(3.1)also satisfies the following Euler–Lagrange equations:(μ
+ λ)φ
1(
x)
=
H+ β
s|φ
1|
2+ |φ
0|
2− |φ
−1|
2φ
1+ β
sφ
20φ
−∗1,
μφ0(
x)
=
H+ β
s|φ
1|
2+ |φ
−1|
2φ0
+
2β
sφ1φ
∗0φ
−1, x∈ R
d,
(μ
− λ)φ
−1(x)
=
H+ β
s|φ
−1|
2+ |φ
0|
2− |φ
1|
2φ
−1+ β
sφ
02φ
∗1,
(3.2)with the constraints
Φ(
·)
2:=
1 j=−1 Rdφ
j(
x)
2 dx
=
1,
Rdφ
1(
x)
2
−
φ
−1(
x)
2 dx
=
M,
(3.3)where
μ
andλ
are the Lagrange multipliers (or chemical potentials) of the time-independent CGPEs(3.2)–(3.3). The CGPEs in(3.2)can be also obtained from its time-dependent counterpart(1.1)with B(
x)
≡
0 by substituting the ansatzψ
±1(x,
t)
=
e−i(μ±λ)tφ
±1(x),
ψ0(
x,
t)
=
e−iμtφ0(
x).
(3.4)Recently, there have been some numerical methods proposed in the literature (see, e.g.,[6,7,19,10]) to compute ground states of spin-1 BECs when B
(
x)
=
0 and all of them directly solve the system of three-component equations. However, we notice that when B(
x)
=
0, the ground states of spin-1 BECs in fact can be described by a single-mode or two-component reduction based on their ferromagnetic or antiferromagnetic characterizations[20]. As a result, we can introduce numerical methods based on the characterization of spin-1 BECs so as to reduce the computational costs in computing ground states of spin-1 BECs. In the following, we will discuss the ferromagnetic and antiferromagnetic system separately. For each system, we will start with reviewing its characterization properties and then propose numerical methods to solve the reduced system.3.1. Ferromagnetic system
Experimental observations[17,28,15]and numerical simulations[34,6,7,19]suggest that in ferromagnetic (
β
s<
0) spin-1 BECs, each component of the ground state is a multiple of one single density function. This is so-called the single-modeapproximation (SMA) in the literature, which has been justified rigorously in mathematics by Lin and Chern in [20]. As a result, in this case one can compute just one density function instead of three. To show it, we denote
ρ
(
x)
0, for x∈ R
d, as a scalar real-valued density function and require it satisfy the normalization conditionρ(
·)
2=
Rd
ρ
2(
x)
dx=
1.
(3.5)Let1
Φ
g
(
x)
= (φ
1,g, φ
0,g, φ
−1,g)
T be the ground state of a ferromagnetic spin-1 BECs. Based on the single-mode approxima-tion, we can setφ
j,g(
x)
=
φ
j,g(
x)
=
γ
jρ(
x),
x∈ R
d,
j=
1,
0,
−
1 (3.6) with constantsγ
j0 (for j=
1,
0,
−
1). Noticing that the ground stateΦ
g(
x)
is defined under the constraints of normaliza-tion and magnetizanormaliza-tion described in(3.3), we obtainγ
12+
γ
02+
γ
−21=
1,
γ
12−
γ
−21=
M.
(3.7)Substituting(3.6)into(1.5)and taking(3.7)into account, we obtain E0(Φg
)
=
Rd 1 2∇
ρ
(
x)
2
+
Vd(
x)ρ
2(
x)
+
β
n 2ρ
4(
x)
+
β
s 2 M2+
2γ
02(γ1
+
γ
−1)2ρ
4(
x)
dx.
(3.8)Notice that the minimization of(3.8) over
ρ
and(
γ1
,
γ0
,
γ
−1)
is separable [20]. Thus, we can take minimization of(3.8)first over
(
γ1
,
γ0
,
γ
−1)
and then overρ
. Sinceβ
s<
0, minimizing(3.8)over(
γ1
,
γ0
,
γ
−1)
is equivalent tomax γ1,γ0,γ−1
M2+
2γ
2 0(γ1
+
γ
−1)2,
subject to (3.7),
which gives the constants
γ0
=
1 2 1−
M2,
γ
± 1=
1 2(
1±
M).
(3.9)Taking(3.9)into account, we define the SMA energy Esma(ρ
)
:=
E0(Φg)
=
Rd 1 2∇
ρ
(
x)
2
+
Vd(
x)ρ
2(
x)
+
κ
2ρ
4(
x)
dx,
(3.10)where the constant
κ
= β
n+ β
s. It is obvious that minimizing E0(Φ
g)
in(3.8)with constraints(3.3)is equivalent to mini-mizing Esma(
ρ
)
with the normalization in(3.5) [20].From the above, we know that for a ferromagnetic system, we can in fact solve for the density function
ρ
and then obtain the ground state of the corresponding spin-1 BECs by combining (3.6) and(3.9). In detail, we solve the following single-component minimization problem:1 For simplicity of notation, in the following sections we will also useΦ
g(x)to represent the ground state of spin-1 BECs with B(x)≡0. To distinguish
Find (
ρ
g(
x)
∈
Ssma) such thatEsma,g
:=
Esma(ρg)
=
min ρ∈SsmaEsma(ρ
)
(3.11)over the set Ssma
=
ρ
(
x)
∈ R
Rd
ρ
2(
x)
dx=
1,
Esma(ρ) <∞
.
The Euler–Lagrange equation corresponding to(3.11)is given by
μ
smaρ
(
x)
= −
1 2
∇
2
ρ
(
x)
+
Vd
(
x)ρ
(
x)
+
κρ
3(
x),
(3.12)with the constraint
ρ
(
·)
2=
1. This is a nonlinear eigenvalue problem with the normalization constraint and the eigenvalueμsma
can be computed fromμsma(ρ)
=
Rd 1 2∇
ρ(
x)
2
+
Vd(
x)ρ
2(
x)
+
κρ
4(
x)
dx.
(3.13)Similar to that in Section2, to find the minimizer of (3.11) we solve a gradient flow with discrete normalization and its discretization will be presented for 1D case for simplicity. Generalizations of the method to higher dimensions are straightforward. For t
∈ [
tn,
tn+1]
, the 1D GFDN corresponding to(3.11)is given by∂ρ
(
x,
t)
∂
t=
1 2∇
2ρ
−
V 1(
x)ρ
−
κρ
3,
x∈ [
a,
b],
t∈ [
tn,
tn+1],
(3.14)ρ
x,
tn++1:=
ρ(
x,
t − n+1)
ρ(
·,
tn−+1)
,
x∈ [
a,
b],
(3.15)where the computational domain
[
a,
b]
is chosen to be sufficiently large and homogeneous Dirichlet boundary conditionsρ
(
a,
t)
=
ρ
(
b,
t)
=
0 are imposed. At t=
0, the initial condition is given byρ(
x,
0)
=
ρ
0(
x),
x∈ [
a,
b],
withρ
0(
·) =
1.
(3.16)To discretize (3.14)–(3.16), we use the sine pseudo-spectral method for spatial derivatives and the backward/forward Euler scheme for linear/nonlinear terms of the time derivative. The detailed scheme is given as below:
ρ
k(1)−
ρ
n kt
=
1 2D s xxρ
(1)x=xk
−
α
nρ
(1) k−
ρ
n k+
Fkn,
k=
1,
2, . . . ,
K−
1,
(3.17)ρ
kn+1=
ρ
(1) kρ
(1),
k=
1,
2, . . . ,
K−
1,
n=
0,
1, . . . ,
(3.18) whereρ
nk is the numerical approximation of
ρ
(
xk,
tn)
,ρ
n is a vector with componentsρ
knandρ
(1)=
xKk=−11
(
ρ
k(1))
2. The term Fkn:=
Fρ
n k= −
V1(
xk)ρ
kn+
κ
ρ
n k 3,
k=
1,
2, . . . ,
K−
1,
n=
0,
1, . . . .
(3.19)The operator Dsxxis defined in(2.9). The stabilization parameter
α
n0 is chosen as[35,5,2]α
n=
1 2 max 1kK−1 V1(xk)
+
κ
ρ
kn2+
min 1kK−1 V1(xk)
+
κ
ρ
kn2.
The initial and boundary conditions are discretized as
ρ
k0=
ρ
0(
xk),
k=
0,
1, . . . ,
K;
ρ
0n=
ρ
nK=
0,
n=
0,
1, . . . ,
(3.20) respectively. This discrete system can be efficiently solved in the same manner as that for(2.8),(2.11)and(2.12).3.2. Antiferromagnetic system
In an antiferromagnetic (
β
s>
0) system, the density distribution of atoms in ground states highly depends on the total magnetization M. It was shown in [20] that if M=
0, the ground states have vanishing zeroth mF=
0 component, i.e.,φ
0,g(
x)
≡
0. As a result, the original three-component system is indeed characterized by a two-component reduction. How-ever, if M=
0, the ground states of three components satisfy the SMA as given in (3.6). But different from ferromagnetic systems, in this case the constantsγ
j are not unique and they are given byγ1
=
γ
−1= ξ
andγ0
=
1
−
2ξ
2 for anyconstant
ξ
∈ [
0,
1/
√
2]
. For the detailed mathematical proof, we refer the readers to[20]. In the following, we review the two-component reduction and discretize it for computing ground state of antiferromagnetic spin-1 BECs when B=
0 andM
=
0. For the case of M=
0, the numerical method is identical to that we described in Section3.1and thus it is omitted here for brevity.Let
Φ
g(
x)
be the ground state of an antiferromagnetic spin-1 BEC with B=
0 and M=
0. Sinceφ
0,g(
x)
≡
0 in this case, the ground state energy E0,greduces toE0(Φg
)
=
Rd 1 2|∇φ
1,g|
2+ |∇φ
−1,g|
2+
Vd(
x)
|φ
1,g|
2+ |φ
−1,g|
2+
χ
2|φ
1,g|
4+ |φ
−1,g|
4+
ν
|φ
1,g|
2|φ
−1,g|
2 dx,
(3.21)where the constants
χ
= β
n+ β
s andν
= β
n− β
s. From the constraints of normalization and magnetization in(3.3), it is easy to obtain Rdφ
1,g(
x)
2 dx
=
1+
M 2,
Rdφ
−1,g(
x)
2 dx
=
1−
M 2.
(3.22)On the other hand, we define a two-component energy
Etca(u1,u2)
:=
Rd2 j=1 1 2
|∇
uj|
2+
V d(
x)
|
uj|
2+
χ
2|
u1|
4+ |
u2|
4+
ν
|
u1|
2|
u2|
2dx
.
(3.23)It is easy to verify that the ground state
(φ
1,g, φ
−1,g)
minimizes the energy Etca under the constraints (3.22). Hence, theminimization problem defined in(3.1)to find the ground state of an antiferromagnetic spin-1 condensate can be reduced to the following two-component minimization problem:
Find
((
u1,g,
u2,g)
∈
Stca)
, such thatEtca,g
:=
Etca(
u1,g,
u2,g)
=
min (u1,u2)∈StcaEtca
(
u1,
u2)
(3.24)over the set Stca
=
(
u1,u2)u12
= (
1+
M)/
2,
u22= (
1−
M)/
2,
Etca(u1,u2) <∞
.
The ground state of the associated antiferromagnetic spin-1 BECs can be obtained by
φ1
,g(
x)
=
u1,g(
x),
φ0
,g(
x)
≡
0,
φ
−1,g(
x)
=
u2,g(
x).
(3.25) The Euler–Lagrange equations corresponding to the minimization problem in(3.24)are given byμ
tca 1 u1(x)
=
−
1 2∇
2+
V d(
x)
+
χ
u1(x
)
2
+
ν
u2(x
)
2 u1(x
),
μ
tca2 u2(x)
=
−
1 2∇
2+
V d(
x)
+
ν
u1(x
)
2
+
χ
u2(x
)
2 u2(x
),
x∈ R
d (3.26)with the constraints
u12= (
1+
M)/
2 andu22= (
1−
M)/
2, where the two-component chemical potentials are defined byμ
tcaj=
Rd 1 2|∇
uj|
2+
V d(
x)
|
uj|
2+
χ
|
uj|
4+
ν
|
uj|
2|
u(3−j)|
2 dx,
j=
1,
2.
(3.27)To find the ground states defined in (3.24), a gradient flow with discrete normalization for two-component system is solved, i.e., over each time interval
[
tn,
tn+1]
(for n=
0,
1, . . .
), we solve∂
uj(
x,
t)
∂
t=
1 2∇
2−
V d(
x)
−
χ
|
uj|
2−
ν
|
u(3−j)|
2 uj(
x,
t),
j=
1,
2,
(3.28) u1(x,
tn+1)=
1+
M 2 u1(x,
t−n+1)
u1(·,
tn−+1)
,
u2(x,
tn+1)=
1−
M 2 u2(x,
tn−+1)
u2(·,
t−n+1)
,
(3.29)and at time t
=
0, the initial conditions are given by uj(
x,
0)
=
u0j(
x),
j=
1,
2 withu01(
·) =
1+
M 2,
u 0 2(
·) =
1−
M 2.
(3.30)In practice, the above gradient flow can be solved in a bounded domain
Ω
with homogeneous Dirichlet boundary conditionsu1(x
,
t)
=
u2(x,
t)
=
0,
x∈ ∂Ω,
t0,
(3.31)due to the confinement of the external trapping potentials.
Next we will give the discretization of(3.28)–(3.31)in the 1D case with
Ω
= [
a,
b]
. Let unj,kbe the numerical approxi-mation of uj(
xk,
tn)
for j=
1,
2 and unj be a vector with component unj,k. Then we haveu(j1,k)
−
unj,kt
=
1 2D s xxu( 1) jx=xk
−
α
n j u(j1,k)−
unj,k+
Gnj,k,
j=
1,
2,
(3.32) un1+,k1=
1+
M 2 u(11,k) u(11),
u n+1 2,k=
1−
M 2 u(21,)k u(21),
1kK−
1,
(3.33) whereu(j1)=
xkK=−11
|
u(j,1k)|
2 for j=
1,
2, and Gnj,k= −
V1(xk)
+
χ
unj,k
2
+
ν
un(3−j),k
2unj,k
,
1kK−
1,
j=
1,
2.
The stabilization parameters
α
nj ( j
=
1,
2) are chosen asα
nj=
1 2 max 1kK−1 V1(xk)
+
χ
unj,k
2
+
ν
un(3−j),k
2
+
min 1kK−1 V1(xk)
+
χ
unj,k
2
+
ν
un(3−j),k
2
.
The operator Dsxxis defined in(2.9). The homogeneous Dirichlet boundary conditions are discretized as
un1,0
=
un1,K=
un2,0=
un2,K=
0,
n=
0,
1, . . . ,
(3.34)and the initial conditions are discretized as u01,k
=
u01(
xk),
u02,k=
u0
2
(
xk),
k=
0,
1, . . . ,
K.
(3.35)For each time step, the discrete system(3.32)–(3.35)can be solved in the same manner as that for(2.8),(2.11)and(2.12).
3.3. Relation between different minimization problems
As we have seen, the ground states of spin-1 BECs are always defined by constrained minimization problems, e.g.,(2.1) for B
=
0 and(3.1)for B≡
0. In the following, we will describe the relation between different minimization problems.When B
=
0, the ground state can be obtained by minimizing the energy functional E0 subject to the conservation ofboth normalization and magnetization. While when B
=
0, the energy E is minimized by requiring only the conservation of normalization. To see the relation between these two minimization problems, we introduceΦ
Mg defined byFind (
Φ
Mg∈
SM), such thatEMg
:=
EΦ
Mg=
minΦ∈SME
(Φ),
(3.36)where for a given M
∈ [−
1,
1]
, the set SM is defined asSM
:=
Φ
= (φ
1, φ0, φ−1)TΦ(
·)
2=
1,
φ1(
·)
2−
φ
−1(·)
2=
M,
E(Φ) <
∞
.
That is,
Φ
gM is a minimizer of the energy functional E subject to the conservation of both normalization and magnetization M, where M is a given fixed constant. To obtain the minimizer of the energy E with only the constraint of normalization,one needs to further minimizer EMg in (3.36) with respect to
−
1M1. Thus, the minimization problem in (2.1)is equivalent toE
(Φ
g)
=
min M∈[−1,1]EΦ
gM=
minM∈[−1,1]Φ∈minSME
(Φ),
(3.37)which decomposes(2.1)into two minimization problems. The inner minimization problem is similar to that defined in(3.1) but with different energy functional, while the outer minimization only involves one variable M
∈ [−
1,
1]
. In the form of (3.37), one can see the effect of the magnetization M on the ground states when B=
0.On the other hand, the minimization problem in(3.1)is over three wave functions
φ
j (for j=
1,
0,
−
1) subject to two constraints. These constraints specify the normalization and total magnetization of the ground states. However, to better understand each component in the ground state, we want to know the norm of individual component. To this end, we defineΦ
0α,g(
x)
as the minimizer of the following nonconvex minimization problem:Find (
Φ
α0,g∈
Sα0), such that Eα0,g:=
E0Φ
0α,g=
minΦ∈Sα0E0
(Φ),
(3.38)where
α
∈ [
0,
1]
is a given constant and Sα0 is a nonconvex set defined as S0α:=
Φ
= (φ
1, φ0, φ−1)Tφ0(
·)
2=
α,
φ
±1(·)
2=
1−
α
±
M 2,
E0(Φ) <∞
.
The minimization problem (3.38) has three unknowns
φ
αj ( j=
1,
0,
−
1) and the same number of constraints, which is much easier than the problem in (3.1). In fact, the minimizerΦ
0α,g can be viewed as a ground state which satisfies the conservation of the total normalization and magnetization as described in (3.3)and also the conservation of the norm for each component.It is easy to see that to find the ground states
Φ
0,g in(3.1), one can first fixα
and obtainΦ
0α,g by minimizing(3.38), and then minimize Eα0,goverα
∈ [
0,
1]
. Thus, the problem in(3.1)can be decomposed into two minimizing processes, i.e.,Find (
Φ
0,g∈
S0), such thatE0(Φ0,g
)
=
min α∈[0,1]E0Φ
0α,g=
minα∈[0,1]Φ∈minSα0 E0(Φ). (3.39)
By minimizing Eα0,g over 0
α
1, one can easily obtain the relation between three components in the ground states.4. Numerical results
In this section, we first test the accuracy of spatial discretization in our numerical methods and then apply them to compute ground states of spin-1 condensates. We remark here that the ground states computed numerically, in general, are independent of time step and time discretization. The ground states are studied when B
(
x)
=
0 and some conjectures are made from our numerical observations. In the case of B(
x)
≡
0, we first compare the performance of our methods with those proposed in the literature[6,7,19]. Then we apply our efficient methods to study the ground states of the ferromagnetic and antiferromagnetic spin-1 BECs. In all simulations, the ground states are obtained by setting the toleranceε
=
10−6.4.1. Numerical accuracy
In order to test the accuracy of the spatial discretization of our methods in different cases, we take d
=
1, V1(
x)
=
x2/
2 andβ
n
=
100 in (1.1)and M=
0.
5 in (1.7) for the cases when B=
0. The ground state is computed on a bounded computational domain chosen asΩ
= [−
32,
32]
. LetΦ
g represent the ‘exact’ ground state which is obtained numerically by using a very fine mesh h=
x=
1/
64, and its energy is denoted as E0,g=
E0(Φ
g)
and Eg=
E(Φ
g)
when B(
x)
=
0 and B(
x)
=
0, respectively. Similarly, letΦ
hg be the ground states computed by using the mesh size h and its energy is denoted as Eh0,g=
E0(Φ
hg
)
and Ehg=
E(Φ
hg)
when B(
x)
=
0 and B(
x)
=
0, respectively.Table 1lists the errors in terms of the ground state, i.e.,Φ
g− Φ
hg, and the ground state energy, i.e.,|
E0,g−
Eh0,g|
or|
Eg−
Ehg|
, under different mesh size h for three different cases: (i) ferromagnetic case with B=
0 andβ
s= −
10√
2; (ii) antiferromagnetic case with B
=
0 andβ
s=
10√
2; and (iii) nonzero external magnetic field with B
=
1+
2i andβ
s= −
10√
2.FromTable 1and additional results not shown here for brevity, it is easy to see that our numerical methods are spectral order accurate in space for computing the ground states and its energy. Therefore when the solution changes rapidly with optical lattice potential and/or high accuracy is required, in general, our methods need much less grid points than those low-order finite difference or finite element methods, and thus the memory cost and/or computational cost will be saved significantly, especially in 2D and 3D.
4.2. Numerical results for B
=
0In the following, we study ground states of spin-1 BECs when B
