# An Efficient Compact Finite Difference Method for the Solution of the Gross-Pitaevskii Equation

Authors: Rongpei Zhang; Jia Liu; Guozhong Zhao
Journal: Advances in Condensed Matter Physics
 (2015)

---

## Abstract

We present an efficient, unconditionally stable, and accurate numerical method for the solution of the Gross-Pitaevskii equation. We begin with an introduction on the gradient flow with discrete normalization (GFDN) for computing stationary states of a nonconvex minimization problem. Then we present a new numerical method, CFDM-AIF method, which combines compact finite difference method (CFDM) in space and array-representation integration factor (AIF) method in time. The key features of our methods are as follows: (i) the fourth-order accuracy in space andrth (r
≥
2) accuracy in time which can be achieved and (ii) the significant reduction of storage and CPU cost because of array-representation technique for efficient handling of exponential
matrices. The CFDM-AIF method is implemented to investigate the ground and first excited state solutions of the Gross-Pitaevskii equation in two-dimensional (2D) and three-dimensional (3D) Bose-Einstein condensates (BECs). Numerical results are presented to demonstrate the validity, accuracy, and efficiency of the CFDM-AIF method.


---

## Body

## 1. Introduction

Bose-Einstein condensates (BECs) were first experimentally observed in 1995 [1, 2], while they were theoretically predicted a long time before by S.N. Bose and A. Einstein. Later on, the study of Bose-Einstein condensates (BECs) has attracted considerable interests both theoretically and experimentally. At temperature T much smaller than the critical condensation temperature T
c, the properties of BECs are well described by a macroscopic wave function ψ
(
x
,
t
) which obeys the Gross-Pitaevskii equation (GPE). More precisely, ψ
(
x
,
t
) is solution to the dimensionless time-dependent GPE [3]:(1)
i
∂
∂
t
ψ
x
,
t
=
-
1
2
∇
2
ψ
x
,
t
+
V
x
ψ
x
,
t
+
β
ψ
x
,
t
2
ψ
x
,
t
,for x
∈
R
d, d
=
2,3, t
>
0. Here ∇
2
=
Δ is the Laplace operator. Function V
(
x
) is the external (usually confining) potential. Parameter β is the nonlinearity strength describing the interaction between atoms of the condensate. Equation (1) is also called a mean-field nonlinear Schrödinger equation (NLS) with local cubic nonlinearity.One of the fundamental problems in numerical simulation of BECs is to find its ground state so as to compare the numerical results with experimental observations and to prepare initial data for studying the dynamics of BEC. To find a stationary solution of (1), we write(2)
ψ
x
,
t
=
e
-
i
μ
t
ϕ
x
,where μ is the chemical potential of the condensate and ϕ a real function independent of time. Inserting the above expression into (1) gives the equation(3)
μ
ϕ
x
=
-
1
2
∇
2
ϕ
x
+
V
x
ϕ
x
+
β
ϕ
x
,
t
2
ψ
x
,
t
,
x
∈
R
d
,under the normalization condition(4)
ϕ
2
=
∫
R
d
ϕ
x
2
d
x
=
1
.Any eigenvalue μ can be computed from its corresponding eigenfunction ϕ by(5)
μ
=
μ
ϕ
=
∫
R
d
1
2
∇
ϕ
x
2
+
V
x
ϕ
x
2
+
β
ϕ
x
4
d
x
=
E
ϕ
+
β
2
∫
R
d
ϕ
x
4
d
x
,where the energy E
(
ϕ
) is given by(6)
E
ϕ
=
∫
R
d
1
2
∇
ϕ
2
+
V
x
ϕ
2
+
β
2
ϕ
4
d
x
.From mathematical point of view, the ground stateϕ
g is defined as the minimizer of E
(
ϕ
) over the unit sphere S
=
ϕ
∣
ϕ
=
1
,
E
ϕ
<
∞. It is easy to show that the ground state ϕ
g is the eigenfunction of (3) associated to the lowest eigenvalue of (5). Other eigenfunctions of (3) whose energies are larger than E
(
ϕ
g
) are called as excited states in the physics literatures.In recent years, a great deal of approaches have been proposed for computing the ground state. In fact, these approaches were divided into two types: one method is either based on solving the nonlinear eigenvalue problem [4]; the other approach is the imaginary time method which can be written as a normalized gradient flow formulation [5]. Choose the time step Δ
t
>
0 and set t
n
=
n
Δ
t, n
=
0,1
,
2
,
…. Applying the steepest decent method to the energy functional without constraint and then projecting the solution back to the unit sphere at the end of each time interval. We obtain the following gradient flow with discrete normalization (GFDN):(7)
∂
∂
t
ϕ
x
,
t
=
1
2
∇
2
ϕ
-
V
x
ϕ
-
β
ϕ
2
ϕ
x
,
t
,
t
n
≤
t
<
t
n
+
1
,
(8)
ϕ
x
,
t
n
+
1
=
ϕ
x
,
t
n
+
1
-
ϕ
x
,
t
n
+
1
-
,
ϕ
x
,
0
=
ϕ
0
x
,
x
∈
R
d
,
with
ϕ
0
=
1
,where ϕ
(
x
,
t
n
±
)
=
l
i
m
t
→
t
n
±
ϕ
(
x
,
t
). In fact, the gradient flow (7) can also be obtained from the GPE (1) by changing time t to imaginary time s
=
i
t. That is why the method is also called imaginary time method.The imaginary time method has been used extensively by the physics community and has proved to be powerful. In this work we choose this approach for computing the ground and excited states of BEC. There are many methods discretizing the normalized gradient flow in the imaginary time. These methods include spectral (pseudospectral) methods [3, 5], finite difference method (FDM) [6, 7], and discontinuous Galerkin (DG) method [8–10].Recently, there has been growing interest in high-order compact methods for solving Gross-Pitaevskii equation [11, 12]. It was shown that the high-order difference methods play an important role in the simulation of high frequency wave phenomena. In this paper we apply the compact finite difference method for the space discretization. To build the cFDM, we adopt implicit compact scheme which equals a combination of the nodal derivatives to a combination of the nodal values of the function. By introducing a compact representation of the discretized differential operator, the nodal derivatives are implicitly evaluated by nodal values of the function.The backward Euler and Crank-Nicolson method were usually applied in the time discretization. However large global nonlinear systems need to be solved at each time step. Therefore, number of operation for the nonlinear scheme may be large. Besides that, these time integration methods are limited to second-order accuracy. In this paper, we use integration factor method for time discretization in which an array representation for the linear differential operators is introduced [13]. The array-representation approach is based on the idea of compact Implicit Integration Factor (cIIF) [14]; that is, when discretizing the terms with partial derivatives, regard the unknown solution as a vector with index connected to corresponding variables, while keeping other indexes fixed with unrelated variables. This new approach yields several discretization matrices of a small size that depend only on the number of derivatives in the continuous operators and the number of spatial discretization points in the direction of each derivative.Finally we combine the compact finite difference method in space and array-representation integration factor method in time to get the CFDM-AIF method for solving the GFDN. This method yields fourth-order accuracy in space andrth order accuracy in time. And moreover, the storage and CPU cost can be significantly reduced. The rest of the paper is organized as follows. In Section 2 we present the compact finite difference method and array-representation integration factor method to solve the normalized gradient flow equation in 3D. Ground states and first excited states in 2D and 3D BECs are reported in Section 3. Finally, we summarize our conclusion in Section 4.


## 2. Numerical Methods

In this section, we first review the compact finite difference method. Then we will propose the array-representation integration factor method for fully discretizing the GFDN. We will construct the scheme in 3D and the procedure in 2D will follow the similar procedure.


### 2.1. Compact Finite Difference Method

To build the compact finite difference method for the gradient flow equation (7) in 3D, we rewrite it into the following reaction-diffusion equation formulation:(9)
∂
∂
t
ϕ
x
,
t
=
1
2
∂
2
ϕ
∂
x
2
+
∂
2
ϕ
∂
y
2
+
∂
2
ϕ
∂
z
2
+
F
ϕ
,where F
(
ϕ
)
=
-
V
(
x
)
ϕ
-
β
|
ϕ
|
2
ϕ and the computation domain is Ω
=
[
a
l
,
a
u
]
×
[
b
l
,
b
u
]
×
[
c
l
,
c
u
]. The computation domain Ω is discretized into grids described by the set (
x
j
,
y
k
,
z
l
)
=
(
a
l
+
j
h
x
,
b
l
+
k
h
y
,
c
l
+
k
h
z
), 0
≤
j
≤
N
x, 0
≤
k
≤
N
y, 0
≤
l
≤
N
z, where h
x
=
(
a
u
-
a
l
)
/
N
x, h
y
=
(
b
u
-
b
l
)
/
N
y, h
z
=
(
c
u
-
c
l
)
/
N
z. Let N
x, N
y, N
z be the number of spatial grid points in each spatial direction and let h
x
,
h
y
,
h
z be the grid size, respectively. Let ϕ
j
,
k
,
l represent the approximate solution at the grid point (
j
h
x
,
k
h
y
,
l
h
z
).Settingu
=
∂
2
ϕ
/
∂
x
2, v
=
∂
2
ϕ
/
∂
y
2, and w
=
∂
2
ϕ
/
∂
z
2, we get the discretization of (9) on the grid point as follows:(10)
d
d
t
ϕ
j
,
k
,
l
=
1
2
u
j
,
k
,
l
+
v
j
,
k
,
l
+
w
j
,
k
,
l
+
F
ϕ
j
,
k
,
l
,
1
≤
j
≤
N
x
,
1
≤
k
≤
N
y
,
1
≤
l
≤
N
z
.Define the the following linear mapping in three-dimensional (11)
δ
x
2
ϕ
j
,
k
,
l
=
ϕ
j
-
1
,
k
,
l
-
2
ϕ
j
,
k
,
l
+
ϕ
j
+
1
,
k
,
l
h
x
2
,
L
x
u
j
,
k
,
l
=
1
+
h
x
2
12
δ
x
2
u
j
,
k
,
l
=
u
j
-
1
,
k
,
l
+
10
u
j
,
k
,
l
+
u
j
+
1
,
k
,
l
12
,The linear mapping δ
y
2
ϕ
j
,
k
,
l, δ
z
2
ϕ
j
,
k
,
l and L
y
v
j
,
k
,
l, L
z
w
j
,
k
,
l are similarly defined.By using a Taylor expansion, we get(12)
δ
x
2
ϕ
j
,
k
,
l
=
u
i
,
j
,
k
+
h
x
2
12
δ
x
2
u
j
,
k
,
l
+
O
h
4
=
L
x
u
j
,
k
,
l
+
O
h
4
,where h
=
max
⁡
h
x
,
h
y
,
h
z. Omitting the small terms O
(
h
4
), we obtain the following approximation:(13)
u
j
,
k
,
l
=
L
x
-
1
δ
x
2
ϕ
j
,
k
,
l
,
v
j
,
k
,
l
=
L
y
-
1
δ
y
2
ϕ
j
,
k
,
l
,
w
j
,
k
,
l
=
L
z
-
1
δ
z
2
ϕ
j
,
k
,
l
.Define the three-dimensional array: Φ
=
ϕ
j
,
k
,
l
∣
j
=
1
,
…
,
N
x
,
k
=
1
,
…
,
N
y
,
l
=
1
,
…
,
N
z. The semidiscretization of fourth-order compact finite difference scheme for (9) takes the form(14)
d
d
t
Φ
=
L
x
-
1
δ
x
2
Φ
+
L
y
-
1
δ
y
2
Φ
+
L
z
-
1
δ
z
2
Φ
+
F
Φ
.


### 2.2. Array-Representation Integration Factor Method

In this subsection, we will give the time discretization method for the ODEs (14). To construct the array-representation integration factor method for (14), we multiply the integration factor, e
(
L
x
-
1
δ
x
2
+
L
y
-
1
δ
y
2
+
L
z
-
1
δ
z
2
)
t, to both sides and integrate over one time step from t
n to t
n
+
1. Then we can derive an implicit integration factor method after approximating the integral. For example, the second order IIF takes the form(15)
Φ
~
n
+
1
=
e
L
x
-
1
δ
x
2
+
L
y
-
1
δ
y
2
+
L
z
-
1
δ
z
2
Δ
t
Φ
n
+
Δ
t
2
F
Φ
n
+
Δ
t
2
F
Φ
~
n
+
1
.We can also apply higher order integration factor method in our scheme. See [14] for the scheme with higher order accuracy.In a typical representation of the linear differential operator, the matrixL
x
-
1
δ
x
2
+
L
y
-
1
δ
y
2
+
L
z
-
1
δ
z
2 has a size of N
x
N
y
N
z
×
N
x
N
y
N
z. Although the matrix itself is sparse, its exponential is usually not, leading to prohibitive storage and computing cost for any fine spatial meshes. To avoid computing the exponential of huge matrix, we adopt the array-representation integration factor method which decomposes the matrix into small matrices based on an array representation. The three-dimensional array Φ can be treated as the collection of all such one-dimensional vector on a two-dimensional array (16)
Φ
=
⨂
k
=
1
,
…
,
N
y
l
=
1
,
…
,
N
z
⁡
Φ
:,
k
,
l
,where Φ
:
,
k
,
l is a vector by fixing the last two indices k,  l,  Φ
:
,
k
,
l = (
ϕ
1
,
k
,
l
,
ϕ
2
,
k
,
l
,
…
,
ϕ
N
x
,
k
,
l
)
T.We define the matricesA
m
=
(
D
/
h
m
2
)
A
N
m
×
N
m and B
m
=
B
N
m
×
N
m, m
=
x
,
y
,
z, where(17)
A
N
×
N
=
-
2
1
0
⋯
1
1
-
2
1
⋯
0
0
1
-
2
1
⋯
⋯
⋯
⋯
⋯
⋯
1
0
⋯
1
-
2
N
×
N
,
B
N
×
N
=
5
6
1
12
0
⋯
1
12
1
12
5
6
1
12
⋯
0
0
1
12
5
6
1
12
⋯
⋯
⋯
⋯
⋯
⋯
1
12
0
⋯
1
12
5
6
N
×
N
.With the definition of matrices A
x and B
x, the exponential of linear mapping L
x
-
1
δ
x
2 in the array representation can be written as matrix forms:(18)
e
L
x
-
1
δ
x
2
U
=
⨂
k
=
1
,
…
,
N
y
l
=
1
,
…
,
N
z
⁡
e
B
x
-
1
A
x
U
:,
k
,
l
.The exponential of linear mapping L
y
-
1
δ
y
2 and L
z
-
1
δ
z
2 has similar array representations.Because the linear mappingL
x
-
1
δ
x
2, L
y
-
1
δ
y
2, and L
z
-
1
δ
z
2 satisfy the commuting property, we get(19)
e
L
x
-
1
δ
x
2
+
L
y
-
1
δ
y
2
+
L
z
-
1
δ
z
2
=
e
L
x
-
1
δ
x
2
e
L
y
-
1
δ
y
2
e
L
z
-
1
δ
z
2
.Direct application of (18) and (19) to (15) results in the following array-representation integration factor CFDM with order O
(
h
4
+
Δ
t
2
):
(20)
Φ
~
n
+
1
=
⨂
k
=
1
,
…
,
N
y
l
=
1
,
…
,
N
z
⁡
e
B
x
-
1
A
x
⨂
j
=
1
,
…
,
N
x
l
=
1
,
…
,
N
z
⁡
e
B
y
-
1
A
y
·
⨂
j
=
1
,
…
,
N
x
k
=
1
,
…
,
N
y
⁡
e
B
z
-
1
A
z
U
j
,
k
,
:
j
,
:
,
l
:
,
k
,
l
+
Δ
t
2
F
Φ
~
n
+
1
,where U
=
Φ
n
+
Δ
t
/
2
F
(
Φ
n
). Then we apply normalization constraint to the solution Φ
~
n
+
1:(21)
Φ
n
+
1
=
Φ
~
n
+
1
Φ
~
n
+
1
,where ϕ
2
=
h
x
h
y
h
z
∑
j
=
1
N
x
∑
k
=
1
N
y
∑
l
=
1
N
z
|
ϕ
j
,
k
,
l
|
2.For the 2D case the second order CFDM-AIF scheme is similar as (20) and (21):(22)
Φ
~
n
+
1
=
⨂
k
=
1
,
…
,
N
y
⁡
e
B
x
-
1
A
x
⨂
j
=
1
,
…
,
N
x
⁡
e
B
y
-
1
A
y
U
j
,
:
:,
k
+
Δ
t
2
F
Φ
~
n
+
1
Φ
n
+
1
=
Φ
~
n
+
1
Φ
~
n
+
1
,where U
=
Φ
n
+
Δ
t
/
2
F
(
Φ
n
) and ϕ
2
=
h
x
h
y
∑
j
=
1
N
x
∑
k
=
1
N
y
|
ϕ
j
,
k
|
2.


## 2.1. Compact Finite Difference Method

To build the compact finite difference method for the gradient flow equation (7) in 3D, we rewrite it into the following reaction-diffusion equation formulation:(9)
∂
∂
t
ϕ
x
,
t
=
1
2
∂
2
ϕ
∂
x
2
+
∂
2
ϕ
∂
y
2
+
∂
2
ϕ
∂
z
2
+
F
ϕ
,where F
(
ϕ
)
=
-
V
(
x
)
ϕ
-
β
|
ϕ
|
2
ϕ and the computation domain is Ω
=
[
a
l
,
a
u
]
×
[
b
l
,
b
u
]
×
[
c
l
,
c
u
]. The computation domain Ω is discretized into grids described by the set (
x
j
,
y
k
,
z
l
)
=
(
a
l
+
j
h
x
,
b
l
+
k
h
y
,
c
l
+
k
h
z
), 0
≤
j
≤
N
x, 0
≤
k
≤
N
y, 0
≤
l
≤
N
z, where h
x
=
(
a
u
-
a
l
)
/
N
x, h
y
=
(
b
u
-
b
l
)
/
N
y, h
z
=
(
c
u
-
c
l
)
/
N
z. Let N
x, N
y, N
z be the number of spatial grid points in each spatial direction and let h
x
,
h
y
,
h
z be the grid size, respectively. Let ϕ
j
,
k
,
l represent the approximate solution at the grid point (
j
h
x
,
k
h
y
,
l
h
z
).Settingu
=
∂
2
ϕ
/
∂
x
2, v
=
∂
2
ϕ
/
∂
y
2, and w
=
∂
2
ϕ
/
∂
z
2, we get the discretization of (9) on the grid point as follows:(10)
d
d
t
ϕ
j
,
k
,
l
=
1
2
u
j
,
k
,
l
+
v
j
,
k
,
l
+
w
j
,
k
,
l
+
F
ϕ
j
,
k
,
l
,
1
≤
j
≤
N
x
,
1
≤
k
≤
N
y
,
1
≤
l
≤
N
z
.Define the the following linear mapping in three-dimensional (11)
δ
x
2
ϕ
j
,
k
,
l
=
ϕ
j
-
1
,
k
,
l
-
2
ϕ
j
,
k
,
l
+
ϕ
j
+
1
,
k
,
l
h
x
2
,
L
x
u
j
,
k
,
l
=
1
+
h
x
2
12
δ
x
2
u
j
,
k
,
l
=
u
j
-
1
,
k
,
l
+
10
u
j
,
k
,
l
+
u
j
+
1
,
k
,
l
12
,The linear mapping δ
y
2
ϕ
j
,
k
,
l, δ
z
2
ϕ
j
,
k
,
l and L
y
v
j
,
k
,
l, L
z
w
j
,
k
,
l are similarly defined.By using a Taylor expansion, we get(12)
δ
x
2
ϕ
j
,
k
,
l
=
u
i
,
j
,
k
+
h
x
2
12
δ
x
2
u
j
,
k
,
l
+
O
h
4
=
L
x
u
j
,
k
,
l
+
O
h
4
,where h
=
max
⁡
h
x
,
h
y
,
h
z. Omitting the small terms O
(
h
4
), we obtain the following approximation:(13)
u
j
,
k
,
l
=
L
x
-
1
δ
x
2
ϕ
j
,
k
,
l
,
v
j
,
k
,
l
=
L
y
-
1
δ
y
2
ϕ
j
,
k
,
l
,
w
j
,
k
,
l
=
L
z
-
1
δ
z
2
ϕ
j
,
k
,
l
.Define the three-dimensional array: Φ
=
ϕ
j
,
k
,
l
∣
j
=
1
,
…
,
N
x
,
k
=
1
,
…
,
N
y
,
l
=
1
,
…
,
N
z. The semidiscretization of fourth-order compact finite difference scheme for (9) takes the form(14)
d
d
t
Φ
=
L
x
-
1
δ
x
2
Φ
+
L
y
-
1
δ
y
2
Φ
+
L
z
-
1
δ
z
2
Φ
+
F
Φ
.


## 2.2. Array-Representation Integration Factor Method

In this subsection, we will give the time discretization method for the ODEs (14). To construct the array-representation integration factor method for (14), we multiply the integration factor, e
(
L
x
-
1
δ
x
2
+
L
y
-
1
δ
y
2
+
L
z
-
1
δ
z
2
)
t, to both sides and integrate over one time step from t
n to t
n
+
1. Then we can derive an implicit integration factor method after approximating the integral. For example, the second order IIF takes the form(15)
Φ
~
n
+
1
=
e
L
x
-
1
δ
x
2
+
L
y
-
1
δ
y
2
+
L
z
-
1
δ
z
2
Δ
t
Φ
n
+
Δ
t
2
F
Φ
n
+
Δ
t
2
F
Φ
~
n
+
1
.We can also apply higher order integration factor method in our scheme. See [14] for the scheme with higher order accuracy.In a typical representation of the linear differential operator, the matrixL
x
-
1
δ
x
2
+
L
y
-
1
δ
y
2
+
L
z
-
1
δ
z
2 has a size of N
x
N
y
N
z
×
N
x
N
y
N
z. Although the matrix itself is sparse, its exponential is usually not, leading to prohibitive storage and computing cost for any fine spatial meshes. To avoid computing the exponential of huge matrix, we adopt the array-representation integration factor method which decomposes the matrix into small matrices based on an array representation. The three-dimensional array Φ can be treated as the collection of all such one-dimensional vector on a two-dimensional array (16)
Φ
=
⨂
k
=
1
,
…
,
N
y
l
=
1
,
…
,
N
z
⁡
Φ
:,
k
,
l
,where Φ
:
,
k
,
l is a vector by fixing the last two indices k,  l,  Φ
:
,
k
,
l = (
ϕ
1
,
k
,
l
,
ϕ
2
,
k
,
l
,
…
,
ϕ
N
x
,
k
,
l
)
T.We define the matricesA
m
=
(
D
/
h
m
2
)
A
N
m
×
N
m and B
m
=
B
N
m
×
N
m, m
=
x
,
y
,
z, where(17)
A
N
×
N
=
-
2
1
0
⋯
1
1
-
2
1
⋯
0
0
1
-
2
1
⋯
⋯
⋯
⋯
⋯
⋯
1
0
⋯
1
-
2
N
×
N
,
B
N
×
N
=
5
6
1
12
0
⋯
1
12
1
12
5
6
1
12
⋯
0
0
1
12
5
6
1
12
⋯
⋯
⋯
⋯
⋯
⋯
1
12
0
⋯
1
12
5
6
N
×
N
.With the definition of matrices A
x and B
x, the exponential of linear mapping L
x
-
1
δ
x
2 in the array representation can be written as matrix forms:(18)
e
L
x
-
1
δ
x
2
U
=
⨂
k
=
1
,
…
,
N
y
l
=
1
,
…
,
N
z
⁡
e
B
x
-
1
A
x
U
:,
k
,
l
.The exponential of linear mapping L
y
-
1
δ
y
2 and L
z
-
1
δ
z
2 has similar array representations.Because the linear mappingL
x
-
1
δ
x
2, L
y
-
1
δ
y
2, and L
z
-
1
δ
z
2 satisfy the commuting property, we get(19)
e
L
x
-
1
δ
x
2
+
L
y
-
1
δ
y
2
+
L
z
-
1
δ
z
2
=
e
L
x
-
1
δ
x
2
e
L
y
-
1
δ
y
2
e
L
z
-
1
δ
z
2
.Direct application of (18) and (19) to (15) results in the following array-representation integration factor CFDM with order O
(
h
4
+
Δ
t
2
):
(20)
Φ
~
n
+
1
=
⨂
k
=
1
,
…
,
N
y
l
=
1
,
…
,
N
z
⁡
e
B
x
-
1
A
x
⨂
j
=
1
,
…
,
N
x
l
=
1
,
…
,
N
z
⁡
e
B
y
-
1
A
y
·
⨂
j
=
1
,
…
,
N
x
k
=
1
,
…
,
N
y
⁡
e
B
z
-
1
A
z
U
j
,
k
,
:
j
,
:
,
l
:
,
k
,
l
+
Δ
t
2
F
Φ
~
n
+
1
,where U
=
Φ
n
+
Δ
t
/
2
F
(
Φ
n
). Then we apply normalization constraint to the solution Φ
~
n
+
1:(21)
Φ
n
+
1
=
Φ
~
n
+
1
Φ
~
n
+
1
,where ϕ
2
=
h
x
h
y
h
z
∑
j
=
1
N
x
∑
k
=
1
N
y
∑
l
=
1
N
z
|
ϕ
j
,
k
,
l
|
2.For the 2D case the second order CFDM-AIF scheme is similar as (20) and (21):(22)
Φ
~
n
+
1
=
⨂
k
=
1
,
…
,
N
y
⁡
e
B
x
-
1
A
x
⨂
j
=
1
,
…
,
N
x
⁡
e
B
y
-
1
A
y
U
j
,
:
:,
k
+
Δ
t
2
F
Φ
~
n
+
1
Φ
n
+
1
=
Φ
~
n
+
1
Φ
~
n
+
1
,where U
=
Φ
n
+
Δ
t
/
2
F
(
Φ
n
) and ϕ
2
=
h
x
h
y
∑
j
=
1
N
x
∑
k
=
1
N
y
|
ϕ
j
,
k
|
2.


## 3. Numerical Tests

In this section, we apply the CFDM-AIF scheme to compute the ground and excited states of BEC with harmonic-plus-optical lattice potential in two and three dimensions. In our computation, the stationary state is reached when(23)
max
i
,
j
,
k
⁡
ϕ
i
,
j
,
k
n
+
1
-
ϕ
i
,
j
,
k
n
Δ
t
<
1
0
-
6
.


### 3.1. The Accuracy Test

The CFDM-AIF method (20) has the global accuracy O
(
h
4
+
Δ
t
2
) in space and time. In order to show the accuracy and the rate of convergence for the CFDM-AIF method, we test our scheme for a linear reaction-diffusion equation with exact solution. In this test, we will investigate the accuracy and efficiency of our new scheme by comparison with other methods, such as second-order Runge-Kutta (RK2) method.Example 1.
We consider the following equation onΩ
=
[
0,2
π
]
3:(24)
∂
ϕ
∂
t
=
0.2
Δ
ϕ
+
0.1
ϕwith periodic boundary conditions. The exact solution of this equation is ϕ
=
e
-
0.1
t
c
o
s
(
x
)
+
c
o
s
(
y
)
+
c
o
s
(
z
). The final computation time is t
=
1 at which the L
∞ error is measured. We list the error, order of accuracy, and CPU time for CFDM-AIF and RK2 method in Table 1. The time step is Δ
t
=
1
/
N
x for CFDM-AIF method and Δ
t
=
h
x
2
/
3 for RK2 method. As shown in Table 1, CFDM-AIF scheme achieves the higher order accuracy and requires less CPU time than the RK2 method which requires a much smaller time step and becomes more expensive. In particular for large grid number, N
=
128, the computation time is too long to be acceptable.Table 1
Error, order of accuracy, and CPU time with CFDM-AIF2 and RK2 scheme for Example1.
N
x × N
y × N
z
CFDM-AIF2,Δ
t
=
1
/
N
x
=
h
x
/
2
π
RK2,Δ
t
=
h
x
2
/
3
L
∞ error
Order
CPU (s)
L
∞ error
Order
CPU (s)
1
6 × 16 × 16
5.50 × 10−5
—
0.05
6.95 × 10−3
—
0.09
32 × 32 × 32
3.59 × 10−6
3.94
0.43
1.74 × 10−3
2.00
2.06
64 × 64 × 64
2.65 × 10−7
3.76
5.93
4.36 × 10−4
2.00
61.9
128 × 128 × 128
2.69 × 10−8
3.30
127
1.09 × 10−4
2.00
2448Remark. In the CFDM-AIF scheme (20), the matrix B
x
-
1
A
x is N
x
×
N
x which has order of magnitude smaller than the size of the matrix L
x
-
1
δ
x
2 which is N
x
N
y
N
z
×
N
x
N
y
N
z matrix. This saving in storage is critical for carrying out simulations with larger numbers of spatial grid points. If we use non-AIF method on the computer with 4 GB-RAM, it will run out of memory when N
x
=
N
y
=
N
z
=
32 because this method needs to store matrices with a size of N
x
N
y
N
z
×
N
x
N
y
N
z. In contrast, we can apply the AIF method on the same machine with larger number, even N
x
=
N
y
=
N
z
=
256.


### 3.2. The Ground and Excited States in 2D

We compute the ground and first excited states with different parameters in the harmonic-plus-optical lattice potential by our method. The potential function in (1) is defined as(25)
V
x
,
y
=
1
2
x
2
+
y
2
+
κ
sin
2
⁡
π
x
d
+
sin
2
⁡
π
y
d
,where d is the distance between neighbor wells (lattice constant). The problem is solved on Ω
=
[
-
16,16
]
2 with grid number N
x
=
N
y
=
256 and time step Δ
t
=
1
d
-
3.Example 2.
We first compute the ground state with different parametersβ and d and fix κ
=
100. The initial data for ground state is chosen as Thomas-Fermi approximation: (26)
ϕ
g
x
,
y
=
ϕ
g
T
F
ϕ
g
T
F
,
ϕ
g
T
F
=
μ
g
T
F
-
V
x
,
y
β
,
V
x
,
y
<
μ
g
T
F
0
,
o
t
h
e
r
w
i
s
e
,where μ
g
T
F
=
β
/
π.For fixedd
=
4, Figure 1 shows the ground state and time evolution of energy with β
=
1000, β
=
3000, and β
=
6000. From them, we see that with the increase of β the number of peaks increases. We can also see that the energy diminishes fast and keeps diminishing from the plot of time evolution of the energy. For fixed β
=
100, Figure 2 shows the numerical results with d
=
4, d
=
3, and d
=
2. From it, we see that the decrease in d will cause more peaks in the lattice, but the density at the condensate center decreases very fast.Figure 1
The ground state (top row) and time evolution of energy (bottom row) with fixedd
=
4 and β
=
1000
,
3000
,
6000 from left to right in Example 2.
(a)
(b)Figure 2
The ground state with fixedβ
=
100 and d
=
4,3 and d
=
2 from left to right in Example 2.Then we compute the first excited state by takingϕ
0
(
x
,
y
)
=
2
x
/
π
1
/
2
e
-
(
x
2
+
y
2
)
/
2 in x-direction, ϕ
0
(
x
,
y
)
=
2
y
/
π
1
/
2
e
-
(
x
2
+
y
2
)
/
2 in y-direction, and ϕ
0
(
x
,
y
)
=
2
x
y
/
π
1
/
2
e
-
(
x
2
+
y
2
)
/
2 in both x-direction and y-direction. Figure 3 visualizes the first excited states for β
=
500, κ
=
50, and d
=
4. The interior layers and multiscale structures are observed in the first excited states.Figure 3
The first excited states inx-direction, in y-direction, and in x and y-direction from left to right in Example 2.


### 3.3. The Ground and Excited States in 3D

Example 3.
Finally, we apply CFDM-AIF method to compute the ground and first excited states in 3D under a combined harmonic and optical lattice potential; that is,V
(
x
,
y
,
z
)
=
1
/
2
(
x
2
+
y
2
+
z
2
)
+
50
s
i
n
2
(
π
x
/
4
)
+
s
i
n
2
(
π
y
/
4
)
+
s
i
n
2
(
π
z
/
4
) in (1). The initial data is chosen as ϕ
0
(
x
,
y
,
z
)
=
1
/
π
3
/
4
e
-
(
x
2
+
y
2
+
z
2
)
/
2 for ground state and ϕ
0
(
x
,
y
,
z
)
=
2
x
/
π
3
/
4
e
-
(
x
2
+
y
2
+
z
2
)
/
2 for the first excited state in x-direction ϕ
100, ϕ
0
(
x
,
y
,
z
)
=
2
x
y
/
π
3
/
4
e
-
(
x
2
+
y
2
+
z
2
)
/
2 for the first excited state in x- and y-direction ϕ
110, and ϕ
0
(
x
,
y
)
=
2
3
/
2
x
y
z
/
π
3
/
4
e
-
(
x
2
+
y
2
+
z
2
)
/
2 for the first excited state in x-, y-, and z-direction ϕ
111. The problem is solved on Ω
=
[
-
8,8
]
3 with grid number N
x
=
N
y
=
128 and time step Δ
t
=
1
d
-
3.Figure4 shows the isosurface of ground state |
ϕ
g
|
=
0.01 with β
=
100, β
=
800, and β
=
6400. From it, the number of peaks increases with the increase in β. Figure 5 depicts the isosurface plots of the first excited state for β
=
800 at values of |
ϕ
100
|
=
0.01, |
ϕ
110
|
=
0.01, and |
ϕ
111
|
=
0.01. Similarly, in 2D case, we observe that interior layers appear in the first excited state when β is large. And multiscale structures are observed in first excited states under an optical lattice potential. The numerical results show that the CFDM-AIF method are capable of capturing the multiscale structures accurately and efficiently.Figure 4
The isosurface plots of ground state in Example3 for β
=
100, β
=
800, and β
=
6400.Figure 5
The isosurface plots of first excited state in Example3 with |
ϕ
100
|
=
0.01, |
ϕ
110
|
=
0.01, and |
ϕ
111
|
=
0.01.


## 3.1. The Accuracy Test

The CFDM-AIF method (20) has the global accuracy O
(
h
4
+
Δ
t
2
) in space and time. In order to show the accuracy and the rate of convergence for the CFDM-AIF method, we test our scheme for a linear reaction-diffusion equation with exact solution. In this test, we will investigate the accuracy and efficiency of our new scheme by comparison with other methods, such as second-order Runge-Kutta (RK2) method.Example 1.
We consider the following equation onΩ
=
[
0,2
π
]
3:(24)
∂
ϕ
∂
t
=
0.2
Δ
ϕ
+
0.1
ϕwith periodic boundary conditions. The exact solution of this equation is ϕ
=
e
-
0.1
t
c
o
s
(
x
)
+
c
o
s
(
y
)
+
c
o
s
(
z
). The final computation time is t
=
1 at which the L
∞ error is measured. We list the error, order of accuracy, and CPU time for CFDM-AIF and RK2 method in Table 1. The time step is Δ
t
=
1
/
N
x for CFDM-AIF method and Δ
t
=
h
x
2
/
3 for RK2 method. As shown in Table 1, CFDM-AIF scheme achieves the higher order accuracy and requires less CPU time than the RK2 method which requires a much smaller time step and becomes more expensive. In particular for large grid number, N
=
128, the computation time is too long to be acceptable.Table 1
Error, order of accuracy, and CPU time with CFDM-AIF2 and RK2 scheme for Example1.
N
x × N
y × N
z
CFDM-AIF2,Δ
t
=
1
/
N
x
=
h
x
/
2
π
RK2,Δ
t
=
h
x
2
/
3
L
∞ error
Order
CPU (s)
L
∞ error
Order
CPU (s)
1
6 × 16 × 16
5.50 × 10−5
—
0.05
6.95 × 10−3
—
0.09
32 × 32 × 32
3.59 × 10−6
3.94
0.43
1.74 × 10−3
2.00
2.06
64 × 64 × 64
2.65 × 10−7
3.76
5.93
4.36 × 10−4
2.00
61.9
128 × 128 × 128
2.69 × 10−8
3.30
127
1.09 × 10−4
2.00
2448Remark. In the CFDM-AIF scheme (20), the matrix B
x
-
1
A
x is N
x
×
N
x which has order of magnitude smaller than the size of the matrix L
x
-
1
δ
x
2 which is N
x
N
y
N
z
×
N
x
N
y
N
z matrix. This saving in storage is critical for carrying out simulations with larger numbers of spatial grid points. If we use non-AIF method on the computer with 4 GB-RAM, it will run out of memory when N
x
=
N
y
=
N
z
=
32 because this method needs to store matrices with a size of N
x
N
y
N
z
×
N
x
N
y
N
z. In contrast, we can apply the AIF method on the same machine with larger number, even N
x
=
N
y
=
N
z
=
256.


## 3.2. The Ground and Excited States in 2D

We compute the ground and first excited states with different parameters in the harmonic-plus-optical lattice potential by our method. The potential function in (1) is defined as(25)
V
x
,
y
=
1
2
x
2
+
y
2
+
κ
sin
2
⁡
π
x
d
+
sin
2
⁡
π
y
d
,where d is the distance between neighbor wells (lattice constant). The problem is solved on Ω
=
[
-
16,16
]
2 with grid number N
x
=
N
y
=
256 and time step Δ
t
=
1
d
-
3.Example 2.
We first compute the ground state with different parametersβ and d and fix κ
=
100. The initial data for ground state is chosen as Thomas-Fermi approximation: (26)
ϕ
g
x
,
y
=
ϕ
g
T
F
ϕ
g
T
F
,
ϕ
g
T
F
=
μ
g
T
F
-
V
x
,
y
β
,
V
x
,
y
<
μ
g
T
F
0
,
o
t
h
e
r
w
i
s
e
,where μ
g
T
F
=
β
/
π.For fixedd
=
4, Figure 1 shows the ground state and time evolution of energy with β
=
1000, β
=
3000, and β
=
6000. From them, we see that with the increase of β the number of peaks increases. We can also see that the energy diminishes fast and keeps diminishing from the plot of time evolution of the energy. For fixed β
=
100, Figure 2 shows the numerical results with d
=
4, d
=
3, and d
=
2. From it, we see that the decrease in d will cause more peaks in the lattice, but the density at the condensate center decreases very fast.Figure 1
The ground state (top row) and time evolution of energy (bottom row) with fixedd
=
4 and β
=
1000
,
3000
,
6000 from left to right in Example 2.
(a)
(b)Figure 2
The ground state with fixedβ
=
100 and d
=
4,3 and d
=
2 from left to right in Example 2.Then we compute the first excited state by takingϕ
0
(
x
,
y
)
=
2
x
/
π
1
/
2
e
-
(
x
2
+
y
2
)
/
2 in x-direction, ϕ
0
(
x
,
y
)
=
2
y
/
π
1
/
2
e
-
(
x
2
+
y
2
)
/
2 in y-direction, and ϕ
0
(
x
,
y
)
=
2
x
y
/
π
1
/
2
e
-
(
x
2
+
y
2
)
/
2 in both x-direction and y-direction. Figure 3 visualizes the first excited states for β
=
500, κ
=
50, and d
=
4. The interior layers and multiscale structures are observed in the first excited states.Figure 3
The first excited states inx-direction, in y-direction, and in x and y-direction from left to right in Example 2.


## 3.3. The Ground and Excited States in 3D

Example 3.
Finally, we apply CFDM-AIF method to compute the ground and first excited states in 3D under a combined harmonic and optical lattice potential; that is,V
(
x
,
y
,
z
)
=
1
/
2
(
x
2
+
y
2
+
z
2
)
+
50
s
i
n
2
(
π
x
/
4
)
+
s
i
n
2
(
π
y
/
4
)
+
s
i
n
2
(
π
z
/
4
) in (1). The initial data is chosen as ϕ
0
(
x
,
y
,
z
)
=
1
/
π
3
/
4
e
-
(
x
2
+
y
2
+
z
2
)
/
2 for ground state and ϕ
0
(
x
,
y
,
z
)
=
2
x
/
π
3
/
4
e
-
(
x
2
+
y
2
+
z
2
)
/
2 for the first excited state in x-direction ϕ
100, ϕ
0
(
x
,
y
,
z
)
=
2
x
y
/
π
3
/
4
e
-
(
x
2
+
y
2
+
z
2
)
/
2 for the first excited state in x- and y-direction ϕ
110, and ϕ
0
(
x
,
y
)
=
2
3
/
2
x
y
z
/
π
3
/
4
e
-
(
x
2
+
y
2
+
z
2
)
/
2 for the first excited state in x-, y-, and z-direction ϕ
111. The problem is solved on Ω
=
[
-
8,8
]
3 with grid number N
x
=
N
y
=
128 and time step Δ
t
=
1
d
-
3.Figure4 shows the isosurface of ground state |
ϕ
g
|
=
0.01 with β
=
100, β
=
800, and β
=
6400. From it, the number of peaks increases with the increase in β. Figure 5 depicts the isosurface plots of the first excited state for β
=
800 at values of |
ϕ
100
|
=
0.01, |
ϕ
110
|
=
0.01, and |
ϕ
111
|
=
0.01. Similarly, in 2D case, we observe that interior layers appear in the first excited state when β is large. And multiscale structures are observed in first excited states under an optical lattice potential. The numerical results show that the CFDM-AIF method are capable of capturing the multiscale structures accurately and efficiently.Figure 4
The isosurface plots of ground state in Example3 for β
=
100, β
=
800, and β
=
6400.Figure 5
The isosurface plots of first excited state in Example3 with |
ϕ
100
|
=
0.01, |
ϕ
110
|
=
0.01, and |
ϕ
111
|
=
0.01.


## 4. Conclusion

In this paper we have presented the CFDM-AIF method for computing the ground and first excited states in BECs. The method is based on the compact finite difference method in space and array-representation integration factor method in time. This method can reduce the computational cost significantly in both storage and CPUs and is an efficient and attractive method for 2D and 3D GPE.



---

*Source: 127580-2015-06-09.xml*