Abstract

We study efficient continuation methods for computing the ground state solution of quasi-2D rotating dipolar Bose-Einstein condensates (BECs). First, the highly accurate spectral collocation method is used to discretize the governing Gross-Pitaevskii equation (GPE). Then, we modify the two-level continuation scheme for 3D dipolar BECs described in Jeng et al. (2014) to develop a single-parameter continuation method for quasi-2D rotating dipolar BECs, where the chemical potential is treated as the continuation parameter. Further, by adding the ratio of dipolar interaction strength to contact interaction strength as the second continuation parameter, we propose an efficient two-parameter continuation method which can effectively show the change of the ground-state vortex structures as the dipolar interaction strength gradually increases. Moreover, we also study linear stability analysis for the GPE. Sample numerical results on quasi-2D rotating dipolar BECs are reported.

1. Introduction

In 2005, the first dipolar Bose-Einstein condensate (BEC) was successfully produced by the group of Griesmaier [1] at the University of Stuttgart, in a gas of 52Cr atoms cooled to 700 nK. Later in 2011, Lu et al. [2] at Stanford University realized a dipolar BEC of 164Dy atoms. The next year, a dipolar BEC of 168Er atoms was also achieved in experiments at Innsbruck University [3]. These successful experiments provided new impetus for the theoretical and numerical studies of dipolar BECs.

In dipolar BECs, since 52Cr, 168Er, and 164Dy atoms have larger magnetic moments, the long-range dipolar interaction is nonnegligible, and it induces various interesting phenomena such as -wave collapse and expansion [4], roton spectrum [5, 6], self-bound dipolar droplets [7, 8], and the supersolid phase [9]. Moreover, vortices are also an important topic in BECs, which can be created by rotating the trap. In this paper, we will develop efficient numerical methods for computing the ground state solution of rotating dipolar BECs and then investigate the ground-state vortex structures.

Using the mean-field theory, the three-dimensional (3D) rotating dipolar BEC at absolute zero temperature is well described by the macroscopic wave function whose evolution is governed by the dimensionless Gross-Pitaevskii equation (GPE) [1012].

where is the imaginary unit; is the spatial variable; is the time variable; denotes the convolution operator with respect to the spatial variable; is the harmonic trapping potential with , , and being the trap frequencies in the -, -, and -direction, respectively; and are constants representing the strength of the contact (or short-range) and dipolar (or long-range) interaction, respectively; is the long-range dipolar interaction potential with the dipolar axis satisfying ; is angular velocity of the laser beam; and is the -component of the angular momentum with the momentum operator . Two important invariants of the GPE are the mass (or normalization) of the wave function

and the energy per particle

Using the equality [13]

where is the Dirac distribution function, , and , we can decouple the convolution into two terms

where the function is defined by

Substituting (5) into (1) and noticing (6), we can transform the GPE (1) for rotating dipolar BECs into a Gross-Pitaevskii-Poisson or Schrödinger-Poisson (SP) type system of the following form

and then the energy (3) can be rewritten as

In many physical experiments of dipolar BECs, the condensates are confined by an anisotropic harmonic potential. When we consider , where is a small parameter describing the strength of confinement in the -direction, the condensates shape will typically like a pancake, and the wave function can be factorized into [14, 15]

Substituting (10) into the 3D SP system (7)–(8) and the normalization constraint (2), we obtain a quasi-two-dimensional (quasi-2D) system of equations

with the constraint

where , , , with , and . To find the stationary states including ground state of a rotating dipolar BEC, we take the ansatz

where is the chemical potential and is a time-independent complex function. Substituting (14) into (11)–(13), we obtain the following nonlinear eigenvalue problem

with the constraint

During the last decade, some significant mathematical theories and numerical methods for the ground state of dipolar BECs have been reported in the literature [1317]. Bao et al. [13] transformed the 3D GPE for dipolar BECs into a Gross-Pitaevskii-Poisson type system and then proved rigorously the existence and uniqueness of its ground state. Using dimension reduction, Cai et al. [14] derived two mean-field equations for quasi-1D and quasi-2D dipolar BECs, respectively, and compared their ground state solutions with those of the 3D GPE. The existence and uniqueness of the ground state solutions of these two mean-field equations were established in [15]. In addition, to compute the ground state of dipolar BECs, various efficient numerical methods have been proposed, such as the backward Euler sine pseudospectral method [13], two-level continuation scheme [16], and normalized gradient flow method with nonuniform fast Fourier transform [17]. However, there were only few studies, e.g., [18], concerning the efficient numerical method for computing the ground state of rotating dipolar BECs. In this paper, we will focus on this issue. First, we modify the two-level continuation scheme [16] for 3D dipolar BECs to develop a single-parameter continuation method for quasi-2D rotating dipolar BECs. Next, by treating the chemical potential and the parameter corresponding to the strength of the dipolar interaction as the continuation parameters simultaneously, we propose an efficient two-parameter continuation method, which can trace the ground state solutions of quasi-2D rotating dipolar BECs with increasing the strength of the dipolar interaction. Thus, the change of the ground-state vortex structure with respect to the dipolar interaction can be easily obtained.

This paper is organized as follows. In Section 2, we briefly describe the spectral collocation method (SCM) for quasi-2D rotating dipolar BECs. In Section 3, we propose single- and two-parameter continuation algorithms for computing the ground state solution of rotating dipolar BECs. In Section 4, we study linear stability analysis for (11). The numerical results are reported in Section 5. Finally, some concluding remarks are given in Section 6.

2. A SCM for Quasi-2D Rotating Dipolar BECs

In this section, we describe the SCM using the sine functions as the basis functions for equations (15)–(17). First, since the wave function rapidly as , we replace the whole space in (15) by a sufficiently large domain and impose the zero Dirichlet boundary condition on the function . Next, since is a complex function, we let , where and are two real-valued functions. Then, equations (15)–(17) can be rewritten as

where and denote the partial derivatives of and with respect to and , respectively. Let be the trial function space, where . Then, and all the functions of satisfy the boundary conditions of (18). We choose uniform grids as the collocation points. The SCM for solving (18) is to find and such that the residuals vanish at the collocation points, that is,

Denote , , , , , and define the matrices

Let the vectors , , , and , where denotes the vectorization of a matrix formed by stacking the columns of into a single column vector. Then, and (21) can be expressed as a nonlinear system of equations involving parameter : where “” denotes the Hadamard product, stands for the -times Hadamard products of , and are defined by

where “” denotes the Kronecker product and is an diagonal matrix with as the diagonal entries. Moreover, the constraint condition (20) can be expressed as

To compute the matrix , we need to efficiently discretize the convolution . Based on the convolution theory, it is natural to consider using the Fourier transform to compute it. To be precise, since , where denotes the 2D Fourier transform operator, we have that

where denotes the 2D inverse Fourier transform operator, and the Fourier transform of the integral is given by [14]

Since the integrand in (28) decays exponentially fast, we can replace the domain of integration by a sufficiently large interval and then evaluate the definite integral via numerical quadratures, e.g., composite trapezoidal rule or Gaussian quadrature. Moreover, in practical computations, since we only know the values of and on the grid points of the domain , the Fourier transform and the inverse Fourier transform mentioned above are replaced by the discrete Fourier transform and the inverse discrete Fourier transform, respectively, and can be efficiently computed via the fast Fourier transform (FFT) and the inverse fast Fourier transform (IFFT).

In our numerical computations, we incorporate the SCM in the continuation method and use single- or two-parameter continuation algorithms, which are described in the next section, to compute the ground state solution of a quasi-2D rotating dipolar BECs. In the continuation method, we need the Jacobian matrix associated with , which is given as

where

and is a diagonal matrix with the components of the vector as the diagonal entries.

3. Efficient Two-Parameter Continuation Algorithms

Jeng et al. [16] proposed a two-level continuation scheme for computing the ground state solution of a 3D dipolar BEC, which is governed by a stationary state SP system. The two major differences between the system (15)–(17) for a quasi-2D rotating dipolar BEC and the stationary state SP system in [16] are that (i) in the former, the wave function is a complex-valued function, while in the latter, the wave function is a real-valued function. (ii) The function in (16) is defined by a convolution, while the corresponding one in the SP system is defined by a Poisson equation. However, the structures of these two systems are similar. We can modify the two-level continuation scheme in [16] to compute the ground state solution of a quasi-2D rotating dipolar BEC which is stated as follows.

First, we rewrite the complex system (15)–(17) into a real system, which is described in (18)–(20). Starting with , (18) is the governing equation for a rotating BEC. By treating the chemical potential as the continuation parameter, we use an efficient continuation algorithm [1922] to trace the ground state solution curve of (18) bifurcating at the first bifurcation point. The constraint condition (20) is set as the target point for the curve-tracking. Next, we compute the convolution in (19) to obtain an approximation of and compute for the next step continuation algorithm. Here, the operators , , and can be efficiently computed via the FFT and the IFFT. Then, we use the continuation algorithm again to trace the ground state solution curve of (18) bifurcating at the first bifurcation point until the constraint (20) is satisfied. Similarly, we use the FFT and IFFT again to compute the approximation of and for the next step continuation algorithm. We repeat the same process as above until two consecutive bifurcation points are close enough. Then, instead of tracing the ground state solution curve of (18), we perform the Newton method to solve (18) and (20) simultaneously and then update the approximation of . We keep repeating the process until the approximation of converges. Then, the desired ground state solution of (18)–(20) can be obtained. A detailed description of this continuation scheme is as follows.

Input:
the -th iterate of the iterative continuation scheme.
accuracy tolerance for two consecutive bifurcation points.
accuracy tolerance for .
, .
Step 1. Compute the first bifurcation point of (18) with .
while ( or ) do
  (i) .
  (ii) Treat as the continuation parameter and use the classical continuation algorithm to trace the ground state solution curve of (18) with until the constraint (20) is satisfied. Set this solution as .
  (iii) Compute and .
  (iv) Compute the first bifurcation point of (18) with .
end
Step 2. Set and , .
while ( or ) do
  (i) .
  (ii) Use the approximate solution as the initial guess, and perform the Newton method to solve (18) with and (20) simultaneously.
  (iii) Set the approximate solution obtained in (ii) by .
  (iv) Compute and .
end
The desired ground state solution .

Note that in Algorithm 1, the purpose of Step 1 is to obtain an appropriate approximate solution as the starting point for Step 2. And the purpose of Step 2 is to correct the function . Moreover, in the while loop of Step 1, when , we obtain the ground state solution of a rotating BEC and then set to continue the next iteration of the while loop. This means that we take as an approximate ground state solution of a rotating dipolar BEC and then use the iterative way to get a better approximation. Hence, although the initial choice , we actually use as an initial approximation of . Additionally, if a suitable approximation for the ground state solution of a rotating dipolar BEC can be easily known in advance (e.g., the Thomas-Fermi approximation is a suitable one when is large), we can choose instead of .

To investigate how the ground-state vortex structures are affected by the strength of the dipolar interaction, we can use Algorithm 1 to compute the ground state solutions for various values of . However, this is obviously time-consuming because we must repeatedly implement Algorithm 1. In the following, we will describe a two-parameter continuation algorithm which can effectively show the change of the ground-state vortex structures as gradually increases.

Let . In the first step of the two-parameter continuation algorithm, we fix and implement Algorithm 1 to compute the ground state solution of (18)–(20) with , where the chemical potential is treated as the continuation parameter. Set this solution as . In the second step, we add as the second continuation parameter and continue to compute the ground state solutions for other values of . Here, the ground state solution is used as the starting point. The procedure of tracing next accepted ground state solution of (18)–(20) with consists of two parts: (I)Find an approximate ground state solution for the next value of :

We fix in (18) and then perform the Euler predictor-Newton corrector process once to obtain a solution of (18) with under the normalization condition (20), where and are treated as the continuation parameters simultaneously. (II)Correct by the iterative way:

We repeatedly update the approximation of and solve (18) and (20) simultaneously until the approximation of converges. Then, the ground state solution of (18)–(20) with is obtained.

Using instead of and repeating (I) and (II), we can obtain the next accepted ground state solution of (18)–(20) with . Similarly, we can continue to trace the ground state solutions , , and so on, until the desired target value of is reached. A detailed description of this two-parameter continuation algorithm is as follows.

Input:
the initial value of the parameter .
the final value of the parameter .
accuracy tolerance for .
Step 1. Use Algorithm 1 to compute the ground state solution of (18)–(20) with . Set this solution as .
Step 2. Add as the second continuation parameter and trace the ground state solution curve of (18)–(20) until :
Use as the starting point and set .
whiledo
  (i) .
  (ii) Compute and .
  (iii) Treat and as the continuation parameters simultaneously, and perform the predictor-corrector process once to obtain a solution of (18) with under the normalization condition (20). Set this solution as .
  (iv) Set and compute , .
  (v) .
  (vi) Use as the initial guess and perform Newton’s method to solve (18) with and (20) simultaneously. Set this solution as .
  (vii) Compute and .
  (viii) Repeat the procedure (v)–(vii) until , then set which is the ground state solution of (18)–(20) with .
end

Note that in practical computations, we suggest choosing , which means that the dipolar interaction is absent, and then (18)–(20) are simplified to the system of equations for rotating BECs. Some efficient continuation algorithms [1921] have been proposed for computing the ground state solution of rotating BECs. We can use one of these algorithms instead of Algorithm 1 in Step 1. In addition, Step 2 has the advantage that it can trace the ground state solutions with increasing the value of . Thus, the evolution of the ground-state vortex structure with respect to the dipolar interaction can be easily observed.

4. Linear Stability Analysis

The aim of this section is to study linear stability analysis for (11). We impose an infinitesimal perturbation on the wave function in (14) by letting

where is the solution of (15) and is a complex function. Substituting (36) into (11) and using the approximation , we obtain

Then taking (15) into account and linearizing (i.e., neglecting terms of and ), we obtain the following linearized equation which models the development of the perturbation :

Since (39) includes both and , we use the technique described in [2325] to decompose the perturbation as

where and are complex functions, and is the unknown coefficient yet to be determined. Inserting (40) into (39), we obtain a linear eigenvalue problem of the form

Let be an operator matrix with and , . Then, (41) and (42) can be written in matrix form as

This eigenvalue problem has the following property.

Proposition 1. The complex eigenvalues of (45) must come in conjugate pairs as .

Proof. Let be a complex eigenvalue-eigenvector pair of (45). From (45), we have Take . Then, , , , and from (46), we have which means that is a complex eigenvalue-eigenvector pair of (45). Therefore, the complex eigenvalues of (45) always appear in conjugate pairs .

When all eigenvalues of (45) have negative or zero real parts; then, the corresponding perturbation in (40) decays exponentially or is bounded; hence, the wave function is linearly stable. In contrast, if at least one eigenvalue has a positive real part; then, grows exponentially; hence, the wave function is linearly unstable. Unfortunately, since (45) includes the unknown wave function , it is difficult to analyze the sign of the real part of its eigenvalues theoretically. However, when the wave function is obtained using numerical computation, we can compute all eigenvalues of (45) numerically and then discuss the linear stability of the wave function, see Example 4 for details.

5. Numerical Results

Algorithms 1 and 2 were implemented to compute the ground state solution of quasi-2D (rapidly) rotating dipolar BECs, where we used the SCM with (or for the rapidly rotating case) as the discretization method. In Examples 13, we chose , , in (15), and the computational domain . The accuracy tolerance for the Newton corrector is . We studied how the ground-state vortex structures were affected by the dipolar direction and the strength of the dipolar interaction. In Example 4, we discussed the linear stability analysis numerically. In Example 5, we considered the case of rapidly rotating dipolar BECs. All computations were executed on an Intel Core i7-2600K PC using Matlab language.

Example 1. We used Algorithm 1 to compute the ground state solution of quasi-2D rotating dipolar BECs (15) with , , and . We chose and in Algorithm 1. Table 1 lists the implementation details, i.e., the first bifurcation points and associated energy levels of (15), the values of , and the total execution time. Table 1 and Figure 1(a) show that we only need to trace the first solution curves twice, and the energy level of the ground state solution is . Figure 1(b) shows the contour plot of the density for the ground state solution, where the hexagonal arrangement of the vortices is clearly visible. Figure 1(c) shows the phase of for the ground state solution.

Example 2. Consider (15) with , , and . We implemented Algorithm 2 with and to compute the ground state solution of (15). Figure 2(a) displays the solution curve of the wave function in two-norm with respect to the chemical potential . The blue solid line and the horizontal red dashed line were obtained by implementing Steps 1 and 2 of Algorithm 2, respectively. Note that in Step 2 of Algorithm 2, we can obtain the ground state solutions for various values of . Figure 2(b) depicts the relationship between and . We observe that the chemical potential increases almost linearly with respect to . Figure 3 displays the contours of the ground state density function for , , , and . From this figure, we can see that (i) when (i.e., the dipolar interaction is absent), the vortex lattice forms a hexagonal structure. (ii) As the parameter increases, the number of vortices also increases and the vortices form straight lines that are parallel to the -axis. Table 2 lists the execution time of implementing Algorithm 2 and the number of iterations in the while loop (NIW) of Step 2. The total NIW of 14 means that from to , we obtain ground state solutions for 14 different values of . The average execution time for each value of is seconds. Hence, using the ground state solution for as the starting point, Step 2 of Algorithm 2 can effectively trace ground state solutions for other values of .

Example 3. We implemented Algorithm 2 to compute the ground state solution of (15) with , , and two different dipolar directions and . Figure 4 shows the contours of the ground state density function for some sample values of . From this figure, we can see that (i) as the parameter increases, the vortices form straight lines and then gradually pinned together to form vortex bands. (ii) For , the vortex band orientation is parallel to the -axis. (iii) For , the vortex band orientation makes an angle of 60 degrees with the -axis.

Example 4 (linear stability analysis). To study the linear stability of (11) numerically, we chose some data in Examples 23 and used the built-in function “eig” in Matlab to compute all eigenvalues of the corresponding equation (45). Figure 5 shows the discrete eigenvalues distribution of (45). From this figure, we can see that (i) when , there is one eigenvalue with a positive real part. (ii) When the angular velocity , some eigenvalues of (45) have positive real parts. Therefore, in these cases, the discrete ground state solutions of (11) are linearly unstable.

Example 5 (rapidly rotating dipolar BECs). In this example, we investigated how the ground-state vortex structures of rapidly rotating dipolar BECs were affected by , , and . To this end, we replaced the harmonic potential in (15) by the harmonic-plus-quartic potential [19, 26]

and set , , , , and consider two cases: and . (1). Algorithm 2 with and was implemented to compute the ground state solution of (15), where we considered two different dipolar directions: (a) and (b) (a). Figure 6(a) shows the contours of the ground state density function for , , , , and . From this figure, we can see that (i) when (i.e., ), six vortices are arranged in a hexagonal pattern. (ii) As the parameter increases, the vortices are gradually rearranged into two rows and form square lattices. And the orientation of the vortex lattice is parallel to the -axis. Moreover, the number of vortices also evolves gradually from 6 to 10(b). Figure 6(b) shows the contours of the ground state density function for , , , , and . We observe that as the parameter increases, the vortices are gradually arranged in three columns. And in the middle column, the number of vortices evolves gradually from 2 to 4(2). Algorithm 2 with and was implemented to compute the ground state solution of (15), where we chose . The execution time of Algorithm 2 and the NIW of Step 2 are listed in Table 3. From this table, we observe that using the ground state solution for as the starting point, Step 2 can effectively trace ground state solutions for other values of , and the average execution time for each value of is seconds. Figure 7 displays vortex lattice structures of the ground state for , , , and . From this figure, we can see that (i) when , one big vortex at the center is surrounded by eight small vortices in a ring. (ii) With the increasing of , the big vortex at the center is gradually split into two small vortices. (iii) When increases further, such as , the vortices are rearranged and aligned in three lines along the direction

From Examples 2, 3, and 5, we observe that as the parameter increases, that is, as the strength of the dipolar interaction increases, the vortices gradually form some vortex bands parallel to the direction . This phenomenon is mainly related to the dipolar term in (15), where involves the unknown wave function . To understand the effect of this dipolar term on the vortex configuration, we take , the ground state solution of (15) with , and numerically compute the function . Figure 8(a) depicts the contour plot of the density function for , while Figures 8(b) and 8(c) depict the contour plots of the function corresponding to with and , respectively. We can see that the contour plot of the function is an ellipse whose minor axis is perpendicular to the direction . This means that the dipolar term is anisotropic even though the density wave function is isotropic. Hence, as the strength of the dipolar interaction increases, the condensates will be squeezed in the direction perpendicular to . Due to this squeezing, the vortices will be rearranged and aligned parallel to the direction . Our numerical results confirm this phenomenon.

6. Conclusions

We have presented an efficient two-parameter continuation algorithm for computing the ground state solution of quasi-2D (rapidly) rotating dipolar BECs, where the SCM was used to discretize the corresponding GPE. The main advantage of this algorithm is that we only need to trace the solution curve once to obtain various ground state solutions associated with different strengths of the dipolar interaction. Thus, the change of the ground-state vortex structure with increasing the dipolar interaction strength can be easily observed. In addition, we have also studied linear stability analysis for the ground state of rotating dipolar BECs. Based on the numerical experiments reported in Section 5, we may give some concluding remarks as follows: (i) the ground-state vortex structure of rotating dipolar BECs is affected by the strength of the dipolar interaction and the dipolar direction. More specifically, as the strength of the dipolar interaction increases, the vortices are gradually aligned in lines along the dipolar direction. When the strength is large enough, some vortices on the vortex lines are pinned together to form vortex bands. Moreover, for the rapidly rotating case, we also observe that the big vortex at the center is gradually split into two small vortices with increasing the dipolar interaction. (ii) In Example 4, we have shown numerically that the discrete ground state solutions of quasi-2D dipolar BECs are linearly unstable, both with and without rotating term.

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This study was funded by Ministry of Science and Technology of R.O.C. (Taiwan) (MOST 103-2115-M-142-003).