Xavier Antoine
email: xavier.antoine@univ-lorraine.fr

Qinglin Tang
email: qinglintang@163.com

Yong Zhang

A Preconditioned Conjugated Gradient Method for Computing Ground States of Rotating Dipolar Bose-Einstein Condensates via Kernel Truncation Method for Dipole-Dipole Interaction Evaluation

Keywords: rotating dipolar BEC, Dipole-Dipole Interaction, Preconditioned Conjugate Gradient method, Kernel Truncation Method, ground state . . . . . . . . .

published or not. The documents may come

A Preconditioned Conjugated Gradient method for computing ground states of rotating dipolar Bose-Einstein condensates via Kernel Truncation method for Dipole-Dipole Interaction evaluation 1 Introduction

The successful realization of Bose-Einstein condensates of gases of 52 Cr [START_REF] Griesmaier | Bose-Einstein condensation of Chromium[END_REF], 164 Dy [START_REF] Kadau | Observing the Rosensweig instability of a quantum ferrofluid[END_REF] and 168 Er [START_REF] Aikawa | Bose-Einstein Condensation of Erbium[END_REF] provides the possibility to study and probe novel interesting phenomenon with dipolar quantum gas. Different from the early BECs of ultra cold atomic gases whose properties are mainly governed by isotropic and short-range interatomic interactions [START_REF] Pitaevskii | Bose-Einstein Condensation[END_REF], in dipolar BEC, the magnetic/electric dipole-dipole interatomic interaction is anisotropic and long-range, and it brings in unique phenomenona, such as roton-maxon spectrum [START_REF] Lahaye | The physics of dipolar bosonic quantum gases[END_REF][START_REF] Ronen | Bogoliubov modes of a dipolar condensate in a cylindrical trap[END_REF], vortex lattice patterns [START_REF] Martin | Vortices and vortex lattices in quantum ferrofluids[END_REF], and the self-bound droplet state [START_REF] Ferrier-Barbut | Observation of quantum droplets in a strongly Dipolar Bose gas[END_REF][START_REF] Schmitt | Self-bound droplets of a dilute magnetic quantum liquid[END_REF]. There is everlasting enthusiasm in studying both the ground state [START_REF] Antoine | GPELab, a Matlab toolbox to solve Gross-Pitaevskii equations I: Computation of stationary solutions[END_REF][START_REF] Antoine | GPELab, a Matlab toolbox to solve Gross-Pitaevskii equations II: dynamics and stochastic simulations[END_REF][START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF][START_REF] Bao | Efficient numerical methods for computing ground states and dynamics of dipolar Bose-Einstein condensates[END_REF][START_REF] Danaila | Computation of ground states of the Gross-Pitaevskii functional via riemannian optimization[END_REF][START_REF] Jiang | Ground state of the dipolar Bose-Einstein condensate[END_REF][START_REF] Santos | Bose-Einstein condensation in trapped dipolar gases[END_REF][START_REF] Wu | A regularized Newton method for computing ground states of Bose-Einstein condensates[END_REF][START_REF] Yi | Trapped condensates of atoms with dipole interactions[END_REF][START_REF] Zeng | Efficiently computing vortex lattices in rapid rotating Bose-Einstein condensates[END_REF] and dynamics [START_REF] Bao | A simple and efficient numerical method for computing dynamics of rotating dipolar Bose-Einstein condensation via a rotating Lagrange coordinate[END_REF][START_REF] Carles | On the Gross-Pitaevskii equation for trapped dipolar quantum gases[END_REF][START_REF] Huang | Numerical simulation of trapped dipolar quantum gases: collapse studies and vortex dynamics[END_REF][START_REF] Lahaye | D-wave collapse and explosion of a dipolar Bose-Einstein condensate[END_REF][START_REF] Parker | Structure formation during the collapse of a dipolar atomic Bose-Einstein condensate[END_REF] of dipolar BECs.

At temperatures T much smaller than the critical temperature T c , the properties of BEC with long-range dipole-dipole interactions (DDI) are well described by the macroscopic complex-valued wave function ψ(x,t) whose evolution is governed by the celebrating three-dimensional (3D) Gross-Pitaevskii equation (GPE) with a DDI term. Moreover, the 3D GPE can be reduced to an effective two-dimensional (2D) version if the external trapping potential is strongly confined in the z-direction [START_REF] Bao | Efficient numerical methods for computing ground states and dynamics of dipolar Bose-Einstein condensates[END_REF][START_REF] Cai | Mean-field regime of trapped dipolar Bose-Einstein condensates in one and two dimensions[END_REF]. In a unified way, the dimensionless GPE with a DDI term in d-dimensions (d = 2 or 3) for modeling a dipolar BEC reads as [START_REF] Bao | Gross-Pitaevskii-Poisson equations for dipolar Bose-Einstein condensate with anisotropic confinement[END_REF][START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF][START_REF] Bao | A simple and efficient numerical method for computing dynamics of rotating dipolar Bose-Einstein condensation via a rotating Lagrange coordinate[END_REF][START_REF] Óral | Bose-Einstein condensation with magnetic dipole-dipole forces[END_REF][START_REF] Yi | Trapped condensates of atoms with dipole interactions[END_REF]:

i∂ t ψ(x,t) = - 1 2
∇ 2 +V(x)+ β|ψ| 2 +λΦ(x,t)-ωL z ψ(x,t), x ∈ R d , t > 0, (1.1)

Φ(x,t) = U dip * |ψ| 2 (x,t), x ∈ R d , t ≥ 0, (1.2) ψ(x,t = 0) = ψ 0 (x), x ∈ R d , (1.3)
where t is the time variable, x = (x,y) T ∈ R 2 or x = (x,y,z) T ∈ R 3 , * represents the convolution operator with respect to the spatial variable, L z = -i(x∂ y -y∂ x ) = -i∂ θ is the z-component of the angular momentum and ω represents the rotating frequency. The dimensionless constant β describes the strength of the short-range two-body interactions in a condensate (positive for repulsive interaction, and resp. negative for attractive interaction), V(x) is a given real-valued external trapping potential which is determined by the type of system under investigation. In most BEC experiments, a harmonic potential is chosen to trap the condensate, i.e.,

V(x) = 1 2 γ 2 x x 2 +γ 2 y y 2 , d = 2, γ 2 x x 2 +γ 2 y y 2 +γ 2 z z 2 , d = 3, (1.4)
where γ x > 0, γ y > 0 and γ z > 0 are dimensionless constants proportional to the trapping frequencies in x-, yand z-direction, respectively. Moreover, λ is a constant characterizing the strength of DDI and U dip (x) is the long-range DDI potential. In 3D, U dip (x) reads as

U dip (x) = 3 4π|x| 3 1- 3(x•n) 2 |x| 2 = -δ(x)-3∂ nn 1 4π|x| , x ∈ R 3 , (1.5)
with n = (n 1 ,n 2 ,n 3 ) T a given unit vector i.e. |n(t)| = n 2 1 +n 2 2 +n 2 3 = 1, representing the dipole axis (or dipole moment), ∂ n = n•∇ and ∂ nn = ∂ n (∂ n ), while in 2D, it is defined as [START_REF] Bao | Gross-Pitaevskii-Poisson equations for dipolar Bose-Einstein condensate with anisotropic confinement[END_REF][START_REF] Cai | Mean-field regime of trapped dipolar Bose-Einstein condensates in one and two dimensions[END_REF]

U dip (x) = - 3 2 ∂ n ⊥ n ⊥ -n 2 3 ∇ 2 ⊥ 1 2π|x| , x ∈ R 2 , (1.6)
where

∇ ⊥ = (∂ x ,∂ y ) T , n ⊥ = (n 1 ,n 2 ) T , ∂ n ⊥ = n ⊥ •∇ ⊥ and ∂ n ⊥ n ⊥ = ∂ n ⊥ (∂ n ⊥ ).
In fact, for smooth densities ρ(x) := |ψ(x)| 2 , the DDI potential can be reformulated via the Coulomb potential [START_REF] Bao | Efficient numerical methods for computing ground states and dynamics of dipolar Bose-Einstein condensates[END_REF][START_REF] Exl | Accurate and efficient computation of nonlocal potentials based on Gaussian-sum approximation[END_REF] whose convolution kernel is

U cou (x)= 1 2 d-1 |x| . The 3D DDI (1.5) is written as follows Φ(x) = -ρ-3 ∂ n ∂ n 1 4π|x| * ρ = -ρ-3 1 4π|x| * (∂ n ∂ n ρ), x ∈ R 3 , (1.7)
and the 2D DDI (1.6) is

Φ(x) = - 3 2 1 2π|x| * [ ∂ n ⊥ n ⊥ -n 2 3 ∇ 2 ⊥ ρ], x ∈ R 2 .
(1.8)

The time dependent GPE (1.1)-(1.3) conserves two important quantities: the total mass (or normalization) of the wave function

N (t) := ψ(x,t) 2 ≡ N (0), (1.9)
and the energy per particle

E (ψ(x,t)) = R d 1 2 |∇ψ| 2 +V(x)|ψ| 2 + λ 2 Φ|ψ| 2 -ωψ * L z ψ dx ≡ E (ψ(x,0)), t ≥ 0. (1.10)
The ground state φ g is defined as follows:

φ g = argmin φ∈S E (φ), where S := {φ(x) | φ 2 := R d |φ(x)| 2 dx = 1, E (φ) < ∞}. (1.11)
There have been extensive mathematical and numerical studies on the ground state of dipolar BEC, for an incomplete list we refer to [11-13, 19, 20, 23, 34, 55].

To compute the ground state, the main difficulties lie in (1) efficient ground state solver especially for fast rotating speed ω and (2) accurate DDI fast solver. During the last decade, there have been several numerical methods to compute the ground state of non-dipolar rotating BEC, which includes the normalized gradient flow method [START_REF] Antoine | Acceleration of the imaginary time method for spectrally computing the stationary states of Gross-Pitaevskii equations[END_REF][START_REF] Bao | Efficient and spectrally accurate numerical methods for computing ground and first excited states in Bose-Einstein condensates[END_REF][START_REF] Bao | Computing the ground state solution of Bose-Einstein condensates by a normalized gradient flow[END_REF], Sobolev gradient method [START_REF] Danaila | A finite element method with mesh adaptivity for computing vortex states in fast-rotating Bose-Einstein condensates[END_REF], projection gradient flow method [START_REF] Wang | A projection gradient method for computing ground state of spin-2 Bose-Einstein condensates[END_REF], preconditioned gradient flow method [START_REF] Antoine | Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose-Einstein Condensates[END_REF] and Conjugate Gradient method [START_REF] Ronen | Bogoliubov modes of a dipolar condensate in a cylindrical trap[END_REF]. We refer to [START_REF] Bao | Efficient numerical methods for computing ground states and dynamics of dipolar Bose-Einstein condensates[END_REF][START_REF] Bao | Accurate and efficient numerical methods for computing ground states and dynamics of dipolar Bose-Einstein condensates via the nonuniform FFT[END_REF] for their adaptations to dipolar BEC case. The efficiency of all these methods suffers severely for fast rotating case. Recently, Antoine et al. [START_REF] Antoine | Efficient spectral computation of the stationary states of rotating Bose-Einstein condensates by preconditioned nonlinear conjugate gradient methods[END_REF] proposed a preconditioned conjugate gradient method (PCG) for the rotating non-dipolar BEC, and its performance for the fast rotating case is stable and extraordinarily efficient. The preconditioner introduced there incorporates almost all the essential information, therefore iteration number drops down dramatically and is independent of domain and mesh size. It is imperative to adapt such mechanism to rotating dipolar BEC.

The second challenge is the accurate and efficient evaluation of nonlocal potentials. As is shown earlier, the DDI can be computed via the Coulomb potential. Starting from the partial differential/pseudo-differential equation reformulation [START_REF] Bao | Dimension reduction of the Schrödinger equation with Coulomb and anisotropic confining potentials[END_REF], on bounded rectangular domains with Dirichlet boundary conditions, the Discrete Sine Transform (DST) method applies directly [START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF][START_REF] Bao | Efficient numerical methods for computing ground states and dynamics of dipolar Bose-Einstein condensates[END_REF]. However, the accuracy of DST method depends linearly on the domain size. Therefore, it requires a quite large domain and presents significant challenges in both memory and CPU time, especially for 3D problems. In the last three years, Jiang et al. [START_REF] Jiang | Fast and accurate evaluation of nonlocal Coulomb and dipole-dipole interactions via the nonuniform FFT[END_REF], Zhang et al. [START_REF] Exl | Accurate and efficient computation of nonlocal potentials based on Gaussian-sum approximation[END_REF], Greengard et al. [START_REF] Vico | Fast convolution with free-space Green's functions[END_REF] respectively proposed the NonUniform Fast Fourier Transform (NUFFT) solver and Gaussian-Summation (GauSum) solver, Kernel Truncation Method (KTM). They are accurate (close to machine precision) and fast (with a O(Nlog(N)) complexity, where N is the total number of grid points) and free of boundary conditions for Φ. We refer to [START_REF] Antoine | On the ground states and dynamics of space fractional nonlinear Schrödinger/Gross-Pitaevskii equations with rotation term and nonlocal nonlinear interactions[END_REF][START_REF] Bao | Computing the ground state and dynamics of the nonlinear Schrödinger equation with nonlocal interactions via the nonuniform FFT[END_REF][START_REF] Bao | Accurate and efficient numerical methods for computing ground states and dynamics of dipolar Bose-Einstein condensates via the nonuniform FFT[END_REF][START_REF] Tang | A robust and efficient numerical method to compute the dynamics of the rotating two-component dipolar Bose-Einstein condensates[END_REF] for extensions and applications of such methods in the context of NLSE. Among them, KTM is the most simple one and has been widely adopted by the phyiscs community to compute the long-range interaction, e.g., [START_REF] Ronen | Bogoliubov modes of a dipolar condensate in a cylindrical trap[END_REF][START_REF] Rozzi | Exact Coulomb cutoff technique for supercell calculations[END_REF]. In this paper, we adapt KTM to 2D/3D DDI with further simplification, especially in efficiency, and generalize it to more general free-space convolutiontype potentials. The combination of PCG and KTM shall provide an accurate and efficient scheme for rotating dipolar BEC.

The main objectives of this paper are threefold:

1. to adapt KTM to 2D/3D DDI via Coulomb potential (1.7)-(1.8), to improve the efficiency by multidimensional discrete convolution with a twofold zero-padded wave function (density) and to discuss its extensions to related quasi-2D convolutiontype potential encountered in dimension reduction of dipolar BEC [START_REF] Cai | Mean-field regime of trapped dipolar Bose-Einstein condensates in one and two dimensions[END_REF].

2. to adapt PCG for computing ground states using KTM as DDI solver, to present de-tailed optimal θ n selection and preconditioners construction, and to verify numerically its performance, especially its efficiency with/without dependence of domain and mesh size.

3. to apply our method to study the ground state of 2D and 3D rotating dipolar BEC with different setups, and to investigate the vortex patterns.

The rest of the paper is organized as follows: In Section 2, we present a detailed introduction of KTM, including its derivation, implementation and acceleration via multidimensional discrete convolution, and discuss its extension to other related nonlocal convolution potentials. In Section 3, we propose an accurate and efficient scheme by combing PCG and KTM, and give three different preconditioners. Extensive numerical results are shown in Section 4 to confirm the performance of our method, together with interesting vortex patterns in 2D and 3D rotating dipolar BEC. Finally, conclusions are drawn in Section 5.

Dipole-Dipole Interaction evaluation via KTM

In presence of a confining potential V(x), the wave function ψ(x) (also the density ρ) is usually smooth and decays exponentially fast. Therefore it is quite convenient to approximate ψ on a truncated rectangular domain D instead of the whole space R d . As is shown by (1.7)-(1.8), the DDI computation boils down to Coulomb potential evaluation. Therefore, in this section, we first give a review of KTM [START_REF] Rozzi | Exact Coulomb cutoff technique for supercell calculations[END_REF][START_REF] Vico | Fast convolution with free-space Green's functions[END_REF] for Coulomb potential, then discuss its extensions to other convolution potentials.

The KTM, also named as "supercell" method [START_REF] Rozzi | Exact Coulomb cutoff technique for supercell calculations[END_REF], is commonly used in the physics community. The basic idea of KTM to compute nonlocal convolution potential on a finitesize domain is to screen unnecessary interaction U(r) at large distance and to compute potential generated by screened interaction which coincides with the originals. To be exact, in order to compute Φ over D, we cut off the interaction outside a larger domain that is usually chosen as a ball B G := {x||x| < G}, by simply setting U(x) =0, ∀x∈R d \B G , then apply trapezoidal quadrature to the modified Fourier integral whose integrand is not singular any more. Very recently, it was analyzed by Vico et al. [START_REF] Vico | Fast convolution with free-space Green's functions[END_REF] and extended to many important mathematical physics kernels including Coulomb, Helmholtz, biharmonic and Laplace kernels.

To be precise, we have

Φ(x) = R d U(x-y)ρ(y)dy = D U(x-y)ρ(y)dy = B G U(y)ρ(x-y)dy (2.1) = R d U G (y)ρ(x-y)dy = 1 (2π) d R d U G (k) ρ(k)e ik•x dk, x ∈ D, (2.2)
where D is the computational domain,

U G (x) = U(x) χ B G (x), χ B G (x)
is the characteristic function on ball B G with radius G not less than the diameter of D.

The Fourier transform of the truncated 3D Coulomb kernel

U G (x) = 1 4π|x| χ B G (x) is given exactly as follows U G (k) = R 3 U G (x)e -ik•x dx = k -2 (1-cos(kG)), k = |k|.
(2.3)

One can see clearly that U G (k) is not singular at the origin and

U G (0) = lim k→0 U G (k) = G 2 2 → ∞ as G → ∞. For the 2D Coulomb kernel, i.e., U(x) = 1 2π 1 |x| , the Fourier transform of U G reads as follows [58] U G (k) = 2π G 0 J 0 (kr)U(r)rdr = G 0 J 0 (kr)dr (2.4) = G 2 π J 1 (kG)SH 0 (kG)+ J 0 (kG) (2-πSH 1 (kG)) ,
where J 0 , J 1 are Bessel functions of the first-kind with index 0 and 1, and SH 0 ,SH 1 are Struve functions of order 0 and 1, respectively [START_REF] Zhang | Handbook of Mathematical Functions[END_REF]. Similarly, there is no singularity in

U G and U G (0) = lim k→0 U G (k) = G.
The optimal radius G is chosen as the diameter of D. For a rectangular domain D

= [-L,L] d , we have G = 2 √ dL, d = 2,3
, which suggests that a threefold and fourfold zeropadding of the density in each spatial direction is sufficient to guarantee spectral accuracy for the 2D and 3D case, respectively. However, direct fourfold zero-padding in 3D brings heavy burden on both memory and CPU time. In fact, with a precomputation, a twofold instead of threefold/fourfold zero-padding is sufficient [START_REF] Vico | Fast convolution with free-space Green's functions[END_REF]. The overall computation boils down to a 2D/3D discrete convolution. Take the 3D case as an example, we have

Φ ijk = ∑ i 1 ,j 1 ,k 1 ∈Λ ρ i 1 j 1 k 1 T i-i 1 ,j-j 1 ,k-k 1 , (i,j,k) ∈ Λ, (2.5)
where the convolution tensor T i,j,k is computed as inverse discrete Fourier transform of U G and is the same size as the discrete density ρ due to radial symmetry property of U G , and the index set is

Λ := {(i,j,k) i = 1,..., N x , j = 1,..., N y , k = 1,..., N z }, where N x , N y , N z
is the number of grid points in x-,yand z-direction, respectively. The multidimensional discrete convolution can be accelerated with discrete Fast Fourier Transform (FFT) [START_REF] Russell | Multidimensional Digital Signal Processing[END_REF] of size 2N tot (N tot := N x N y N z ), thus making KTM optimal in terms of efficiency. The acceleration algorithm is purely algebraic, which allows for an automatic implementation as long as T is available.

It is worthwhile to point out that KTM applies readily once the Fourier transform of the truncated kernel is available analytically or numerically. For example, the effective 2D convolution kernel that is derived from dimension reduction of Schr ödinger system [START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF][START_REF] Bao | A simple and efficient numerical method for computing dynamics of rotating dipolar Bose-Einstein condensation via a rotating Lagrange coordinate[END_REF] reads as follows

U(r) := 1 (2π) 3/2 R e -s 2 /2 √ r 2 +s 2 ε 2 ds = 1 ε(2π) 3/2 e r 2 4ε 2 K 0 ( r 2 4ε 2 ), 0 < ε < 1, (2.6)
where K 0 (x) is the modified Bessel function of the second-kind. The Fourier transform of its truncated kernel:

U G (k) = 2π G 0 U(r)J 0 (kr)rdr, k ∈ R, (2.7)
can be evaluated via adaptive Gauss-Kronrod quadrature accurately. The above method is readily applied to nonlocal potential generated by such kernel. The ground state of quasi-2D dipolar BEC is to be studied later in another paper.

The accuracy of KTM when applied to 2D/3D DDI and 2D Coulomb potential is confirmed numerically by the following two examples. Here, we consider the DDI with a general dipole interaction kernel U dip (x) with two different dipole orientations n and m as [START_REF] Baranov | Theoretical progress in many body physics of dipolar gases[END_REF][START_REF] Jiang | Fast and accurate evaluation of nonlocal Coulomb and dipole-dipole interactions via the nonuniform FFT[END_REF][START_REF] Kawaguchi | Spinor Bose-Einstein condensates[END_REF][START_REF] O'dell | Exact hydrodynamics of a trapped dipolar Bose-Einstein condensate[END_REF]]

U dip (x) =    -(n•m)δ(x)-3∂ nm 1 4π|x| , x ∈ R 3 , -3 2 ∂ n ⊥ m ⊥ -n 3 m 3 ∇ 2 ⊥ 1 2π|x| , x ∈ R 2 .
(2.8)

Here n = (n 1 ,n 2 ,n 3 ) T and m = (m 1 ,m 2 ,m 3 ) T are two given unit vectors representing the two dipole orientations,

m ⊥ = (m 1 ,m 2 ) T , ∂ m ⊥ = m ⊥ •∇ ⊥ and ∂ n ⊥ m ⊥ = ∂ n ⊥ (∂ m ⊥ ).
In this section, unless stated otherwise, we choose the mesh size as h x = h y = h z = h and denote the mesh grid as M h . The numerical accuracy is measured in relative maximum norm over M h as follows

e h := Φ h -Φ exact l ∞ Φ exact l ∞ = max x∈M h |Φ h -Φ exact | max x∈M h |Φ exact | , (2.9)
where Φ h is the numerical solution obtained with mesh size h and Φ exact is the exact solution.

Example 2.1. Dipole-Dipole Interaction in 3D.

Choose the source density ρ(x) = e -|x| 2 /σ 2 , with σ > 0. The 3D DDI with two dipole orientations n and m in (2.8) can be given explicitly as

Φ(x) = -(n•m)ρ(x)-3 ∂ nm σ 2 √ π 4 Erf(r/σ) r/σ = -(n•m)ρ(x)-3 n T B(x)m, (2.10)
where the matrix B(x) = (b jl (x)) 3 j,l=1 is given as

b jl (x) = σ 2 2r 2 e -r 2 σ 2 - σ 3 √ π 4r 3 Erf r σ δ jl +x j x l - 3 σ 2 2 r 4 e -r 2 σ 2 - 1 r 2 e -r 2 σ 2 + 3 σ 3 √ π 4 r 5 Erf r σ ,
(2.11) with δ jl the Kronecker delta, x = (x 1 ,x 2 ,x 3 ) T and Erf(r) = 2 √ π r 0 e -s 2 ds the error function. We choose σ = 1.4 and compute Φ(x) on a uniform mesh grid of domain D := [-L,L] 3 with mesh size h. Table 1 shows the relative errors of 3D DDI with different dipole axis, i.e., n = (0.82778,0.41505,-0.37751) T and m = (0.31180,0.93780,-0.15214) T . Choose the source density as ρ(x) = e -|x| 2 /σ 2 , with σ > 0. The exact 2D Coulomb potential is

Φ(x) = √ π σ 2 I 0 |x| 2 2σ 2 e -|x| 2 2σ 2 , (2.12)
and the exact DDI with two dipole orientations n ⊥ and m ⊥ in (2.8) is given below [START_REF] Jiang | Fast and accurate evaluation of nonlocal Coulomb and dipole-dipole interactions via the nonuniform FFT[END_REF]:

Φ(x) = 3 √ π e -r 4σ (n ⊥ •m ⊥ )(I 0 (r)-I 1 (r))- 2(x•n ⊥ )(x•m ⊥ ) σ 2 I 0 (r) - 1+2r 2r I 1 (r) + 3 √ π n 3 m 3 re -r σ I 0 (r)-I 1 (r)- I 0 (r) 2r , (2.13)
where r = |x| 2 2σ 2 , I 0 and I 1 are the modified Bessel functions of order 0 and 1, respectively [START_REF] Zhang | Handbook of Mathematical Functions[END_REF]. Here, we choose σ = √ 1.2 and dipole axis as n ⊥ = (0,-0.896) T , n 3 = 0.44404, m ⊥ = (0,-0.52476) T and m 3 = 0.85125. Table 2 shows the errors of 2D Coulomb potential and DDI obtained with different mesh size h on domain

D := [-L,L] 2 .
[-L,L] 2 . Coulomb h = 2 h = 1 h = 1/2 h = 1/4 h = 1/8 L =

Numerical method

In this section, we propose an efficient and accurate numerical method for computing the ground state by combining PCG method and KTM for DDI evaluation.

Discretization

To numerically solve the minimization problem, the wave function φ ∈ L 2 (R d ) shall be discretized. High spatial resolution schemes are required due to the presence of vortices in fast rotating system. High order finite difference/element schemes were proposed to compute the ground states of the rotating BEC without DDI [START_REF] Danaila | A finite element method with mesh adaptivity for computing vortex states in fast-rotating Bose-Einstein condensates[END_REF][START_REF] Danaila | A new Sobolev gradient method for direct minimization of the Gross-Pitaevskii energy with rotation[END_REF]. Here we apply the standard Fourier pseudo-spectral discretization [START_REF] Bao | Efficient numerical methods for computing ground states and dynamics of dipolar Bose-Einstein condensates[END_REF][START_REF] Bao | Efficient and spectrally accurate numerical methods for computing ground and first excited states in Bose-Einstein condensates[END_REF]. To this end, we first truncate the wave function to a bounded rectangular domain D with periodic B.C. To simplify the notation, we only present the 2D discetization, extension to higher dimension is straightforward. Let us set D = [L x ,R x ]×[L y ,R y ] and let M and N be two positive even integers. We fix

h x = R x -L x M and h y = R y -L y N
as the mesh sizes in the xand y-directions, respectively. Define the following sets of indices and grids points

O MN = (m, n) ∈ N 2 |0 ≤ m ≤ M, 0 ≤ n ≤ N , O MN = (p, q) ∈ N 2 | -M/2 ≤ p ≤ M/2-1, -N/2 ≤ q ≤ N/2-1 , Q xy = (x m , y n ) =: (L x +mh x ,L y +nh y ), (m,n) ∈ O MN .
Moreover, we introduce the following functions T pq (x,y) = e iν x p (x-L x ) e iν y q (y-L y ) , with

ν x p = 2π p R x -L x , ν y q = 2πq R y -L y (p,q) ∈ O MN . (3.1)
Define f mn =: f (x m , y n ) as the value of an abstract function f (x,y) at grid (x m , y n ) ∈ Q xy and f as the vector with components { f mn , (m,n) ∈ O MN }. The Fourier pseudo-spectral approximations of φ at Q xy and the operators ∆ & L z read as:

φ mn ≈ φ mn = M/2-1 ∑ p=-M/2 N/2-1 ∑ q=-N/2 ( φ) pq T pq (x m ,y n ), (m, n) ∈ O MN , (3.2) (∆φ) mn ≈ ( ∆ φ) mn = M/2-1 ∑ p=-M/2 N/2-1 ∑ q=-N/2 (ν x p ) 2 +(ν y q ) 2 ( φ) pq T pq (x m ,y n ), (3.3)
(L z φ) mn ≈ ( L z φ) mn = ω M/2-1 ∑ p=-M/2 N/2-1 ∑ q=-N/2 (x m ν y q -y n ν x p ) ( φ) pq T pq (x m ,y n ), (3.4)
where ( φ) pq is the discrete Fourier transform coefficient of the vector φ. By further defining the vector V , | φ| 2 and Φ as

( V ) mn = V mn , ( | φ| 2 ) mn = | φ mn | 2 , ( Φ ) mn = Φ mn , (3.5)
the Hamiltonian operator H φ and the total energy are discretized as following

H φ ≈ H φ = - 1 2 ∆ + V + β | φ| 2 +λ Φ -ω L z , (3.6)
E (φ) ≈ E ( φ) = φ, H φ φ - 1 2 β | φ| 2 2 - 1 2 λ Φ , | φ| 2 . (3.7)
Here, •, • and • are the standard C MN inner product and l 2 norm, respectively. Moreover, Φ is the approximation of Φ with density ρ =|φ| 2 replaced by ρ =| φ| 2 and is numerically evaluated by KTM. Hereafter, to simplify the presentation, we omit the "∼" and drop the brackets in the operators. Moreover, we denote f n as the value of f at the n-th (n ≥ 0) step in an iterative algorithm.

The preconditioned conjugate gradient (PCG) method

The classical approaches coming from PDE theory to compute the ground state of a BEC are mainly based on directly solving the corresponding imaginary-time equation (ITE) [START_REF] Chiofalo | Ground state of trapped interacting Bose-Einstein condensates by an explicit imaginary-time algorithm[END_REF]. As such, these methods lack the insight of minimization algorithms. The spherical constraint makes minimization problem (1.11) non-convex and hinders the application of classical minimization algorithms. Nevertheless, the general theory of optimization algorithm on Riemannian manifolds has been developed in [START_REF] Absil | Optimization Algorithms on Matrix Manifolds[END_REF][START_REF] Edelman | The geometry of algorithms with orthogonality constraints[END_REF]. Recently, following the approaches in [START_REF] Absil | Optimization Algorithms on Matrix Manifolds[END_REF][START_REF] Edelman | The geometry of algorithms with orthogonality constraints[END_REF], Antoine et al. proposed a (constrained analogues) preconditioned conjugate gradient algorithm on the manifold S to compute the ground states of rotating BEC without DDI [START_REF] Antoine | Efficient spectral computation of the stationary states of rotating Bose-Einstein condensates by preconditioned nonlinear conjugate gradient methods[END_REF]. Here, we extend this approach by integrating KTM (in section 2) to compute the ground states of rotating dipolar BEC. Assume P is an arbitrary symmetric definite preconditioner which will be determined later, the preconditioned conjugate gradient method to solve the minimization problem (1.11) is listed in Algorithm 1.

Algorithm 1

The preconditioned conjugate gradient method n = 0, given initial data φ n while not converged do

µ n = H φ n φ n ,φ n r n = H φ n φ n -µ n φ n β cg n = Re r n -r n-1 ,Pr n / r n-1 ,Pr n-1 β cg n = max(β cg n ,0) d n = -Pr n + β cg n p n-1 p n = d n -Re d n ,φ n φ n θ n = argmin θ E (cos(θ)φ n +sin(θ)p n / p n ) φ n+1 = cos(θ n )φ n +sin(θ n )p n / p n n = n+1 end while
Here, Re f denotes the real part of f . To find θ n , one can generally perform the linesearch

θ n = argmin θ E (cos(θ)φ n +sin(θ)p n / p n ), (3.8)
which is a one-dimensional minimization problem. Since E (θ) is not a quadratic function, this nonlinear minimization problem usually requires a number of evaluations. However, one can explicitly write out the formula of E (θ) in terms of sin(θ) and cos(θ). All the resulted coefficients that require FFTs could be pre-computed. Since the FFT step dominates the computation cost when computing the energy, the evaluation of E (θ) at many points does not increase the cost significantly than the evaluation at a single point. Hence, it is feasible to use a standard one-dimensional minimization routine to solve (3.8).

Alternatively, we can also obtain a simple and cheaper approximation of θ n . To this end, expanding E (θ n ) up to second order in θ n , we have

E (θ n ) =: E (φ n+1 ) = E (cos(θ n )φ n +γ n sin(θ n )p n ) ≈ E (φ n )+2γ n θ n Re p n , H φ n φ n -θ 2 n µ n -γ 2 n p n , H φ n p n -γ 2 n Re g n , p n . (3.9)
Here, γ n = 1/ p n , and g n = 2(βρ n pφ +λU dip * ρ n pφ ) φ n with ρ n pφ = Re(φ n pn ). First, we check if the coefficient of the first-order term is positive. If so, we set β cg n = 0, i.e., using the preconditioned steepest descent (PSD) method [START_REF] Antoine | Efficient spectral computation of the stationary states of rotating Bose-Einstein condensates by preconditioned nonlinear conjugate gradient methods[END_REF], and choose a sufficiently small value for θ n . With this choice, noticing that P is positive definite and φ n ,r n = 0, one can easily obtain that

2γ n Re p n , H φ n φ n = 2 Re p n ,r n +µ n φ n p n = 2 Re d n -Re d n ,φ n φ n ,r n +µ n φ n d n -Re d n ,φ n φ n = - 2 Re r n ,Pr n d n -Re d n ,φ n φ n < 0. (3.10)
Hence, the energy decreases for a sufficiently small θ n . Otherwise, we check further if the second-order term is positive, i.e.,

µ n < γ 2 n p n , H φ n p n +Re ρ n , p n . ( 3.11)
If so, minimizing (3.9) with respect to θ n yields

θ opt n = γ n Re p n , H φ n φ n µ n -γ 2 n p n , H φ n p n +Re g n , p n . ( 3.12)
If (3.11) is not fulfilled, again we set β cg n = 0 and choose a sufficiently small θ n . This approach guarantees the energy decay and hence the convergence of the algorithm.

In our numerical test, we found that this precautions of checking and using a step size control mechanism are useful in the first stage of the algorithm to locate the neighborhood of a minimum. Once a minimum is approximately located, both (3.10) & (3.11) are always fulfilled and the step size choice (3.12) always decreases the energy.

The discretization of the minimization problem (1.11) on a fixed fine mesh usually leads to a problem of huge size. Besides, starting for any commonly used initial data (Gaussian-type functions with phase imprinting , Thomas-Fermi approximation or their combinations), which are not a good initial guess especially for fast rotating systems, it usually takes relatively long time for the PCG scheme to converge. These together make the computation expensive, especially for high dimensional case. A useful technique is to apply the cascadic multi-grid method [START_REF] Bornemann | The Cascadic multigrid method for elliptic problems[END_REF]. The approach is as follows: starting from any commonly used initial data, employ the PCG algorithm to compute the ground state on a coarsest mesh with grid number N p × N p (N p = 2 p ). Denote the corresponding ground state as φ A ground state at coarser level usually will provide a much better initial data for the refined level than the commonly-used initial guess, hence a faster convergence is expected. This is actually confirmed by the numerical test in Example 4.3, which shows numerically that cascadic multigrid approach is more likely to converge faster to a ground state with lower energy than the fix-grid approach.

Stopping criteria and preconditioner

As discussed in [START_REF] Antoine | Efficient spectral computation of the stationary states of rotating Bose-Einstein condensates by preconditioned nonlinear conjugate gradient methods[END_REF], there are three types of stopping criteria one can use. A commonly used one is based on

φ ∞ err := φ n+1 -φ n ∞ ≤ ε. (3.13)
The other possible two are based on the norm of the symmetric-covariant residual

r n,∞ err := r n ∞ ≤ ε, (3.14)
or the symmetry-invariant energy difference

E n,∞ err := |E (φ n+1 )-E (φ n )| ≤ ε. (3.15)
Stopping criteria (3.13) are commonly used in the numerical schemes based on solving the ITE [START_REF] Chiofalo | Ground state of trapped interacting Bose-Einstein condensates by an explicit imaginary-time algorithm[END_REF]. However, this could be problematic if the minima are not isolated but form a continuum due to symmetries. For instance, when the dipole orientation is parallel to the rotating axis, the energy has rotational invariance. The third one (3.15) converges faster than the other two: assume φ * is a minimum, then

E (φ n )-E (φ * ) = O( φ n -φ * 2 ).
Extensive numerical results show that a criteria based on (3.13) or (3.14) would lead to long computational times, especially for large rotation ω, even without changing the total energy. Hence, in this paper, we always use the energy based stopping criteria (3.15).

The preconditioner P, which acts on the descent direction r n , is to make r n point closer to the minimum. It should be an approximation of the inverse of the Hessian matrix of the problem. If no preconditioner is applied, we denote P as

P I = I, (3. 16
)
where I is the identity matrix. Similar to [START_REF] Antoine | Efficient spectral computation of the stationary states of rotating Bose-Einstein condensates by preconditioned nonlinear conjugate gradient methods[END_REF], there are three other types of preconditioner.

Here we simply list their formula and properties. For a detailed analysis, we refer to [START_REF] Antoine | Efficient spectral computation of the stationary states of rotating Bose-Einstein condensates by preconditioned nonlinear conjugate gradient methods[END_REF].

Kinetic energy preconditioner.

The first one is to use only the kinetic energy term

P = P ∆ = (α ∆ -∆/2) -1 , (3.17)
where α ∆ is a positive shifting constant to construct an invertible operator. Here we choose as

α ∆ = μn := - 1 2 φ n , ∆φ n + V , |φ n | 2 + β |φ n | 2 , |φ n | 2 + λ Φ n ,|φ n | 2 . (3.18)
On a fixed computational domain, the effect of this preconditioner is to make the number of iterations independent of the mesh size h. However, this operator is not bounded in the whole space R d . In fact, the number of iterations increases as L increases. In addition, when β or λ is large, the nonlinear terms become dominant. The kinetic energy preconditioner becomes inefficient if any of the parameters β, λ or L is large.

Poten-Int energy preconditioner. The second approach is to use the potential and interaction energy terms for the preconditioner:

P = P V = α V +V + β|φ n | 2 +|λ|(1+sign(Φ n ))Φ n /2 -1 . ( 3.19)
Similarly, α V is a positive shifting constant to construct an invertible operator. We found that it is efficient to use α V = μn , and we always use this parameter for the corresponding numerical tests. Dual to the kinetic energy preconditioner, in terms of iteration number, the performance of this preconditioner deteriorates as the spatial resolution is increased (i.e., mesh size h is decreased). However, it is stable when L, β or λ is increased.

Combined preconditioner.

A third approach is to combine the two previous ones:

P = P C = P 1/2 V P ∆ P 1/2 V . (3.20)
This preconditioner has a stable performance independent of both the domain size L and/or the spatial resolution (mesh size h).

The performance of these preconditioners are confirmed in the numerical Example 4.1. The Combined preconditioner outperforms the other two, especially for small mesh sizes h and large domain sizes L, which is indeed necessary for computing ground states of the BECs with a fast rotation ω and large nonlinearity β. Therefore, in the practical computations, we always choose the Combined preconditioner.

Numerical results

Firstly, we denote the method by PCG ν (ν = I,∆,V) when the preconditioner P ν (3.16)-(3.20) is integrated in the PCG algorithm (Algorithm 1). Similarly, we denote by Multi-PCG ν if the cascadic multigrid algorithm is applied (Algorithm 2). We also define #iter as the number of iterations for PCG ν to converge, and shall always adopt the stopping criteria (3.15).

In the following numerical test, unless stated otherwise, we fix ε = 10 -14 and let

V(x) = |x| 2 /2, x ∈ D =: [-L,L] d . ( 4.1)
For the initial data, there are 10 types of frequently used initial data:

(a) φ a (x) = 1 √ π e -|x| 2 2 , (b) φ b (x) = (x+iy)φ a (x), ( b) φ b (x) = φb (x), (4.2)
(c) φ c (x) = φ a (x)+φ b (x) φ a (x)+φ b (x) , (c) φ c(x) = φc (x), (4.3)
(d) φ d (x) = (1-ω)φ a (x)+ωφ b (x) (1-ω)φ a (x)+ωφ b (x) , (d) φ d(x) = φd (x), (4.4
)

(e) φ e (x) = ωφ a (x)+(1-ω)φ b (x) ωφ a (x)+(1-ω)φ b (x) , (e) φ ē(x) = φe (x), ( f ) φ TF g (x) φ TF g (x) , (4.5)
where

φ TF g = (µ TF g -V(x))/β, V(x) < µ TF g , 0,
otherwise, Here, we compare the performance of PCG ν (ν = I,∆,V) with respect to different domain sizes L and mesh sizes h x = h y = h = 2 -p . To this end, we take d = 2, ω = 0, β = 250, λ = 100 and n = (1,0,0) T . Fig. 1 shows the total iteration number of these schemes with different h and L. From this figure and additional numerical results not shown here for brevity, we can see that: (i). Preconditioned algorithms outperform its unpreconditioned counterparts. (ii) The Poten-Int energy preconditioner P V enables the solver to stabilise the performance as L increases for a fixed h. On the contrary, the kinetic energy preconditioner P ∆ prevents the deterioration as h decreases for a fixed L. However, the performance of PCG V (PCG ∆ ) worsens as h decreases for a fixed L (as L increases for a fixed h). These observations confirm our expectation in subsection 3.3. (iii) The Combined preconditioner P C makes the solvers performance almost independent of both L and h, and outperforms all the others. In the following examples, we always use P C . Case 3. Fix the γ = λ β , vary both λ and β. Here, we consider γ = 0.25, 0.5 & 0.95. Fig. 2 illustrates the number of iterations of PCG C for Cases 1-3. From this figure and other numerical experiments now shown here, we can see that: (i) When λ increases (cf. Fig. 2 (a) &(c)), the iteration number also increases. Although the increment is small, it still implies a weak dependence on λ. Hence, the Poten-Int energy preconditioner P V should be analyzed and investigated further to remedy this dependence. (ii) For a fixed λ (cf. Fig. 2 (b)), the iteration number for PCG C actually oscillates in a small regime, which shows the very slight dependence on β. Nevertheless, the dependence on β (also on λ) would generally become stronger for a rotating system, i.e., ω = 0. It would be quite interesting trying to construct a robust preconditioner that also includes the rotating effects to get a weaker ω-dependence. We leave it as future work. Tables 3 and4 list respectively the energies obtained by PCG C and Multi-PCG C for different initial data and ω. The underline ones are the lowest energies of the converged states φ g (x) and the corresponding CPU times are listed in the same table. Meanwhile, the contour plots of the corresponding density |φ g (x)| 2 are shown in Fig. 3. From these tables and figures, we can see that: (i) Usually, algorithms with initial data of type (d) or (d) converges to the stationary states with the lowest energy. (ii) The cascadic multi-grid approach is more robust than the fix-grid one. It outperforms the fix-grid approach in terms of CPU time and provides the possibility to obtain a stationary state with lower energy. 4.

µ TF g = 1 2 (4βγ x γ y /π) 1/2 , d = 2, (15βγ x γ y γ z /4π) 2/5 , d = 3. ( 4
respectively. We set n = (cosθ,sinθ,0) T , β = (250-λ) √ 5/π and consider two cases: (i) Case 1: let λ = 175, ω = 0.7 and/or 0.9, vary θ. (ii) Case 2: let λ = 250, ω = 0.7 and/or 0.9, vary θ. Fig. 4 shows the contour plots of the ground states |φ g | 2 for Cases 1 & 2. From this figure, we can see that the dipole-dipole interaction affects the ground states essentially. The Abrikosov vortex lattice for rotating non-dipolar BEC is no longer valid in BEC with strong DDI, instead the rectangular vortex lattice appear. The vortex lattice orientation are parallel to the dipole orientation. Indeed, the presence of DDI brings about a much more rich phase diagram for ground states. However, to our best knowledge, a thorough study of the structure of the vortex lattice is still lacking to date. This might be tackled in future with the help of our efficient numerical scheme.

Example 4.5. Here, we apply the Multi-PCG C to solve a more difficult problem. We compute the ground state of BEC with very large nonlinearity and fast rotation. To this end, we let ω = 3.5, β = 5000, λ = 3000 and n = (cosθ,sinθ,0) T . For comparison, we also compute the one with λ =0 (other parameters are kept unchanged). The trapping potential is chosen as the harmonic plus quartic potential [START_REF] Antoine | Efficient spectral computation of the stationary states of rotating Bose-Einstein condensates by preconditioned nonlinear conjugate gradient methods[END_REF][START_REF] Danaila | A finite element method with mesh adaptivity for computing vortex states in fast-rotating Bose-Einstein condensates[END_REF][START_REF] Zeng | Efficiently computing vortex lattices in rapid rotating Bose-Einstein condensates[END_REF] V(x) = (3|x| 4 -4|x| 2 )/40. Example 4.6. Here, we compute the ground states in 3D. To this end, we take L = 8 and apply the Multi-PCG C with coarsest grids number N =32 and finest one N =128. We take β = 401.432, and consider the following six cases: (i) ω = 0.7, λ = 0.9β and n = (0,0,1) T ; (ii) ω = 0.7, λ = 0.9β and n = (1,0,0) T ; (iii) ω = 0.8, λ = 0.9β and n = (0,0,1) T ; (iv) ω = 0.8, λ = 0.9β and n = (1,0,0) T ; (v) ω = 0.9, λ = 0.5β and n = (0,0,1) T ; (vi) ω = 0.9, λ = 0.5β and n = (1,0,0) T .

Conclusion

We proposed an original preconditioned conjugated gradient (PCG) algorithm integrated with a newly developed kernel truncation method (KTM) for nonlocal potential evaluation to compute the ground states of rotating dipolar BECs. The KTM is numerically proved to be spectrally accurate and efficient. The resulted PCG-KTM algorithm, which is simple to implement, turns out to be robust and accurate. It is very efficient to compute the ground states of GPE with high nonlinearities in both the contact interaction and DDI, even for realistic 3D problems with fast rotation. Our numerical methods provide a useful tool to investigate the complicated phase diagram of the ground states of various fast rotating dipolar systems.

pg. 2 Algorithm 2

 22 Next refine the mesh by half, i.e., with N p+1 ×N p+1 mesh grids. Interpolate φ p g on the refined grid points to get a new initial data and call PCG again to obtain φ p+1 g , keep repeating this process until the finest grid N×N is reached. The approach is summarized in Algorithm The cascadic multigrid PCG method given an initial grid number N = 2 p and initial mesh Q p xy given an initial data φ p 0 on Q p xy while not reach the finest mesh do launch Algorithm 1 to compute the GS at level p: φ

. 6 )

 6 In the following examples, except Example 4.3, we choose the initial data as type (d) if β < 500 and type ( f ) otherwise. The algorithms were implemented in Matlab (Release 8.1.0), and run on a 2.27GH Intel(R) Xeon(R) CPU E5520 with a 8 MB cache in Debian GNU/Linux.Example 4.1.

Figure 1 :Example 4 . 2 .

 142 Figure 1: Example 4.1. Number of iterations for PCG ν (ν = I,∆,V) vs. mesh size h.

Figure 2 : 3 . 4 . 3 .

 2343 Figure 2: Example 4.2. Number of iterations to converge in PCG C for Cases 1-3.

Figure 3 :

 3 Figure 3: Contour plots of |φ g (x)| 2 corresponding to the lowest energies in Table4.

(4. 7 )

 7 The domain size, coarsest and finest grids number are chosen as L = 12, N = 64 and N = 512, respectively. Fig. 5 depicts the contour plot of the density function |φ g (x)| 2 for different λ & θ. From this figure, we can see that: (i) When DDI is absent, the structure of the density profile

Figure 4 :

 4 Figure 4: Example 4.4. Contour plot of the density |φ g | 2 of the ground states of the 2D rotating dipolar BEC with different ω in Case 1 (upper two rows) and Case 2 (lower two rows) for θ = 0, π 6 , π 2 ,-π 6 (left to right).

Figure 5 :

 5 Figure 5: Example 4.5. Contour plot of the density |φ g | 2 of the ground states of the 2D rotating BEC with ω = 3.5, β = 5000 and diffierent λ & n = (cosθ,sinθ,0) T . The CPU cost for these three cases are respectively: 9478s (left), 77500s (middle) and 58795s (right).

Fig. 6

 6 Fig.6shows the corresponding isosurface |φ g (x)| 2 = 4E-4 and the surface plot of |φ g (x,y,z = 0)| 2 for these six cases. The CPU times for these six cases are respectively 5418s, 6845s, 5073s, 4404s, 23408s, 5289s. From the figure, we can see that the DDI modifies the density profile of the ground state significantly. The density profile are elongated along the dipole orientation n (cf. Fig.6(i) & (iii)). Besides, the vortex lines are always parallel to the rotation axis, with the vortex lattice orientated also along the dipole axis n (cf. Fig.6(iii) & (iv)).

Figure 6 :

 6 Figure 6: Example 4.6. Isosurface |φ g (x)| 2 =4e-4 ((i) & (iii)) and surface plot of |φ g (x,y,z= 0)| 2 ((ii) & (iv)) for the six cases (top to bottom then left to right). The CPU cost for these six cases are respectively 5418s, 6845s, 5073s, 4404s, 23408s, 5289s.

Table 1 :

 1 Errors of 3D DDI with different dipole axis for different h and L.

		h = 2	h = 1	h = 1/2	h = 1/4	h = 1/8
	L = 6	9.5386E-01	5.3366E-02	1.3839E-06	3.2303E-09	2.4007E-09
	L = 8	9.5256E-01	5.6188E-02	1.1920E-06	4.6919E-16	6.8574E-16
	L =10	9.3653E-01	5.7520E-02	1.0087E-06	4.6919E-16	6.8574E-16
	Example 2.2. Dipole-Dipole Interaction and Coulomb potential in 2D.

Table 2 :

 2 Errors of 2D Coulomb potential and DDI with mesh size h on

Table 3 :

 3 Example 4.3. Fixed grid PCG C (with N = 2 8 ): converged energies and the CPU times (seconds) for the solution with lowest energy (which is underlined). Here, we apply the Multi-PCG C algorithm to compute the ground states of 2D rotating BEC with different n and ω. To this end, the domains are chosen as D = [-12,12] 2 and the coarsest and finest grids number are set to N = 64 & N = 256,

	ω	(a)	(b)	(b)	(c)	(c)	(d)	(d)	(e)	(e)	CPU
	0.6 5.7753 5.7762 5.7826 5.783 6.7723 5.783 5.7723 5.783 5.7723 1045
	0.7 5.3678 5.3448 5.3444 5.3415 5.3415 5.3415 5.3415 5.3415 5.3444 3270
	0.8 4.7282 4.7282 4.7285 4.7285 4.7285 4.7285 4.7282 4.72852 4.7282 1368
	0.9 3.773 3.773 3.773 3.773 3.773 3.773	3.773	3.773	3.773 1167
	Example 4.4.

Table 4 :

 4 Example 4.3. Cascadic multigrid Multi-PCG C (starting from the coarsest grid N = 26 to the finest grid N = 2 8 ): converged energies and the CPU times (seconds) for the solution with lowest energy (underlined).

Acknowledgements

Xavier Antoine thanks the support from the LIASFMA (University of Lorraine), ANR-12-MONU-0007-02 BECASIM and the hosting of the Beijing CSRC. Qinglin Tang thanks the support from the Singapore Ministry of Education Academic Research Fund Tier 2 R-146-000-223-112. Yong Zhang gratefully acknowledges the support of the Austrian Science Foundation FWF under grant No F41 (SFB "VICOM"), grant No F65 (SFB Com-