Advances in Condensed Matter Physics

Volume 2015 (2015), Article ID 127580, 7 pages

http://dx.doi.org/10.1155/2015/127580

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

^{1}School of Sciences, Liaoning Shihua University, Fushun 113001, China^{2}School of Foreign Language, Liaoning Shihua University, Fushun 113001, China^{3}Faculty of Mathematics, Baotou Teachers College, Baotou 014030, China

Received 25 March 2015; Accepted 20 May 2015

Academic Editor: Sergei Sergeenkov

Copyright © 2015 Rongpei Zhang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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 *r*th () 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 [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 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.