# Efficient Continuation Methods for Computing Ground States of Quasi-2D Rotating Dipolar Bose-Einstein Condensates

Authors: B.-W. Jeng; P.-H. Su; M.-I. Char
Journal: Advances in Mathematical Physics
 (2022)

---

## 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.


---

## Body

## 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, since52Cr, 168Er, and 164Dy atoms have larger magnetic moments, the long-range dipolar interaction is nonnegligible, and it induces various interesting phenomena such as d-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Ψx,t whose evolution is governed by the dimensionless Gross-Pitaevskii equation (GPE) [10–12].
(1)i∂tΨx,t=−12Δ+Vx+βΨ2+λUdipx∗Ψ2−ωLzΨx,t,x∈ℝ3,t>0,wherei=−1 is the imaginary unit; x=x,y,z is the spatial variable; t is the time variable; ∗ denotes the convolution operator with respect to the spatial variable; Vx=1/2γx2x2+γy2y2+γz2z2 is the harmonic trapping potential with γx, γy, and γz being the trap frequencies in the x-, y-, and z-direction, respectively; β and λ are constants representing the strength of the contact (or short-range) and dipolar (or long-range) interaction, respectively; Udipx=3/4π1−3x⋅n2/x2/x3 is the long-range dipolar interaction potential with the dipolar axis n=n1,n2,n3∈ℝ3 satisfying n=n12+n22+n32=1; ω is angular velocity of the laser beam; and Lz=xpy−ypx=−ix∂y−y∂x is the z-component of the angular momentum L=x×P with the momentum operator P=−i∇=px,py,pz. Two important invariants of the GPE are the mass (or normalization) of the wave function
(2)Ψx,t2≔∫ℝ3Ψx,t2dx≡∫ℝ3Ψx,02dx=1,t≥0,and the energy per particle(3)EΨx,t≔∫ℝ312∇Ψ2+VxΨ2+β2Ψ4+λ2Udipx∗Ψ2Ψ2−ωReΨ¯LzΨdx≡EΨx,0,t≥0.Using the equality [13]
(4)Udipx=34π1−3x⋅n2/x2x3=−δx−3∂nn14πx,whereδx is the Dirac distribution function, ∂n=n⋅∇, and ∂nn=∂n∂n, we can decouple the convolution Udipx∗Ψ2 into two terms
(5)Udipx∗Ψx,t2=−Ψx,t2−3∂nnφx,t,where the functionφ is defined by
(6)φx,t≔14πx∗Ψx,t2,x∈ℝ3,t>0.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
(7)i∂tΨx,t=−12Δ+Vx+β−λΨ2−3λ∂nnφx,t−ωLzΨx,t,x∈ℝ3,t>0,(8)−Δφx,t=Ψx,t2,limx⟶∞φx,t=0,and then the energy (3) can be rewritten as
(9)EΨx,t=∫ℝ312∇Ψ2+VxΨ2+β−λ2Ψ4+3λ2∂n∇φ2−ωReΨ¯LzΨdx.In many physical experiments of dipolar BECs, the condensates are confined by an anisotropic harmonic potential. When we considerVx,y,z=1/2γx2x2+γy2y2+1/2ε4z2, where 0<ε≪1 is a small parameter describing the strength of confinement in the z-direction, the condensates shape will typically like a pancake, and the wave function Ψx,t can be factorized into [14, 15]
(10)Ψx,t=e−it/2ε2ψx,y,twεz,wherewεz=ε−1/2π−1/4e−z2/2ε2.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
(11)i∂tψx,y,t=−12Δ+V2x,y+β+λ3n32−12πεψ2−3λ2∂n⊥n⊥−n32Δφ2x,y,t−ωLzψx,y,t,x,y∈ℝ2,t≥0,(12)φ2x,y,t=U2x,y∗ψx,y,t2,U2x,y=122π3/2∫ℝe−s2/2x2+y2+ε2s2ds,with the constraint(13)ψx,y,t2=∫ℝ2ψx,y,t2dxdy=1,whereΔ=∂xx+∂yy, V2x,y=1/2γx2x2+γy2y2, n⊥=n1,n2, ∂n⊥=n⊥⋅∇⊥ with ∇⊥=∂x,∂y, and ∂n⊥n⊥=∂n⊥∂n⊥. To find the stationary states including ground state of a rotating dipolar BEC, we take the ansatz
(14)ψx,y,t=e−iμtϕx,y,t≥0,whereμ∈ℝ is the chemical potential and ϕx,y is a time-independent complex function. Substituting (14) into (11)–(13), we obtain the following nonlinear eigenvalue problem
(15)−12Δ+V2x,y+β+λ3n32−12πεϕ2−3λ2∂n⊥n⊥−n32Δφ~2x,y−ωLzϕx,y=μϕx,y,x,y∈ℝ2,(16)φ~2x,y=U2x,y∗ϕx,y2,with the constraint(17)ϕ2=∫ℝ2ϕx,y2dxdy=1.During the last decade, some significant mathematical theories and numerical methods for the ground state of dipolar BECs have been reported in the literature [13–17]. 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 Section2, 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 ϕx,y⟶0 rapidly as x,y⟶∞, we replace the whole space ℝ2 in (15) by a sufficiently large domain Ω=−ℓ,ℓ2 and impose the zero Dirichlet boundary condition on the function ϕx,y. Next, since ϕx,y is a complex function, we let ϕx,y=ux,y+ivx,y, where ux,y and vx,y are two real-valued functions. Then, equations (15)–(17) can be rewritten as
(18)−12Δu+V2x,yu+β+λ3n32−12πεu2+v2u−3λ2∂n⊥n⊥−n32Δφ~2x,yu−ωxvy−yvx=μu−12Δv+V2x,yv+β+λ3n32−12πεu2+v2v−3λ2∂n⊥n⊥−n32Δφ~2x,yv+ωxuy−yux=μvinΩ,ux,y=vx,y=0on∂Ω,(19)φ~2x,y=U2x,y∗u2+v2,(20)constraint:∫Ωu2+v2dxdy=1,whereux,vx and uy,vy denote the partial derivatives of u and v with respect to x and y, respectively. Let WN=spansjxsky:j,k=1,2,⋯,N−1 be the trial function space, where snx=sinnπx+ℓ/2ℓ. Then, dimWN=N−12 and all the functions of WN satisfy the boundary conditions of (18). We choose uniform grids xn,ym=−ℓ+n2ℓ/N,−ℓ+m2ℓ/N:n,m=1,2,⋯,N−1 as the collocation points. The SCM for solving (18) is to find uNx,y=∑j,k=1N−1ξjksjxsky∈WN and vNx,y=∑j,k=1N−1ηjksjxsky∈WN such that the residuals vanish at the collocation points, that is,
(21)−12ΔuNxn,ym+V2xn,ymuNxn,ym+β+λ3n32−12πεuN2xn,ym+vN2xn,ymuNxn,ym−3λ2∂n⊥n⊥−n32Δφ~2xn,ymuNxn,ym−ωxnvNyxn,ym−ymvNxxn,ym−μuNxn,ym=0,n,m=1,2,⋯,N−1,−12ΔvNxn,ym+V2xn,ymvNxn,ym+β+λ3n32−12πεuN2xn,ym+vN2xn,ymvNxn,ym−3λ2∂n⊥n⊥−n32Δφ~2xn,ymvNxn,ym+ωxnuNyxn,ym−ymuNxxn,ym−μvNxn,ym=0,n,m=1,2,⋯,N−1.Denotesij=sjxi, sij′=sj′xi, sij″=sj″xi, v2nm=V2xn,ym, qnm=∂n⊥n⊥−n32Δφ~2xn,ym, and define the N−1×N−1 matrices
(22)S=sij1≤i,j≤N−1,S′=sij′1≤i,j≤N−1,S″=sij″1≤i,j≤N−1,V2=v2nm1≤n,m≤N−1,Q=qnm1≤n,m≤N−1,U=ξjk1≤j,k≤N−1,V=ηjk1≤j,k≤N−1.Let the vectorsu=vecU, v=vecV, v2=vecV2, and q=vecQ, where vecW denotes the vectorization of a matrix W formed by stacking the columns of W into a single column vector. Then, u,v,v2,q∈ℝN−12 and (21) can be expressed as a nonlinear system of equations involving parameter μ:
(23)H1u,v,μ=−12Au+v2∘Bu+β+λ3n32−12πεBu2+Bv2∘Bu−3λ2q∘Bu−ωCv−μBu=0,(24)H2u,v,μ=−12Av+v2∘Bv+β+λ3n32−12πεBu2+Bv2∘Bv−3λ2q∘Bv+ωCu−μBv=0,where “∘” denotes the Hadamard product, wk=w∘⋯∘w stands for the k-times Hadamard products of w, and A,B,C∈ℝN−12×N−12 are defined by
(25)A=S⊗S″+S″⊗S,B=S⊗S,C=S′⊗diagx1,⋯,xN−1S−diagy1,⋯,yN−1S⊗S′,where “⊗” denotes the Kronecker product and diagd1,⋯,dN−1 is an N−1×N−1 diagonal matrix with d1,⋯,dN−1 as the diagonal entries. Moreover, the constraint condition (20) can be expressed as
(26)ℓ2uTu+vTv=1.To compute the matrixQ, we need to efficiently discretize the convolution φ~2. Based on the convolution theory, it is natural to consider using the Fourier transform to compute it. To be precise, since Fφ~2=FU2∗u2+v2=FU2⋅Fu2+v2, where F denotes the 2D Fourier transform operator, we have that
(27)φ~2=F−1FU2⋅Fu2+v2,whereF−1 denotes the 2D inverse Fourier transform operator, and the Fourier transform FU2 of the integral U2 is given by [14]
(28)FU2kx,ky=12π2∫ℝe−ε2s2/2kx2+ky2+s2ds.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 u and v 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 withH=H1,H2T, which is given as
(29)DHu,v,μ=DuH1u,v,μDvH1u,v,μDμH1u,v,μDuH2u,v,μDvH2u,v,μDμH2u,v,μ∈ℝ2N−12×2N−12+1,where(30)DuH1u,v,μ=−12A+diagv2B+β+λ3n32−12πεdiag3Bu2+Bv2B−3λ2diagqB−μB,(31)DvH1u,v,μ=β+λ3n32−12πεdiag2Bv∘BuB−ωC,(32)DμH1u,v,μ=−Bu,(33)DuH2u,v,μ=β+λ3n32−12πεdiag2Bu∘BvB+ωC,(34)DvH2u,v,μ=−12A+diagv2B+β+λ3n32−12πεdiagBu2+3Bv2B−3λ2diagqB−μB,(35)DμH2u,v,μ=−Bv,anddiagc is a diagonal matrix with the components of the vector c 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 φ~2 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 φ~2x,y=0, (18) is the governing equation for a rotating BEC. By treating the chemical potential μ as the continuation parameter, we use an efficient continuation algorithm [19–22] 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 φ~2 and compute ∂n⊥n⊥−n32Δφ~2 for the next step continuation algorithm. Here, the operators ∗, ∂n⊥n⊥, 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 φ~2 and ∂n⊥n⊥−n32Δφ~2 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 φ~2. We keep repeating the process until the approximation of φ~2 converges. Then, the desired ground state solution of (18)–(20) can be obtained. A detailed description of this continuation scheme is as follows.Algorithm 1: A single-parameter continuation algorithm for computing the ground state of a quasi-2D rotating dipolar BEC.
Input:i≔ the i-th iterate of the iterative continuation scheme.τ1≔ accuracy tolerance for two consecutive bifurcation points.τ2≔ accuracy tolerance for φ~2i,j.φ~20≔0, i≔0.Step 1. Compute the first bifurcation point0,0,μ1 of (18) with φ~2=φ~20.while (i=0 or μi−μi+1>τ1) do(i)i=i+1.(ii) Treatμ as the continuation parameter and use the classical continuation algorithm to trace the ground state solution curve of (18) with φ~2=φ~2i−1 until the constraint (20) is satisfied. Set this solution as ui∗,vi∗,μi∗.(iii) Computeφ~2i=U2∗ui∗2+vi∗2 and ∂n⊥n⊥−n32Δφ~2i.(iv) Compute the first bifurcation point0,0,μi+1 of (18) with φ~2=φ~2i.endStep 2. Setj=0 and ui,0∗,vi,0∗,μi,0∗=ui∗,vi∗,μi∗, φ~2i,0=φ~2i.while (j=0 or φ~2i,j−φ~2i,j−1∞>τ2) do(i)j=j+1.(ii) Use the approximate solutionui,j−1∗,vi,j−1∗,μi,j−1∗ as the initial guess, and perform the Newton method to solve (18) with φ~2=φ~2i,j−1 and (20) simultaneously.(iii) Set the approximate solution obtained in (ii) byui,j∗,vi,j∗,μi,j∗.(iv) Computeφ~2i,j=U2∗ui,j∗2+vi,j∗2 and ∂n⊥n⊥−n32Δφ~2i,j.endThe desired ground state solutionu∗,v∗,μ∗=ui,j∗,vi,j∗,μi,j∗.Note that in Algorithm1, 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 φ~2. Moreover, in the while loop of Step 1, when i=1, we obtain the ground state solution ϕ1∗=u1∗+iv1∗ of a rotating BEC and then set φ~21=U2∗ϕ1∗2 to continue the next iteration of the while loop. This means that we take ϕ1∗ 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 φ~20=0, we actually use φ~21=U2∗ϕ1∗2 as an initial approximation of φ~2. Additionally, if a suitable approximation ϕapp 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 φ~20=U2∗ϕapp2 instead of φ~20=0.To investigate how the ground-state vortex structures are affected by the strength of the dipolar interaction, we can use Algorithm1 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 γ=γ0 and implement Algorithm 1 to compute the ground state solution of (18)–(20) with λ=γ0β, where the chemical potential μ is treated as the continuation parameter. Set this solution as u0,v0,μ0. 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 y0=u0,v0,μ0,γ0=γ0 is used as the starting point. The procedure of tracing next accepted ground state solution y1=u1,v1,μ1,γ1 of (18)–(20) with λ=γ1β consists of two parts:
(I)
Find an approximate ground state solution for the next value ofγ:We fixφ~2=φ~20≔U2∗u02+v02 in (18) and then perform the Euler predictor-Newton corrector process once to obtain a solution of (18) with φ~2=φ~20 under the normalization condition (20), where μ and γ are treated as the continuation parameters simultaneously.
(II)
Correctφ~2 by the iterative way:We repeatedly update the approximation ofφ~2 and solve (18) and (20) simultaneously until the approximation of φ~2 converges. Then, the ground state solution y1 of (18)–(20) with λ=γ1β is obtained.Usingy1 instead of y0 and repeating (I) and (II), we can obtain the next accepted ground state solution y2=u2,v2,μ2,γ2 of (18)–(20) with λ=γ2β. Similarly, we can continue to trace the ground state solutions y3, y4, and so on, until the desired target value of γ is reached. A detailed description of this two-parameter continuation algorithm is as follows.Algorithm 2: A two-parameter continuation algorithm for computing the ground state of a quasi-2D rotating dipolar BEC.
Input:γ0≔ the initial value of the parameter γ.γend≔ the final value of the parameter γ.τ≔ accuracy tolerance for φ~2ij.Step 1. Use Algorithm1 to compute the ground state solution of (18)–(20) with λ=γ0β. Set this solution as u0,v0,μ0.Step 2. Addγ as the second continuation parameter and trace the ground state solution curve of (18)–(20) until γ≥γend:Useu0,v0,μ0,γ0=u0,v0,μ0,γ0 as the starting point and set i=0.whileγi<γenddo(i)i=i+1.(ii) Computeφ~2i−1=U2∗ui−12+vi−12 and ∂n⊥n⊥−n32Δφ~2i−1.(iii) Treatμ and γ as the continuation parameters simultaneously, and perform the predictor-corrector process once to obtain a solution of (18) with φ~2=φ~2i−1 under the normalization condition (20). Set this solution as ui,0,vi,0,μi,0,γi,0.(iv) Setj=0 and compute φ~2i,0=U2∗ui,02+vi,02, ∂n⊥n⊥−n32Δφ~2i,0.(v)j=j+1.(vi) Useui,j−1,vi,j−1,μi,j−1,γi,j−1 as the initial guess and perform Newton’s method to solve (18) with φ~2=φ~2i,j−1 and (20) simultaneously. Set this solution as ui,j,vi,j,μi,j,γi,j.(vii) Computeφ~2i,j=U2∗ui,j2+vi,j2 and ∂n⊥n⊥−n32Δφ~2i,j.(viii) Repeat the procedure (v)–(vii) untilφ~2i,j−φ~2i,j−1∞<τ, then set ui,vi,μi,γi=ui,j,vi,j,μi,j,γi,j which is the ground state solution of (18)–(20) with λ=γiβ.endNote that in practical computations, we suggest choosingγ0=0, 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 [19–21] 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 rx,y,t on the wave function ψx,y,t in (14) by letting
(36)ψx,y,t=e−iμtϕx,y+rx,y,t,whereμ,ϕx,y is the solution of (15) and rx,y,t is a complex function. Substituting (36) into (11) and using the approximation φ2=U2∗ψ2≈U2∗ϕ2=φ~2, we obtain
(37)e−iμtμϕ+r+i∂r∂t≈e−iμt−12Δ+V2+β+λ3n32−12πεϕ2+ϕr¯+ϕ¯r+r2−3λ2∂n⊥n⊥−n32Δφ~2−ωLzϕ+r,(38)⇒i∂r∂t≈−12Δ+V2+β+λ3n32−12πεϕ2−3λ2∂n⊥n⊥−n32Δφ~2−ωLz−μϕ+−12Δ+V2−3λ2∂n⊥n⊥−n32Δφ~2−ωLz−μr+β+λ3n32−12πε2ϕ2r+ϕ2r¯+2ϕr2+ϕ¯r2+r2r.Then taking (15) into account and linearizing (i.e., neglecting terms of Ork and k≥2), we obtain the following linearized equation which models the development of the perturbation rx,y,t:
(39)i∂r∂t=−12Δ+V2x,y−3λ2∂n⊥n⊥−n32Δφ~2x,y−ωLz−μr+β+λ3n32−12πε2ϕ2r+ϕ2r¯.Since (39) includes both r and r¯, we use the technique described in [23–25] to decompose the perturbation rx,y,t as
(40)rx,y,t=px,yeαt+qx,y¯eα¯t,wherepx,y and qx,y 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
(41)−12Δ+V2x,y+2β+λ3n32−12πεϕ2−3λ2∂n⊥n⊥−n32Δφ~2x,y−ωLz−μpx,y+β+λ3n32−12πεϕ2qx,y=iαpx,y,(42)−12Δ+V2x,y+2β+λ3n32−12πεϕ2−3λ2∂n⊥n⊥−n32Δφ~2x,y−ωLz¯−μqx,y+β+λ3n32−12πεϕ2¯px,y=−iαqx,y.LetAϕ=f11ϕf12ϕf12ϕ¯f22ϕ¯ be an operator matrix with
(43)f11ϕ=−12Δ+V2x,y+2β+λ3n32−12πεϕ2−3λ2∂n⊥n⊥−n32Δφ~2x,y−ωLz−μ,(44)f12ϕ=β+λ3n32−12πεϕ2,and B=diagi,−i∈ℂ2×2, r=px,y,qx,yT. Then, (41) and (42) can be written in matrix form as
(45)Aϕr=αBr.This eigenvalue problem has the following property.Proposition 1.
The complex eigenvalues of (45) must come in conjugate pairs as α,α¯.Proof.
Letα,r be a complex eigenvalue-eigenvector pair of (45). From (45), we have
(46)Aϕ¯r¯=α¯B¯r¯.
TakeQ=0110. Then, Q−1=Q, QAϕQ−1=Aϕ¯, QBQ−1=B¯, and from (46), we have
(47)AϕQ−1r¯=α¯BQ−1r¯,
which means thatα¯,Q−1r¯ 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 rx,y,t 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, rx,y,t 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

Algorithms1 and 2 were implemented to compute the ground state solution of quasi-2D (rapidly) rotating dipolar BECs, where we used the SCM with N=51 (or N=65 for the rapidly rotating case) as the discretization method. In Examples 1–3, we chose V2x,y=1.003x2+0.997y2/2, ε=1/5, ω=0.9 in (15), and the computational domain Ω=−12,122. The accuracy tolerance for the Newton corrector is 5×10−7. 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 CoreTM i7-2600K PC using Matlab language.Example 1.
We used Algorithm1 to compute the ground state solution of quasi-2D rotating dipolar BECs (15) with n=1,0,1/2, β=500, and λ=200. We chose τ1=0.1 and τ2=10−6 in Algorithm 1. Table 1 lists the implementation details, i.e., the first bifurcation points μi and associated energy levels μi∗ of (15), the values of φ~2i,j−φ~2i,j−1∞, 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 μ∗=6.796239445. Figure 1(b) shows the contour plot of the density ϕ2 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 n=1,0,1/2, β=500, and λ=γβ. We implemented Algorithm 2 with γ0=0 and γend=2.2 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 ϕ2 for γ=0, 0.5123, 1.1017, and 2.0856. From this figure, we can see that (i) when γ=0 (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 x-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 γ=0 to γ=2.0856, we obtain ground state solutions for 14 different values of γ. The average execution time for each value of γ is 551.28/14≈39.38 seconds. Hence, using the ground state solution for γ=0 as the starting point, Step 2 of Algorithm 2 can effectively trace ground state solutions for other values of γ.Example 3.
We implemented Algorithm2 to compute the ground state solution of (15) with β=500, λ=γβ, and two different dipolar directions n=0,1,1/2 and n=1,1,1/3. Figure 4 shows the contours of the ground state density function ϕ2 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 n=0,1,1/2, the vortex band orientation is parallel to the y-axis. (iii) For n=1,1,1/3, the vortex band orientation makes an angle of 60 degrees with the x-axis.Example 4 (linear stability analysis).
To study the linear stability of (11) numerically, we chose some data in Examples 2–3 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 ω=0, there is one eigenvalue with a positive real part. (ii) When the angular velocity ω=0.9, 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ω, n, and λ. To this end, we replaced the harmonic potential V2x,y=1.003x2+0.997y2/2 in (15) by the harmonic-plus-quartic potential [19, 26]
(48)V~2x,y=1.003x2+0.997y22+1.003x2+0.997y224,Table 1
The first bifurcation points and associated energy levels of (15) withV2x,y=1.003x2+0.997y2/2, ε=1/5, n=1,0,1/2, β=500, λ=200, and ω=0.9.
Step 1 of Algorithm 1Step 2 of Algorithm 1iμiμi−μi−1μi∗jμ2,j∗φ~22,j−φ~22,j−1∞10.999994079—6.94881042116.8098209323.05e-0420.8131896990.1868043796.83548357726.8107178134.63e-0530.8008194360.012370264—36.8102291104.02e-0546.8093999743.87e-0556.8082228093.97e-0566.8066319554.82e-0576.8046212075.55e-0586.8024252775.53e-0596.8004561244.69e-05106.7989653563.57e-05116.7979543512.64e-05126.7973086161.95e-05136.7969066621.40e-05146.7966574709.84e-06156.7965016816.89e-06166.7964028914.80e-06176.7963392583.33e-06186.7962976752.30e-06196.7962701681.58e-06206.7962518001.09e-06216.7962394457.44e-07μ∗=6.796239445Total execution time: 723.39 secFigure 1
The solution curves, the contour plot ofϕ2, and the phase of ϕ for the ground state solution of (15), where V2x,y=1.003x2+0.997y2/2, ε=1/5, n=1,0,1/2, β=500, λ=200, and ω=0.9.
(a)
Solution curves(b)
ϕ2(c)
Phase ofϕFigure 2
Implementing Algorithm2 for the ground state solution branch of (15), where V2x,y=1.003x2+0.997y2/2, ε=1/5, n=1,0,1/2, ω=0.9, β=500, and λ=γβ.
(a)
Solution curve(b)
Relationship betweenγ and μFigure 3
The contours of the ground state density functionϕ2 of (15), where V2x,y=1.003x2+0.997y2/2, ε=1/5, n=1,0,1/2, ω=0.9, β=500, and λ=γβ.Table 2
The execution time (in seconds) for computing the ground state solution of Example2 by implementing Algorithm 2. NIW: number of iterations in the while loop of Step 2.
γNIWTime (s)Step 10—349.10Step 20–0.51234158.570.5123–1.1017485.511.1017–2.08566307.20Total time of Step 2: 551.28Total14900.38Figure 4
The contours of the ground state density functionϕ2 of (15), where V2x,y=1.003x2+0.997y2/2, ε=1/5, ω=0.9, β=500, and λ=γβ.
(a)
n=0,1,1/2(b)
n=1,1,1/3Figure 5
The discrete eigenvalues distribution of (45) for quasi-2D (rotating) dipolar BEC, where V2x,y=1.003x2+0.997y2/2, ε=1/5, β=500, and λ=γβ.
(a)
ω=0, n=1,0,1/2, γ=0.8(b)
ω=0.9, n=1,0,1/2, γ=1.1017(c)
ω=0.9, n=0,1,1/2, γ=1.0848(d)
ω=0.9, n=1,1,1/3, γ=0.9722and setε=1/5, β=500, λ=γβ, Ω=−6,62, and consider two cases: ω=1.7 and ω=1.9.
(1)
ω=1.7. Algorithm 2 with γ0=0 and γend=1.65 was implemented to compute the ground state solution of (15), where we considered two different dipolar directions: (a) n=1,0,1/2 and (b) n=0,1,1/2(a)
n=1,0,1/2. Figure 6(a) shows the contours of the ground state density function ϕ2 for γ=0, 0.3410, 0.7196, 1.1740, and 1.6068. From this figure, we can see that (i) when γ=0 (i.e., λ=0), 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 x-axis. Moreover, the number of vortices also evolves gradually from 6 to 10(b)
n=0,1,1/2. Figure 6(b) shows the contours of the ground state density function ϕ2 for γ=0, 0.3422, 1.0160, 1.3246, and 1.6460. 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)
ω=1.9. Algorithm 2 with γ0=0 and γend=1.5 was implemented to compute the ground state solution of (15), where we chose n=1,0,1/2. 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 γ=0 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 4958.92/27≈183.66 seconds. Figure 7 displays vortex lattice structures of the ground state for γ=0, 0.7928, 0.8076, and 1.4735. From this figure, we can see that (i) when γ=0, 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 γ=1.4735, the vortices are rearranged and aligned in three lines along the direction n⊥Figure 6
The contours of the ground state density functionϕ2 of (15), where V2x,y=1.003x2+0.997y2/2+1.003x2+0.997y22/4, ε=1/5, ω=1.7, β=500, and λ=γβ.
(a)
n=1,0,1/2(b)
n=0,1,1/2Table 3
The execution time (in seconds) for computing the ground state solution of Example5(2) by implementing Algorithm 2. NIW: number of iterations in the while loop of Step 2.
γNIWTime (s)Step 10—1907.07Step 20–0.8076162767.420.8076–1.4735112191.50Total time of Step 2: 4958.92Total276865.99Figure 7
The contours of the ground state density functionϕ2 of (15), where V2x,y=1.003x2+0.997y2/2+1.003x2+0.997y22/4, ε=1/5, n=1,0,1/2, ω=1.9, β=500, and λ=γβ.From Examples2, 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 n⊥. This phenomenon is mainly related to the dipolar term −3λ/2∂n⊥n⊥−n32Δφ~2 in (15), where φ~2=U2∗ϕ2 involves the unknown wave function ϕ. To understand the effect of this dipolar term on the vortex configuration, we take ϕ=ϕg, the ground state solution of (15) with λ=0, and numerically compute the function ∂n⊥n⊥−n32Δφ~2. Figure 8(a) depicts the contour plot of the density function ϕg2 for ω=0, while Figures 8(b) and 8(c) depict the contour plots of the function ∂n⊥n⊥−n32Δφ~2 corresponding to ϕ=ϕg with n=1,0,1/2 and 0,1,1/2, respectively. We can see that the contour plot of the function ∂n⊥n⊥−n32Δφ~2 is an ellipse whose minor axis is perpendicular to the direction n⊥. This means that the dipolar term is anisotropic even though the density wave function ϕg2 is isotropic. Hence, as the strength of the dipolar interaction increases, the condensates will be squeezed in the direction perpendicular to n⊥. Due to this squeezing, the vortices will be rearranged and aligned parallel to the direction n⊥. Our numerical results confirm this phenomenon.Figure 8
The contour of the density functionϕg2 of the ground state of (15) with ω=0 and λ=0, and the contours of the function ∂n⊥n⊥−n32Δφ~2 corresponding to ϕg.
(a)
ϕg2(b)
∂n⊥n⊥−n32Δφ~2, n=1,0,1/2(c)
∂n⊥n⊥−n32Δφ~2, n=0,1,1/2


## 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 Section5, 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.



---

*Source: 9485124-2022-08-04.xml*