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

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 nonconvexminimization problem.Thenwe present a new numerical method, CFDM-AIFmethod, 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 and rth (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.


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  much smaller than the critical condensation temperature   , the properties of BECs are well described by a macroscopic wave function (x, ) which obeys the Gross-Pitaevskii equation (GPE).More precisely, (x, ) is solution to the dimensionless time-dependent GPE [3]: ∇ 2  (x, ) +  (x)  (x, ) +       (x, )     2  (x, ) , for x ∈ R  ,  = 2, 3,  > 0.Here ∇ 2 = Δ is the Laplace operator.Function (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 where  is the chemical potential of the condensate and  a real function independent of time.Inserting the above expression into (1) gives the equation where the energy () is given by From mathematical point of view, the ground state   is defined as the minimizer of () over the unit sphere  = { | ‖‖ = 1, () < ∞}.It is easy to show that the ground state   is the eigenfunction of (3) associated to the lowest eigenvalue of (5).Other eigenfunctions of (3) whose energies are larger than (  ) 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 Δ > 0 and set   = Δ,  = 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): where (x,  ±  ) = lim  →  ±  (x, ).In fact, the gradient flow (7) can also be obtained from the GPE (1) by changing time  to imaginary time  = .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][9][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 and th 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.

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.

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: where () = −(x) − || Define the the following linear mapping in three-dimensional The linear mapping  2   ,, ,  2   ,, and L  V ,, , L   ,, are similarly defined.
By using a Taylor expansion, we get where ℎ = max{ℎ  , ℎ y , ℎ  }.Omitting the small terms O(ℎ 4 ), we obtain the following approximation: Define the three-dimensional array: Φ = ( ,, |  = 1, . . .,   ,  = 1, . . .,   ,  = 1, . . .,   ).The semidiscretization of fourth-order compact finite difference scheme for (9) takes the form 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, ) , to both sides and integrate over one time step from   to  +1 .Then we can derive an implicit integration factor method after approximating the integral.For example, the second order IIF takes the form 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 matrix ) has a size of       ×       .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 twodimensional array where Φ(:, , ) is a vector by fixing the last two indices , , Φ(:, , ) = ( 1,, ,  2,, , . . .,    ,, )  .We define the matrices A  = (/ℎ 2  )A   ×  and B  = B   ×  ,  = , , , where ) × , With the definition of matrices A  and B  , the exponential of linear mapping L −1   2  in the array representation can be written as matrix forms: The exponential of linear mapping L −1   2  and L −1   2  has similar array representations.

Numerical Tests
In this section, we apply the CFDM-AIF scheme to compute the ground and excited states of BEC with harmonic-plusoptical lattice potential in two and three dimensions.In our computation, the stationary state is reached when max with periodic boundary conditions.The exact solution of this equation is  =  −0.1 (cos() + cos() + cos()).The final computation time is  = 1 at which the  ∞ 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 Δ = 1/  for CFDM-AIF method and Δ = ℎ 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,  = 128, the computation time is too long to be acceptable.
Remark.In the CFDM-AIF scheme (20), the matrix B −1  A  is   ×   which has order of magnitude smaller than the size of the matrix L −1   2  which is       ×       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   =   =   = 32 because this method needs to store matrices with a size of       ×       .In contrast, we can apply the AIF method on the same machine with larger number, even   =   =   = 256.

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 where  TF  = √/.
For fixed  = 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  = 4,  = 3, and  = 2. From it, we see that the decrease in  will cause more peaks in the lattice, but the density at the condensate center decreases very fast.
Figure 4 shows the isosurface of ground state |  | = 0.01 with  = 100,  = 800, and  = 6400.From it, the number of peaks increases with the increase in .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.

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.

Figure 1 :Example 2 .
Figure 1: The ground state (top row) and time evolution of energy (bottom row) with fixed  = 4 and  = 1000, 3000, 6000 from left to right in Example 2.

Figure 2 :Figure 3 :
Figure 2: The ground state with fixed  = 100 and  = 4, 3 and  = 2 from left to right in Example 2.