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

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.

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 Section 2, 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 𝑖ℏ𝜕𝜑₁(𝑟,⃗𝑡)=−ℏ𝜕𝑡2𝑚²₁∇²+𝑉₁+𝑈₁𝜑₁²+𝑈₁₂𝜑₂²+𝑉_dip∗𝜑₁²𝜑₁,𝑖ℏ𝜕𝜑₂(𝑟,⃗𝑡)=−ℏ𝜕𝑡2𝑚²₂∇²+𝑉₂+𝑈₂𝜑₂²+𝑈₁₂𝜑₁²𝜑₂, (1) where 𝜑₁, 𝜑₂ are the wave functions of components one and two, respectively. The interatomic and the intercomponent s-wave scattering interactions are described by 𝑈_𝑗(𝑗=1,2) and 𝑈₁₂, respectively, with the following expressions [11]: 𝑈_𝑗=4𝜋ℏ²𝑎_𝑗𝑚_𝑗,𝑈₁₂=2𝜋ℏ²𝑎_𝑗𝑚₁𝑚₂/𝑚₁+𝑚₂,𝑗=1,2, (2) where 𝑎_𝑗 is the scattering length of component 𝑗 and 𝑎₁₂ is that between components 1 and 2. Here ℏ is the Planck constant, 𝑚_𝑗 is the mass of the atom of component 𝑗, 𝑗=1,2, and 𝑉_𝑗 is the external trapping potential confining the gas. Generally, that is, 𝑉_𝑗(𝑟)⃗=(𝑚_𝑗/2)(𝜔²_𝑗𝑥𝑥²+𝜔²_𝑗𝑦𝑦²+𝜔²_𝑗𝑧𝑧²) with 𝜔_𝑗𝑝(𝑝=𝑥,𝑦,𝑧) representing the trap frequency in 𝑥, 𝑦, and 𝑧 directions, respectively. The local mean-field 𝑈_𝑗|𝜑_𝑗|² represents the s-wave interaction. 𝑉_dip is the long-range isotropic dipolar interaction potential between two dipoles and it is defined by 𝑉_dip=𝜇₀𝜇²_dip⋅(⃗𝑛)4𝜋1−3𝑟⋅²/|𝑟|⃗²|𝑟|⃗³=𝜇₀𝜇²_dip⋅4𝜋1−3cos²𝜃|𝑟|⃗³, (3) where 𝜃 is the angle between the polarization axis 𝑛 and the relative of two atoms (i.e., cos𝜃=𝑛⋅𝑟/⃗𝑟, 𝑟=|𝑟|⃗=𝑥²+𝑦²+𝑧²). The wave function 𝜑_𝑖(𝑥,𝑡) is normalized according to ‖𝜑_𝑖‖²=∫_𝑅3|𝜑_𝑖(𝑟,⃗𝑡)|²𝑑𝑟⃗=𝑁_𝑖, 𝑖=1,2, where 𝑁_𝑖 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 𝑎₀, 1/𝜔₀, and ℏ𝜔₀, respectively, with 𝜔₀=min{𝜔_𝑗𝑥,𝜔_𝑗𝑦,𝜔_𝑗𝑧}, 𝑎₀=ℏ/𝑚₁𝜔₀ [12].

By introducing the dimensionless variables 𝑡^′=𝑡/𝜔₀,  𝑟^′=𝑟/⃗𝑎₀, 𝜑^′_𝑗=𝑎₀^3/2𝜑_𝑗, we obtain the dimensionless GP equations in 3D from (1) as follows: 𝑖𝜕𝜑^′₁𝜕𝑡^′=−12∇²+𝑉^′₁(𝑟)⃗+𝛽₁₁𝜑^′₁²+𝛽₁₂𝜑^′₂²+𝜆𝑉_dip∗𝜑^′₁²𝜑^′₁,𝑖𝜕𝜑^′₂𝜕𝑡^′=−12∇²+𝑎_𝑚𝑉^′₂(𝑟)⃗+𝛽₂₁𝜑^′₁²+𝛽₂₂𝜑^′₂²𝜑^′₂, (4) where 𝛽₁₁=4𝜋𝑎₁𝑁₁, 𝛽₁₂=((1+𝑎_𝑚)/𝑎_𝑚)2𝜋𝑎₁₂𝑁₂, 𝛽₂₁=((1+𝑎_𝑚)/𝑎_𝑚)2𝜋𝑎₁₂𝑁₁, 𝛽₂₂=(4𝜋𝑎₂/𝑎_𝑚)𝑁₂, 𝑎_𝑚=𝑚₂/𝑚₁, 𝜆=𝑚₁𝑁₁𝜇₀𝜇²_dip/3ℏ²𝑎₀, and 𝑉_dip=(3/4𝜋)⋅((1−3cos²𝜃)/𝑟³). In addition the wave functions in (4) satisfy ∫_𝑅3|𝜑^′₁|²=1, ∫_𝑅3|𝜑^′₂|²=1. By using the following formula [13] 1𝑟³3(𝑟)⃗1−𝑛⋅²𝑟²3=−𝛿(𝑟)⃗4𝜋−𝜕_𝑛𝑛1𝑟, (5) where 𝛿(𝑟)⃗ is the Dirac delta function and 𝜕_𝑛𝑛=𝑛⋅∇(𝑛⋅∇), we can get 𝑉_dip∗𝜑^′₁²=−𝜑^′₁²−3𝜕_𝑛𝑛(𝜉), (6) where 𝜉(⃗)=1𝑟,𝑡4𝜋_𝑅³1𝑟⃗−𝑟^′⋅𝜑^′₁𝑟^′,𝑡²𝑑𝑟^′. (7) And it is easy to see that ∇²𝜉=−𝜑^′₁(𝑟,⃗𝑡)². (8) Plugging (6) into (4) and noticing (7) and (8), we can reformulate GPE (4) into the Schrödinger-Poisson type system 𝑖𝜕𝜑₁=−1𝜕𝑡2∇²+𝑉₁+𝛽₁₁−𝜆𝜑₁²+𝛽₁₂𝜑₂²−3𝜆𝜕_𝑛𝑛(𝜉)𝜑₁,𝑖𝜕𝜑₂=−1𝜕𝑡2∇²+𝑉₂+𝛽₂₁𝜑₁²+𝛽₂₂𝜑₂²𝜑₂,∇²𝜉=−𝜑₁(𝑟,⃗𝑡)². (9) In practice, the whole space problem is usually truncated into a bounded computational domain Ω=[𝑎,𝑏]×[𝑐,𝑑]×[𝑒,𝑓] with the homogeneous Dirichlet boundary condition. Let Ω_𝑀𝐾𝐿=},Ω𝑗,𝑘,𝑙∣𝑗=1,2,…,𝑀−1,𝑘=1,2,…,𝐾−1,𝑙=1,2,…,𝐿−1_{0𝑀𝐾𝐿}}.𝑗,𝑘,𝑙∣𝑗=0,1,…,𝑀,𝑘=0,1,…,𝐾,𝑙=0,1,…,𝐿 (10) Choose the spatial mesh size as ℎ_𝑥=(𝑏−𝑎)/𝑀,  ℎ_𝑦=(𝑑−𝑐)/𝐾 and ℎ_𝑧=(𝑓−𝑒)/𝐿 and define 𝑥_𝑗=𝑎+𝑗ℎ_𝑥, 𝑦_𝑘=𝑐+𝑘ℎ_𝑦, and 𝑧_𝑙=𝑒+𝑙ℎ_𝑧, 𝑗,𝑘,𝑙∈Ω⁰_𝑀𝐾𝐿. Then denote the space 𝑌_𝑀𝐾𝐿=span{𝜙_𝑝𝑞𝑠(𝑟)⃗,(𝑝,𝑞,𝑠)∈Ω_𝑀𝐾𝐿} with 𝜙_𝑝𝑞𝑠(𝑟)=sin𝑢_𝑝(𝑥−𝑎)sin𝑢_𝑞𝑦−𝑐sin𝑢_𝑠(𝑧−𝑒), (11) where 𝑢_𝑝=𝑝𝜋/(𝑏−𝑎), 𝑢_𝑞=𝑞𝜋/(𝑑−𝑐), and 𝑢_𝑠=𝑠𝜋/(𝑓−𝑒).

We propose a time-splitting sine pseudo-spectral method for computing the dynamics of the BEC [12].

From 𝑡=𝑡_𝑘 to =𝑡_𝑘+1, the GP equation (9) is solved by three steps. First, we solve 𝑖𝜕𝜑₁(𝑟,⃗𝑡)1𝜕𝑡=−2∇²𝜑₁(⃗),𝑖𝑟,𝑡𝜕𝜑₂(𝑟,⃗𝑡)1𝜕𝑡=−2∇²𝜑₂(⃗),𝜑𝑟,𝑡₁(𝑟,⃗𝑡)∣_{𝑟∈𝜕Ω}=0,𝜑₂(𝑟,⃗𝑡)∣_{𝑟∈𝜕Ω}=0 (12) from 𝑡_𝑘 to 𝑡_𝑘+1/2, followed by solving the nonlinear ODE 𝑖𝜕𝜑₁=𝑉𝜕𝑡¹+𝛽₁₁−𝜆𝜑₁²+𝛽₁₂𝜑₂²−3𝜆𝜕_𝑛𝑛(𝜉)𝜑₁,𝑖𝜕𝜑₂=𝑉𝜕𝑡²+𝛽₂₁𝜑₁²+𝛽₂₂𝜑₂²𝜑₂,∇²𝜉=−𝜑₁(𝑟,⃗𝑡)²,𝜑₁(𝑟,⃗𝑡)∣_{𝑟∈𝜕Ω}=0,𝜑₂(𝑟,⃗𝑡)∣_{𝑟∈𝜕Ω}⃗=0,𝑟∈Ω𝑡_𝑛≤𝑡≤𝑡_𝑛+1 (13) for one time step. Again, we solve (12) from 𝑡_𝑘+1/2 to 𝑡_𝑘+1.

Suppose the exact solutions are 𝜑₁𝑥_𝑗,𝑦_𝑘,𝑧_𝑙,𝑡_𝑛=𝜑^𝑛_{1𝑗𝑘𝑙}=^𝑀−1_𝑝=1^𝐾−1_𝑞=1^𝐿−1_𝑠=1𝜑_{1𝑝𝑞𝑠}𝑡_𝑛sin𝑗𝑝𝜋𝑀sin𝑘𝑞𝜋𝐾sin𝑙𝑠𝜋𝐿,𝜑₂𝑥_𝑗,𝑦_𝑘,𝑧_𝑙,𝑡_𝑛=𝜑^𝑛_{2𝑗𝑘𝑙}=^𝑀−1_𝑝=1^𝐾−1_𝑞=1^𝐿−1_𝑠=1𝜑_{2𝑝𝑞𝑠}𝑡_𝑛sin𝑗𝑝𝜋𝑀sin𝑘𝑞𝜋𝐾sin𝑙𝑠𝜋𝐿. (14)(15) Substitute (14), (15) into (12); we can find that 𝑖𝑑𝜑𝑑𝑡^{1𝑝𝑞𝑠}𝑡_𝑛=𝑢²_𝑝+𝑢²_𝑞+𝑢²_𝑠2𝜑_{1𝑝𝑞𝑠}𝑡_𝑛,𝑖𝑑𝜑𝑑𝑡^{2𝑝𝑞𝑠}𝑡_𝑛=𝑢²_𝑝+𝑢²_𝑞+𝑢²_𝑠2𝜑_{2𝑝𝑞𝑠}𝑡_𝑛,𝑝,𝑞,𝑠∈Ω⁰_𝑀𝐾𝐿, (16) which can be solved exactly, and we obtain 𝜑⁽¹⁾_{1𝑗𝑘𝑙}=^𝑀−1_𝑝=1^𝐾−1_𝑞=1^𝐿−1_𝑠=1𝑒^{−𝑖𝜏((𝑢2𝑝+𝑢2𝑞+𝑢2𝑠)/4)}×𝛼_𝑝𝑞𝑠sin𝑗𝑝𝜋𝑀sin𝑘𝑞𝜋𝐾sin𝑙𝑠𝜋𝐿,𝜑⁽¹⁾_{2𝑗𝑘𝑙}=^𝑀−1_𝑝=1^𝐾−1_𝑞=1^𝐿−1_𝑠=1𝑒^{−𝑖𝜏((𝑢2𝑝+𝑢2𝑞+𝑢2𝑠)/4)}×𝛼^′_𝑝𝑞𝑠sin𝑗𝑝𝜋𝑀sin𝑘𝑞𝜋𝐾sin𝑙𝑠𝜋𝐿,𝑗,𝑘,𝑙∈Ω⁰_𝑀𝐾𝐿, (17) where 𝛼_𝑝𝑞𝑠=8𝑀𝐾𝐿_𝑝𝑞𝑠𝜑^𝑛_{1𝑗𝑘𝑙}sin𝑗𝑝𝜋𝑀sin𝑘𝑞𝜋𝐾sin𝑙𝑠𝜋𝐿,𝛼^′_𝑝𝑞𝑠=8𝑀𝐾𝐿_𝑝𝑞𝑠𝜑^𝑛_{2𝑗𝑘𝑙}sin𝑗𝑝𝜋𝑀sin𝑘𝑞𝜋𝐾sin𝑙𝑠𝜋𝐿. (18) 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, for 𝑡∈[𝑡_𝑛,𝑡_𝑛+1], multiplying (13) by the conjugation of 𝜑(𝑟,𝑡), that is, 𝜑(𝑟,𝑡), we get 𝑖𝜕𝜑₁𝜕𝑡𝜑₁=𝑉₁+𝛽₁₁−𝜆𝜑₁²+𝛽₁₂𝜑₂²−3𝜆𝜕_𝑛𝑛(𝜉)𝜑₁𝜑₁,𝑖𝜕𝜑₂𝜕𝑡𝜑₂=𝑉₂+𝛽₂₁𝜑₁²+𝛽₂₂𝜑₂²𝜑₂𝜑₂, (19) and we also have 𝜕−𝑖𝜑₁𝜑𝜕𝑡¹=𝑉₁+𝛽₁₁−𝜆𝜑₁²+𝛽₁₂𝜑₂²−3𝜆𝜕_𝑛𝑛(𝜉)𝜑₁𝜑₁,𝜕−𝑖𝜑₂𝜑𝜕𝑡²=𝑉₂+𝛽₂₁𝜑₁²+𝛽₂₂𝜑₂²𝜑₂𝜑₂. (20) Therefore, subtracting (13) from (14), one obtains 𝑖𝑑𝜑_𝑖(𝑟,⃗𝑡)²𝑑𝑡=0,𝑖=1,2, (21) which implies that 𝜑_𝑖(𝑟,⃗𝑡)²=𝜑_𝑖𝑟,⃗𝑡_𝑘²,𝑖=1,2∀𝑘,𝑡_𝑘≤𝑡≤𝑡_𝑘+1. (22) Substituting (22) into (13), we get a linear ODE 𝑖𝜕𝜑₁(𝑟,⃗𝑡)=𝑉𝜕𝑡¹+𝛽₁₁−𝜆𝜑₁𝑟,⃗𝑡_𝑘²+𝛽₁₂𝜑₂𝑟,⃗𝑡_𝑘²−3𝜆𝜕_𝑛𝑛(𝜉)𝜑₁(⃗),𝑖𝑟,𝑡𝜕𝜑₂(𝑟,⃗𝑡)=𝑉𝜕𝑡¹+𝛽₂₁𝜑₁𝑟,⃗𝑡_𝑘²+𝛽₂₂𝜑₂𝑟,⃗𝑡_𝑘²𝜑₂(𝑟,⃗𝑡) (23) which can be solved exactly. Integrating (19) from 𝑡_𝑘 to 𝑡, one gets 𝜑⁽²⁾_{1𝑗𝑘𝑙}=𝜑⁽¹⁾_{1𝑗𝑘𝑙}×𝑒^{−𝑖𝜏[𝑉1+(𝛽11−𝜆)|𝜑1(⃗𝑟,𝑡𝑛)|2+𝛽12|𝜑2(⃗𝑟,𝑡𝑛)|2−3𝜆𝜕𝑛𝑛(𝜉𝑛)]},𝜑⁽²⁾_{2𝑗𝑘𝑙}=𝜑⁽¹⁾_{2𝑗𝑘𝑙}×𝑒^{−𝑖𝜏[𝑉1+(𝛽11−𝜆)|𝜑1(⃗𝑟,𝑡𝑛)|2+𝛽12|𝜑2(⃗𝑟,𝑡𝑛)|2−3𝜆𝜕𝑛𝑛(𝜉𝑛)]},𝑗,𝑘,𝑙∈Ω⁰_𝑀𝐾𝐿, (24) where 𝜕_𝑛𝑛𝜉^𝑛(𝑟)∣_𝑗𝑘𝑙=−∑_𝑝𝑞𝑠(𝑢²_𝑝+𝑢²_𝑞+𝑢²_𝑠)𝛾_𝑝𝑞𝑠sin(𝑗𝑝𝜋/𝑀)sin(𝑘𝑞𝜋/𝐾)sin(𝑙𝑠𝜋/𝐿).

The discrete sine transform coefficients of the vector 𝜙_𝑝𝑞𝑠(𝑟)∣_{(𝑥𝑗,𝑦𝑘,𝑧𝑙)} for (𝑝,𝑞,𝑠)∈Ω_𝑀𝐾𝐿 are 𝛾_𝑝𝑞𝑠=1𝑢²_𝑝+𝑢²_𝑞+𝑢²_𝑠8×𝑀𝐾𝐿_{𝑀−1𝑗=1𝐾−1𝑘=1𝐿−1𝑙=1}𝜑⁽¹⁾_{1𝑗𝑘𝑙}²sin𝑗𝑝𝜋𝑀sin𝑘𝑞𝜋𝐾sin𝑙𝑠𝜋𝐿. (25) Let 𝜑^𝑛_{1𝑗𝑘𝑙} and 𝜑^𝑛_{2𝑗𝑘𝑙} be the approximations of 𝜑₁(𝑥_𝑗,𝑦_𝑘,𝑧_𝑙,𝑡_𝑛) and 𝜑₂(𝑥_𝑗,𝑦_𝑘,𝑧_𝑙,𝑡_𝑛), respectively, which are the solution of (9). For 𝑛=1,2,…, a second-order time-splitting and sine spectral method for solving (9) via the standard Strang splitting is [14–16] 𝜑⁽¹⁾_{1𝑗𝑘𝑙}=^𝑀−1_𝑝=1^𝐾−1_𝑞=1^𝐿−1_𝑠=1𝑒^{−𝑖𝜏((𝑢2𝑝+𝑢2𝑞+𝑢2𝑠)/4)}𝛼_𝑝𝑞𝑠sin𝑗𝑝𝜋𝑀×sin𝑘𝑞𝜋𝐾sin𝑙𝑠𝜋𝐿,𝜑⁽¹⁾_{2𝑗𝑘𝑙}=^𝑀−1_𝑝=1^𝐾−1_𝑞=1^𝐿−1_𝑠=1𝑒^{−𝑖𝜏((𝑢2𝑝+𝑢2𝑞+𝑢2𝑠)/4)}𝛼^′_𝑝𝑞𝑠×sin𝑗𝑝𝜋𝑀sin𝑘𝑞𝜋𝐾sin𝑙𝑠𝜋𝐿,𝜑⁽²⁾_{1𝑗𝑘𝑙}=𝜑⁽¹⁾_{1𝑗𝑘𝑙}𝑒^{−𝑖𝜏[𝑉1+(𝛽11−𝜆)|𝜑1(𝑟,𝑡𝑛)|2+𝛽12|𝜑2(𝑟,𝑡𝑛)|2−3𝜆𝜕𝑛𝑛(𝜉𝑛)]},𝜑⁽²⁾_{2𝑗𝑘𝑙}=𝜑⁽¹⁾_{2𝑗𝑘𝑙}𝑒^{−𝑖𝜏[𝑉1+(𝛽11−𝜆)|𝜑1(𝑟,𝑡𝑛)|2+𝛽12|𝜑2(𝑟,𝑡𝑛)|2−3𝜆𝜕𝑛𝑛(𝜉𝑛)]},𝜑^𝑛+1_{1𝑗𝑘𝑙}=^𝑀−1_𝑝=1^𝐾−1_𝑞=1^𝐿−1_𝑠=1𝑒^{−𝑖𝜏((𝑢2𝑝+𝑢2𝑞+𝑢2𝑠)/4)}𝛽_𝑝𝑞𝑠×sin𝑗𝑝𝜋𝑀sin𝑘𝑞𝜋𝐾sin𝑙𝑠𝜋𝐿,𝜑^𝑛+1_{2𝑗𝑘𝑙}=^𝑀−1_𝑝=1^𝐾−1_𝑞=1^𝐿−1_𝑠=1𝑒^{−𝑖𝜏((𝑢2𝑝+𝑢2𝑞+𝑢2𝑠)/4)}𝛽^′_𝑝𝑞𝑠×sin𝑗𝑝𝜋𝑀sin𝑘𝑞𝜋𝐾sin𝑙𝑠𝜋𝐿,𝑗,𝑘,𝑙∈Ω⁰_𝑀𝐾𝐿, (26) where 𝛽^′_𝑝𝑞𝑠=8𝑀𝐾𝐿_𝑝𝑞𝑠𝜑⁽²⁾_{2𝑗𝑘𝑙}sin𝑗𝑝𝜋𝑀sin𝑘𝑞𝜋𝐾sin𝑙𝑠𝜋𝐿,𝛽_𝑝𝑞𝑠=8𝑀𝐾𝐿_𝑝𝑞𝑠𝜑⁽²⁾_{1𝑗𝑘𝑙}sin𝑗𝑝𝜋𝑀sin𝑘𝑞𝜋𝐾sin𝑙𝑠𝜋𝐿. (27)

3. Numerical Results

3.1. Example with the Same Initial Condition

The confining potential is a cigar potential with 𝑉(𝑟)=(1/2)(𝑥²+𝑦²+0.04𝑧²). Consider the dynamics of the BEC in the cigar trap. The initial condition is given as follows: 𝜑_1,0=𝜑₁𝑥,𝑦,𝑧,0=𝜋^−3/4𝛾_𝑥^1/2𝛾_𝑧^1/4𝑒^{−(1/2)(𝛾𝑥(𝑥2+𝑦2)+𝛾𝑧𝑧2)},𝜑_2,0=𝜑₂𝑥,𝑦,𝑧,0=𝜋^−3/4𝛾_𝑥^1/2𝛾_𝑧^1/4𝑒^{−(1/2)(𝛾𝑥(𝑥2+𝑦2)+𝛾𝑧𝑧2)},𝑥,𝑦,𝑧∈𝑅³. (28) Here 𝛾_𝑥=𝜔_𝑥/2𝜔₀, 𝛾_𝑧=𝜔_𝑧/2𝜔₀, and 𝜔₀=min{𝜔_𝑗𝑥,𝜔_𝑗𝑦,𝜔_𝑗𝑧}. We solve this system on [−8,8]³ with ℎ_𝑥=ℎ_𝑦=ℎ_𝑧=1/2 and Δ𝑡=0.01.

Figure 1 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 to 𝑡.

Figure 2 

The wave function evolution according to time. Surface plots for |𝜑1(𝑥,0,𝑧,𝑡)|2 at different times.

Figure 3 

The wave function evolution according to time. Surface plots for |𝜑2(𝑥,0,𝑧,𝑡)|2 at different times.

3.2. Example with the Different Initial Condition

The confining cigar trap potential is 𝑉(𝑟)=(1/2)(𝑥²+𝑦²+0.04𝑧²). Solve the dynamics problem for a dipolar BEC with the cigar trap. The initial condition is 𝜑_1,0=𝜑₁𝑥,𝑦,𝑧,0=𝜋^−3/4𝛾_𝑥^1/4𝛾_𝑦^1/4𝛾_𝑧^1/4𝑒^{−(1/2)(𝛾1𝑥𝑥2+𝛾1𝑦𝑦2+𝛾1𝑧𝑧2)},𝜑_2,0=𝜑₂=1𝑥,𝑦,𝑧,02𝜋^−3/8𝑒^{−(1/2)(𝛾2𝑥𝑥2+𝛾2𝑦𝑦2+𝛾2𝑧𝑧2)},𝑥,𝑦,𝑧∈𝑅³. (29) 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 to 𝑡.

Figure 5 

The wave function evolution according to time. Surface plots for |𝜑1(𝑥,0,𝑧,𝑡)|2 at different times.

Figure 6 

The wave function evolution according to time. Surface plots for |𝜑2(𝑥,0,𝑧,𝑡)|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.

Authors’ Contribution

All authors conceived the study, participated in its design and coordination, drafted the paper, participated in the sequence alignment, and read and approved the final paper.

Acknowledgments

This work was supported by National Science Foundation of China (no. 11171032) and Beijing Municipal Education Commission (no. KM201110772017).