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 and rth () 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.

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 much smaller than the critical condensation temperature , the properties of BECs are well described by a macroscopic wave function which obeys the Gross-Pitaevskii equation (GPE). More precisely, is solution to the dimensionless time-dependent GPE [3]:for , , . Here is the Laplace operator. Function 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 writewhere is the chemical potential of the condensate and a real function independent of time. Inserting the above expression into (1) gives the equationunder the normalization conditionAny eigenvalue can be computed from its corresponding eigenfunction bywhere the energy is given by

From mathematical point of view, the ground state is defined as the minimizer of over the unit sphere . 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 and set , . 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 . 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 [810].

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.

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:where and the computation domain is . The computation domain is discretized into grids described by the set , , , , where , , . Let , , be the number of spatial grid points in each spatial direction and let be the grid size, respectively. Let represent the approximate solution at the grid point .

Setting , , and , we get the discretization of (9) on the grid point as follows:Define the the following linear mapping in three-dimensional The linear mapping , and , are similarly defined.

By using a Taylor expansion, we get where . Omitting the small terms , we obtain the following approximation:Define the three-dimensional array: . 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 . Then we can derive an implicit integration factor method after approximating the integral. For example, the second order IIF takes the formWe 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 two-dimensional array where is a vector by fixing the last two indices ,  ,   = .

We define the matrices and , , whereWith the definition of matrices and , the exponential of linear mapping in the array representation can be written as matrix forms:The exponential of linear mapping and has similar array representations.

Because the linear mapping , , and satisfy the commuting property, we getDirect application of (18) and (19) to (15) results in the following array-representation integration factor CFDM with order :where . Then we apply normalization constraint to the solution :where .

For the 2D case the second order CFDM-AIF scheme is similar as (20) and (21):where and .

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

3.1. The Accuracy Test

The CFDM-AIF method (20) has the global accuracy 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 :with periodic boundary conditions. The exact solution of this equation is . The final computation time is 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 for CFDM-AIF method and 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, , the computation time is too long to be acceptable.

Remark. In the CFDM-AIF scheme (20), the matrix is which has order of magnitude smaller than the size of the matrix 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 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 .

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 aswhere is the distance between neighbor wells (lattice constant). The problem is solved on with grid number and time step .

Example 2. We first compute the ground state with different parameters and and fix . The initial data for ground state is chosen as Thomas-Fermi approximation: where .

For fixed , Figure 1 shows the ground state and time evolution of energy with , , and . 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 , Figure 2 shows the numerical results with , , and . 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.

Then we compute the first excited state by taking in -direction, in -direction, and in both -direction and -direction. Figure 3 visualizes the first excited states for , , and . The interior layers and multiscale structures are observed in the first excited states.

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, in (1). The initial data is chosen as for ground state and for the first excited state in -direction , for the first excited state in - and -direction , and for the first excited state in -, -, and -direction . The problem is solved on with grid number and time step .

Figure 4 shows the isosurface of ground state with , , and . From it, the number of peaks increases with the increase in . Figure 5 depicts the isosurface plots of the first excited state for at values of , , and . 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.

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.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

Rongpei Zhang and Guozhong Zhao’s work was supported by the National Nature Science Foundation of China (11261035). The authors are grateful to the referees for their invaluable suggestions and comments regarding the improvement of the paper.