# A Time-Splitting and Sine Spectral Method for Dynamics of Dipolar Bose-Einstein Condensate

Authors: Si-Qi Li; Xiang-Gui Li; Dong-Ying Hua
Journal: Advances in Mathematical Physics
 (2013)

---

## Abstract

A two-component Bose-Einstein condensate (BEC) described by two coupled a three-dimension  Gross-Pitaevskii (GP) equations  is considered, where one equation has dipole-dipole interaction while the other one has only the usual s-wave contact interaction, in a cigar trap. The time-splitting and  sine spectral method in space is proposed to discretize the time-dependent equations for computing the dynamics of dipolar BEC. The singularity in the dipole-dipole interaction brings significant difficulties both in mathematical analysis and in numerical simulations. Numerical results are given to show the efficiency of this method.


---

## Body

## 1. Introduction

The achievement of the Bose-Einstein condensation of dilute gases in 1995 marked the beginning of a new era in atomic, molecular, and optical physics. That has attracted much attention both theoretically and experimentally. Most of their properties of these trapped quantum gases are governed by the interactions between particles in the condensate [1]. In the last several years, there has been an investigation for realizing a new kind of quantum gases with the dipolar interaction, acting between particles having a permanent magnetic or electric dipole moment. The experimental realization of a BEC of 52Cr atoms [2, 3] at the University of Stuttgart in 2005 gave new impetus to the theoretical and the numerical investigations on these novel dipolar quantum gases at low temperature. Recently more detailed and controlled experimental results have been obtained, illustrating the effects of phase separation in a multicomponent BEC [4–6]. In these papers, the studies of the binary condensates were limited to the case of s-wave interaction, while a great deal of attention has been drawn recently to the dipolar BEC.In this work, a numerical method for computing the dynamics of the two-component dipolar BEC is considered, where one equation has dipole-dipole interaction and the other has only the usual s-wave contact interaction. However, since the dipole-dipole interaction is of long range, anisotropic, and partially attractive and the computational cost in three dimensions high, the nontrivial task of achieving and controlling the dipolar BEC is thus particularly challenging.This paper is organized as follows. In Section2, a numerical method for computing ground states is presented. In Section 3, numerical results are reported to verify the efficiency of this numerical method. Finally, some concluding remarks are drawn in Section 4.


## 2. Numerical Method for Computing the Dynamics

### 2.1. The Nonlocal Gross-Pitaevskii Equation

The two-component dipolar BEC, confined in a cigar trap, is described by two coupled Gross-Pitaevskii equations. As far as the dipolar interaction is concerned, a convolution term is introduced [7–9] to modify the classical Gross-Pitaevskii equation, which results in the following differential-integral equations (1). Since the transition metal has a magnetic dipole interaction while the alkali metal does not have, we take into account this factor in this system. We take Cr as component 1 and Rb as component 2 [10]. Then the GP equations for this system can be written as
(1)iℏ∂φ1(r→,t)∂t=[-ℏ22m1∇2+V1+U1|φ1|2+U12|φ2|2+Vdip*|φ1|2ℏ22m]φ1,iℏ∂φ2(r→,t)∂t=[-ℏ22m2∇2+V2+U2|φ2|2+U12|φ1|2]φ2,
where φ1, φ2 are the wave functions of components one and two, respectively.  The interatomic and the intercomponent s-wave scattering interactions are described by Uj(j=1,2) and U12, respectively, with the following expressions [11]:
(2)Uj=4πℏ2ajmj,U12=2πℏ2ajm1m2/(m1+m2),j=1,2,
where aj is the scattering length of component j and a12 is that between components 1 and 2. Here ℏ is the Planck constant, mj is the mass of the atom of component j, j=1,2, and Vj is the external trapping potential confining the gas. Generally, that is, Vj(r→)=(mj/2)(ωjx2x2+ωjy2y2+ωjz2z2) with ωjp(p=x,y,z) representing the trap frequency in x, y, and z directions, respectively. The local mean-field Uj|φj|2 represents the s-wave interaction. Vdip is the long-range isotropic dipolar interaction potential between two dipoles and it is defined by
(3)Vdip=μ0μdip24π·1-3(r→·n→)2/|r→|2|r→|3=μ0μdip24π·1-3cos2θ|r→|3,
where θ is the angle between the polarization axis n→ and the relative of two atoms (i.e., cosθ=n→·r→/r, r=|r→|=x2+y2+z2). The wave function φi(x,t) is normalized according to ∥φi∥2=∫R3|φi(r→,t)|2dr→=Ni, i=1,2, where Ni is the number of the atoms in the dipolar BEC.


### 2.2. Time-Splitting and Sine Spectral Numerical Method for Dynamics

The system (1) can be made dimensionless and simplified by adopting a unit system where the units for length, time, and energy are given by a0, 1/ω0, and ℏω0, respectively, with ω0=min{ωjx,ωjy,ωjz}, a0=ℏ/m1ω0 [12].By introducing the dimensionless variablest′=t/ω0,  r→′=r→/a0, φj′=a03/2φj, we obtain the dimensionless GP equations in 3D from (1) as follows:
(4)i∂φ1′∂t′=[-12∇2+V1′(r→)+β11|φ1′|2+β12|φ2′|2+λ(Vdip*|φ1′|2)12]φ1′,i∂φ2′∂t′=[-12∇2+amV2′(r→)+β21|φ1′|2+β22|φ2′|2]φ2′,
where β11=4πa1N1, β12=((1+am)/am)2πa12N2, β21=((1+am)/am)2πa12N1, β22=(4πa2/am)N2, am=m2/m1, λ=m1N1μ0μdip2/3ℏ2a0, and Vdip=(3/4π)·((1-3cos2θ)/r3). In addition the wave functions in (4) satisfy ∫R3|φ1′|2=1, ∫R3|φ2′|2=1. By using the following formula [13]
(5)1r3(1-3(n→·r→)2r2)=-34πδ(r→)-∂nn(1r),
where δ(r→) is the Dirac delta function and ∂nn=n→·∇(n→·∇), we can get
(6)Vdip*|φ1′|2=-|φ1′|2-3∂nn(ξ),
where
(7)ξ(r→,t)=14π∫R31|r→-r→′|·|φ1′(r→′,t)|2dr→′.
And it is easy to see that
(8)∇2ξ=-|φ1′(r→,t)|2.
Plugging (6) into (4) and noticing (7) and (8), we can reformulate GPE (4) into the Schrödinger-Poisson type system
(9)i∂φ1∂t=[-12∇2+V1+(β11-λ)|φ1|2+β12|φ2|2-3λ∂nn(ξ)12]φ1,i∂φ2∂t=[-12∇2+V2+β21|φ1|2+β22|φ2|2]φ2,∇2ξ=-|φ1(r→,t)|2.
In practice, the whole space problem is usually truncated into a bounded computational domain Ω=[a,b]×[c,d]×[e,f] with the homogeneous Dirichlet boundary condition. Let
(10)ΩMKL={(j,k,l)∣j=1,2,…,M-1,k=1,2,…,K-1,l=1,2,…,L-1},ΩMKL0{(j,k,l)∣j=0,1,…,M,k=0,1,…,K,l=0,1,…,L}.
Choose the spatial mesh size as hx=(b-a)/M,  hy=(d-c)/K and  hz=(f-e)/L  and define xj=a+jhx, yk=c+khy, and  zl=e+lhz, j,k,l∈ΩMKL0. Then denote the space YMKL=span{ϕpqs(r→),(p,q,s)∈ΩMKL} with
(11)ϕpqs(r)=sin(up(x-a))sin(uq(y-c))sin(us(z-e)),
where up=pπ/(b-a),  uq=qπ/(d-c), and  us=sπ/(f-e).We propose a time-splitting sine pseudo-spectral method for computing the dynamics of the BEC [12].Fromt=tk to =tk+1, the GP equation (9) is solved by three steps. First, we solve
(12)i∂φ1(r→,t)∂t=-12∇2φ1(r→,t),i∂φ2(r→,t)∂t=-12∇2φ2(r→,t),φ1(r→,t)∣r∈∂Ω=0,φ2(r→,t)∣r∈∂Ω=0
from tk to tk+1/2, followed by solving the nonlinear ODE
(13)i∂φ1∂t=[V1+(β11-λ)|φ1|2+β12|φ2|2-3λ∂nn(ξ)]φ1,i∂φ2∂t=[V2+β21|φ1|2+β22|φ2|2]φ2,∇2ξ=-|φ1(r→,t)|2,φ1(r→,t)∣r∈∂Ω=0,φ2(r→,t)∣r∈∂Ω=0,r→∈Ωtn≤t≤tn+1
for one time step. Again, we solve (12) from tk+1/2 to tk+1.Suppose the exact solutions are(14)φ1(xj,yk,zl,tn)=φ1jkln=∑p=1M-1∑q=1K-1∑s=1L-1φ1pqs(tn)sinjpπMsinkqπKsinlsπL,(15)φ2(xj,yk,zl,tn)=φ2jkln=∑p=1M-1∑q=1K-1∑s=1L-1φ2pqs(tn)sinjpπMsinkqπKsinlsπL.
Substitute (14), (15) into (12); we can find that
(16)iddtφ1pqs(tn)=up2+uq2+us22φ1pqs(tn),iddtφ2pqs(tn)=up2+uq2+us22φ2pqs(tn),(p,q,s)∈ΩMKL0,
which can be solved exactly,  and we obtain
(17)φ1jkl(1)=∑p=1M-1∑q=1K-1∑s=1L-1e-iτ((up2+uq2+us2)/4)×αpqssinjpπMsinkqπKsinlsπL,φ2jkl(1)=∑p=1M-1∑q=1K-1∑s=1L-1e-iτ((up2+uq2+us2)/4)×αpqs′sinjpπMsinkqπKsinlsπL,(j,k,l)∈ΩMKL0,
where
(18)αpqs=8MKL∑pqsφ1jklnsinjpπMsinkqπKsinlsπL,αpqs′=8MKL∑pqsφ2jklnsinjpπMsinkqπKsinlsπL.
Equations (12) will be discretized in space by sine pseudo-spectral method and integrated in time [14]. Next, we will show that (13) can be solved exactly.In fact, fort∈[tn,tn+1], multiplying (13) by the conjugation of φ(r,t), that is, φ(r,t)¯, we get
(19)i∂φ1∂tφ1¯=[V1+(β11-λ)|φ1|2+β12|φ2|2-3λ∂nn(ξ)|φ2|2]φ1φ1¯,i∂φ2∂tφ2¯=[V2+β21|φ1|2+β22|φ2|2]φ2φ2¯,
and we also have
(20)-i∂φ1¯∂tφ1=[V1+(β11-λ)|φ1|2+β12|φ2|2-3λ∂nn(ξ)(β11-λ)|φ1|2]φ1¯φ1,-i∂φ2¯∂tφ2=[V2+β21|φ1|2+β22|φ2|2]φ2¯φ2.
Therefore, subtracting (13) from (14), one obtains
(21)id|φi(r→,t)|2dt=0,i=1,2,
which implies that
(22)|φi(r→,t)|2=|φi(r→,tk)|2,i=1,2∀k,tk≤t≤tk+1.
Substituting (22) into (13), we get a linear ODE
(23)i∂φ1(r→,t)∂t=[V1+(β11-λ)|φ1(r→,tk)|2+β12|φ2(r→,tk)|2-3λ∂nn(ξ)]φ1(r→,t),i∂φ2(r→,t)∂t=[V1+β21|φ1(r→,tk)|2+β22|φ2(r→,tk)|2]φ2(r→,t)
which can be solved exactly. Integrating (19) from tk to t, one gets
(24)φ1jkl(2)=φ1jkl(1)×e-iτ[V1+(β11-λ)|φ1(r→,tn)|2+β12|φ2(r→,tn)|2-3λ∂nn(ξn)],φ2jkl(2)=φ2jkl(1)×e-iτ[V1+(β11-λ)|φ1(r→,tn)|2+β12|φ2(r→,tn)|2-3λ∂nn(ξn)],(j,k,l)∈ΩMKL0,
where ∂nnξn(r)∣jkl=-∑pqs(up2+uq2+us2)γpqssin(jpπ/M)sin(kqπ/K)sin(lsπ/L).The discrete sine transform coefficients of the vectorϕpqs(r)∣(xj,yk,zl)  for (p,q,s)∈ΩMKL are
(25)γpqs=1up2+uq2+us28MKL×∑j=1M-1∑k=1K-1∑l=1L-1|φ1jkl(1)|2sin(jpπM)sin(kqπK)sin(lsπL).
Let φ1jkln and φ2jkln be the approximations of φ1(xj,yk,zl,tn) and φ2(xj,yk,zl,tn), respectively, which are the solution of (9). For n=1,2,…, a second-order time-splitting and sine spectral method for solving (9) via the standard Strang splitting is [14–16]
(26)φ1jkl(1)=∑p=1M-1∑q=1K-1∑s=1L-1e-iτ((up2+uq2+us2)/4)αpqssinjpπM×sinkqπKsinlsπL,φ2jkl(1)=∑p=1M-1∑q=1K-1∑s=1L-1e-iτ((up2+uq2+us2)/4)αpqs′×sinjpπMsinkqπKsinlsπL,φ1jkl(2)=φ1jkl(1)e-iτ[V1+(β11-λ)|φ1(r,tn)|2+β12|φ2(r,tn)|2-3λ∂nn(ξn)],φ2jkl(2)=φ2jkl(1)e-iτ[V1+(β11-λ)|φ1(r,tn)|2+β12|φ2(r,tn)|2-3λ∂nn(ξn)],φ1jkln+1=∑p=1M-1∑q=1K-1∑s=1L-1e-iτ((up2+uq2+us2)/4)βpqs×sinjpπMsinkqπKsinlsπL,φ2jkln+1=∑p=1M-1∑q=1K-1∑s=1L-1e-iτ((up2+uq2+us2)/4)βpqs′×sinjpπMsinkqπKsinlsπL,(j,k,l)∈ΩMKL0,
where
(27)βpqs′=8MKL∑pqsφ2jkl(2)sinjpπMsinkqπKsinlsπL,βpqs=8MKL∑pqsφ1jkl(2)sinjpπMsinkqπKsinlsπL.


## 2.1. The Nonlocal Gross-Pitaevskii Equation

The two-component dipolar BEC, confined in a cigar trap, is described by two coupled Gross-Pitaevskii equations. As far as the dipolar interaction is concerned, a convolution term is introduced [7–9] to modify the classical Gross-Pitaevskii equation, which results in the following differential-integral equations (1). Since the transition metal has a magnetic dipole interaction while the alkali metal does not have, we take into account this factor in this system. We take Cr as component 1 and Rb as component 2 [10]. Then the GP equations for this system can be written as
(1)iℏ∂φ1(r→,t)∂t=[-ℏ22m1∇2+V1+U1|φ1|2+U12|φ2|2+Vdip*|φ1|2ℏ22m]φ1,iℏ∂φ2(r→,t)∂t=[-ℏ22m2∇2+V2+U2|φ2|2+U12|φ1|2]φ2,
where φ1, φ2 are the wave functions of components one and two, respectively.  The interatomic and the intercomponent s-wave scattering interactions are described by Uj(j=1,2) and U12, respectively, with the following expressions [11]:
(2)Uj=4πℏ2ajmj,U12=2πℏ2ajm1m2/(m1+m2),j=1,2,
where aj is the scattering length of component j and a12 is that between components 1 and 2. Here ℏ is the Planck constant, mj is the mass of the atom of component j, j=1,2, and Vj is the external trapping potential confining the gas. Generally, that is, Vj(r→)=(mj/2)(ωjx2x2+ωjy2y2+ωjz2z2) with ωjp(p=x,y,z) representing the trap frequency in x, y, and z directions, respectively. The local mean-field Uj|φj|2 represents the s-wave interaction. Vdip is the long-range isotropic dipolar interaction potential between two dipoles and it is defined by
(3)Vdip=μ0μdip24π·1-3(r→·n→)2/|r→|2|r→|3=μ0μdip24π·1-3cos2θ|r→|3,
where θ is the angle between the polarization axis n→ and the relative of two atoms (i.e., cosθ=n→·r→/r, r=|r→|=x2+y2+z2). The wave function φi(x,t) is normalized according to ∥φi∥2=∫R3|φi(r→,t)|2dr→=Ni, i=1,2, where Ni is the number of the atoms in the dipolar BEC.


## 2.2. Time-Splitting and Sine Spectral Numerical Method for Dynamics

The system (1) can be made dimensionless and simplified by adopting a unit system where the units for length, time, and energy are given by a0, 1/ω0, and ℏω0, respectively, with ω0=min{ωjx,ωjy,ωjz}, a0=ℏ/m1ω0 [12].By introducing the dimensionless variablest′=t/ω0,  r→′=r→/a0, φj′=a03/2φj, we obtain the dimensionless GP equations in 3D from (1) as follows:
(4)i∂φ1′∂t′=[-12∇2+V1′(r→)+β11|φ1′|2+β12|φ2′|2+λ(Vdip*|φ1′|2)12]φ1′,i∂φ2′∂t′=[-12∇2+amV2′(r→)+β21|φ1′|2+β22|φ2′|2]φ2′,
where β11=4πa1N1, β12=((1+am)/am)2πa12N2, β21=((1+am)/am)2πa12N1, β22=(4πa2/am)N2, am=m2/m1, λ=m1N1μ0μdip2/3ℏ2a0, and Vdip=(3/4π)·((1-3cos2θ)/r3). In addition the wave functions in (4) satisfy ∫R3|φ1′|2=1, ∫R3|φ2′|2=1. By using the following formula [13]
(5)1r3(1-3(n→·r→)2r2)=-34πδ(r→)-∂nn(1r),
where δ(r→) is the Dirac delta function and ∂nn=n→·∇(n→·∇), we can get
(6)Vdip*|φ1′|2=-|φ1′|2-3∂nn(ξ),
where
(7)ξ(r→,t)=14π∫R31|r→-r→′|·|φ1′(r→′,t)|2dr→′.
And it is easy to see that
(8)∇2ξ=-|φ1′(r→,t)|2.
Plugging (6) into (4) and noticing (7) and (8), we can reformulate GPE (4) into the Schrödinger-Poisson type system
(9)i∂φ1∂t=[-12∇2+V1+(β11-λ)|φ1|2+β12|φ2|2-3λ∂nn(ξ)12]φ1,i∂φ2∂t=[-12∇2+V2+β21|φ1|2+β22|φ2|2]φ2,∇2ξ=-|φ1(r→,t)|2.
In practice, the whole space problem is usually truncated into a bounded computational domain Ω=[a,b]×[c,d]×[e,f] with the homogeneous Dirichlet boundary condition. Let
(10)ΩMKL={(j,k,l)∣j=1,2,…,M-1,k=1,2,…,K-1,l=1,2,…,L-1},ΩMKL0{(j,k,l)∣j=0,1,…,M,k=0,1,…,K,l=0,1,…,L}.
Choose the spatial mesh size as hx=(b-a)/M,  hy=(d-c)/K and  hz=(f-e)/L  and define xj=a+jhx, yk=c+khy, and  zl=e+lhz, j,k,l∈ΩMKL0. Then denote the space YMKL=span{ϕpqs(r→),(p,q,s)∈ΩMKL} with
(11)ϕpqs(r)=sin(up(x-a))sin(uq(y-c))sin(us(z-e)),
where up=pπ/(b-a),  uq=qπ/(d-c), and  us=sπ/(f-e).We propose a time-splitting sine pseudo-spectral method for computing the dynamics of the BEC [12].Fromt=tk to =tk+1, the GP equation (9) is solved by three steps. First, we solve
(12)i∂φ1(r→,t)∂t=-12∇2φ1(r→,t),i∂φ2(r→,t)∂t=-12∇2φ2(r→,t),φ1(r→,t)∣r∈∂Ω=0,φ2(r→,t)∣r∈∂Ω=0
from tk to tk+1/2, followed by solving the nonlinear ODE
(13)i∂φ1∂t=[V1+(β11-λ)|φ1|2+β12|φ2|2-3λ∂nn(ξ)]φ1,i∂φ2∂t=[V2+β21|φ1|2+β22|φ2|2]φ2,∇2ξ=-|φ1(r→,t)|2,φ1(r→,t)∣r∈∂Ω=0,φ2(r→,t)∣r∈∂Ω=0,r→∈Ωtn≤t≤tn+1
for one time step. Again, we solve (12) from tk+1/2 to tk+1.Suppose the exact solutions are(14)φ1(xj,yk,zl,tn)=φ1jkln=∑p=1M-1∑q=1K-1∑s=1L-1φ1pqs(tn)sinjpπMsinkqπKsinlsπL,(15)φ2(xj,yk,zl,tn)=φ2jkln=∑p=1M-1∑q=1K-1∑s=1L-1φ2pqs(tn)sinjpπMsinkqπKsinlsπL.
Substitute (14), (15) into (12); we can find that
(16)iddtφ1pqs(tn)=up2+uq2+us22φ1pqs(tn),iddtφ2pqs(tn)=up2+uq2+us22φ2pqs(tn),(p,q,s)∈ΩMKL0,
which can be solved exactly,  and we obtain
(17)φ1jkl(1)=∑p=1M-1∑q=1K-1∑s=1L-1e-iτ((up2+uq2+us2)/4)×αpqssinjpπMsinkqπKsinlsπL,φ2jkl(1)=∑p=1M-1∑q=1K-1∑s=1L-1e-iτ((up2+uq2+us2)/4)×αpqs′sinjpπMsinkqπKsinlsπL,(j,k,l)∈ΩMKL0,
where
(18)αpqs=8MKL∑pqsφ1jklnsinjpπMsinkqπKsinlsπL,αpqs′=8MKL∑pqsφ2jklnsinjpπMsinkqπKsinlsπL.
Equations (12) will be discretized in space by sine pseudo-spectral method and integrated in time [14]. Next, we will show that (13) can be solved exactly.In fact, fort∈[tn,tn+1], multiplying (13) by the conjugation of φ(r,t), that is, φ(r,t)¯, we get
(19)i∂φ1∂tφ1¯=[V1+(β11-λ)|φ1|2+β12|φ2|2-3λ∂nn(ξ)|φ2|2]φ1φ1¯,i∂φ2∂tφ2¯=[V2+β21|φ1|2+β22|φ2|2]φ2φ2¯,
and we also have
(20)-i∂φ1¯∂tφ1=[V1+(β11-λ)|φ1|2+β12|φ2|2-3λ∂nn(ξ)(β11-λ)|φ1|2]φ1¯φ1,-i∂φ2¯∂tφ2=[V2+β21|φ1|2+β22|φ2|2]φ2¯φ2.
Therefore, subtracting (13) from (14), one obtains
(21)id|φi(r→,t)|2dt=0,i=1,2,
which implies that
(22)|φi(r→,t)|2=|φi(r→,tk)|2,i=1,2∀k,tk≤t≤tk+1.
Substituting (22) into (13), we get a linear ODE
(23)i∂φ1(r→,t)∂t=[V1+(β11-λ)|φ1(r→,tk)|2+β12|φ2(r→,tk)|2-3λ∂nn(ξ)]φ1(r→,t),i∂φ2(r→,t)∂t=[V1+β21|φ1(r→,tk)|2+β22|φ2(r→,tk)|2]φ2(r→,t)
which can be solved exactly. Integrating (19) from tk to t, one gets
(24)φ1jkl(2)=φ1jkl(1)×e-iτ[V1+(β11-λ)|φ1(r→,tn)|2+β12|φ2(r→,tn)|2-3λ∂nn(ξn)],φ2jkl(2)=φ2jkl(1)×e-iτ[V1+(β11-λ)|φ1(r→,tn)|2+β12|φ2(r→,tn)|2-3λ∂nn(ξn)],(j,k,l)∈ΩMKL0,
where ∂nnξn(r)∣jkl=-∑pqs(up2+uq2+us2)γpqssin(jpπ/M)sin(kqπ/K)sin(lsπ/L).The discrete sine transform coefficients of the vectorϕpqs(r)∣(xj,yk,zl)  for (p,q,s)∈ΩMKL are
(25)γpqs=1up2+uq2+us28MKL×∑j=1M-1∑k=1K-1∑l=1L-1|φ1jkl(1)|2sin(jpπM)sin(kqπK)sin(lsπL).
Let φ1jkln and φ2jkln be the approximations of φ1(xj,yk,zl,tn) and φ2(xj,yk,zl,tn), respectively, which are the solution of (9). For n=1,2,…, a second-order time-splitting and sine spectral method for solving (9) via the standard Strang splitting is [14–16]
(26)φ1jkl(1)=∑p=1M-1∑q=1K-1∑s=1L-1e-iτ((up2+uq2+us2)/4)αpqssinjpπM×sinkqπKsinlsπL,φ2jkl(1)=∑p=1M-1∑q=1K-1∑s=1L-1e-iτ((up2+uq2+us2)/4)αpqs′×sinjpπMsinkqπKsinlsπL,φ1jkl(2)=φ1jkl(1)e-iτ[V1+(β11-λ)|φ1(r,tn)|2+β12|φ2(r,tn)|2-3λ∂nn(ξn)],φ2jkl(2)=φ2jkl(1)e-iτ[V1+(β11-λ)|φ1(r,tn)|2+β12|φ2(r,tn)|2-3λ∂nn(ξn)],φ1jkln+1=∑p=1M-1∑q=1K-1∑s=1L-1e-iτ((up2+uq2+us2)/4)βpqs×sinjpπMsinkqπKsinlsπL,φ2jkln+1=∑p=1M-1∑q=1K-1∑s=1L-1e-iτ((up2+uq2+us2)/4)βpqs′×sinjpπMsinkqπKsinlsπL,(j,k,l)∈ΩMKL0,
where
(27)βpqs′=8MKL∑pqsφ2jkl(2)sinjpπMsinkqπKsinlsπL,βpqs=8MKL∑pqsφ1jkl(2)sinjpπMsinkqπKsinlsπL.


## 3. Numerical Results

### 3.1. Example with the Same Initial Condition

The confining potential is a cigar potential withV(r)=(1/2)(x2+y2+0.04z2). Consider the dynamics of the BEC in the cigar trap. The initial condition is given as follows:
(28)φ1,0=φ1(x,y,z,0)=π-3/4γx1/2γz1/4e-(1/2)(γx(x2+y2)+γzz2),φ2,0=φ2(x,y,z,0)=π-3/4γx1/2γz1/4e-(1/2)(γx(x2+y2)+γzz2),(x,y,z)∈R3.
Here γx=ωx/2ω0, γz=ωz/2ω0, and ω0=min{ωjx,ωjy,ωjz}. We solve this system on [-8,8]3 with hx=hy=hz=1/2 and Δt=0.01.Figure1 shows energy evolutions of dipole BEC. And the energy is conserved. Figures 2 and 3 show the wave function evolutions according to time.Figure 1
The energy evolution according tot.Figure 2
The wave function evolution according to time. Surface plots for|φ1(x,0,z,t)|2 at different times.Figure 3
The wave function evolution according to time. Surface plots for|φ2(x,0,z,t)|2 at different times.


### 3.2. Example with the Different Initial Condition

The confining cigar trap potential isV(r)=(1/2)(x2+y2+0.04z2). Solve the dynamics problem for a dipolar BEC with the cigar trap. The initial condition is
(29)φ1,0=φ1(x,y,z,0)=π-3/4γx1/4γy1/4γz1/4e-(1/2)(γ1xx2+γ1yy2+γ1zz2),φ2,0=φ2(x,y,z,0)=12π-3/8e-(1/2)(γ2xx2+γ2yy2+γ2zz2),(x,y,z)∈R3.
Figure 4 shows energy evolutions of dipole BEC. And the energy is conserved. Figures 5 and 6 show the wave function evolutions according to time.Figure 4
The energy evolution according tot.Figure 5
The wave function evolution according to time. Surface plots for|φ1(x,0,z,t)|2 at different times.Figure 6
The wave function evolution according to time. Surface plots for|φ2(x,0,z,t)|2 at different times.


## 3.1. Example with the Same Initial Condition

The confining potential is a cigar potential withV(r)=(1/2)(x2+y2+0.04z2). Consider the dynamics of the BEC in the cigar trap. The initial condition is given as follows:
(28)φ1,0=φ1(x,y,z,0)=π-3/4γx1/2γz1/4e-(1/2)(γx(x2+y2)+γzz2),φ2,0=φ2(x,y,z,0)=π-3/4γx1/2γz1/4e-(1/2)(γx(x2+y2)+γzz2),(x,y,z)∈R3.
Here γx=ωx/2ω0, γz=ωz/2ω0, and ω0=min{ωjx,ωjy,ωjz}. We solve this system on [-8,8]3 with hx=hy=hz=1/2 and Δt=0.01.Figure1 shows energy evolutions of dipole BEC. And the energy is conserved. Figures 2 and 3 show the wave function evolutions according to time.Figure 1
The energy evolution according tot.Figure 2
The wave function evolution according to time. Surface plots for|φ1(x,0,z,t)|2 at different times.Figure 3
The wave function evolution according to time. Surface plots for|φ2(x,0,z,t)|2 at different times.


## 3.2. Example with the Different Initial Condition

The confining cigar trap potential isV(r)=(1/2)(x2+y2+0.04z2). Solve the dynamics problem for a dipolar BEC with the cigar trap. The initial condition is
(29)φ1,0=φ1(x,y,z,0)=π-3/4γx1/4γy1/4γz1/4e-(1/2)(γ1xx2+γ1yy2+γ1zz2),φ2,0=φ2(x,y,z,0)=12π-3/8e-(1/2)(γ2xx2+γ2yy2+γ2zz2),(x,y,z)∈R3.
Figure 4 shows energy evolutions of dipole BEC. And the energy is conserved. Figures 5 and 6 show the wave function evolutions according to time.Figure 4
The energy evolution according tot.Figure 5
The wave function evolution according to time. Surface plots for|φ1(x,0,z,t)|2 at different times.Figure 6
The wave function evolution according to time. Surface plots for|φ2(x,0,z,t)|2 at different times.


## 4. Conclusion

An efficient numerical method is presented for computing the dynamics of the dipolar Bose-Einstein condensates based on two coupled three-dimensional Gross-Pitaevskii equations where one equation has a dipole-dipole interaction potential and the other one has only the usual s-wave contact interaction. Using equality (5), we can reformulate the GPE for dipolar BEC into a Grosss-Pitaevskii-Poisson type system. Numerical examples are given to show the efficiency of our method. The figures show the evolution of the wave function with time. And in all cases, total energy is conserved. The results agree with the previous work [16]. Numerical results are given to demonstrate the efficiency of our numerical method.



---

*Source: 517395-2013-08-01.xml*