A numerical scheme for the ground state of rotating spin-1 Bose–Einstein condensates

We study the existence of nontrivial solution branches of three-coupled Gross–Pitaevskii equations (CGPEs), which are used as the mathematical model for rotating spin-1 Bose–Einstein condensates (BEC). The Lyapunov–Schmidt reduction is exploited to test the branching of nontrivial solution curves from the trivial one in some neighborhoods of bifurcation points. A multilevel continuation method is proposed for computing the ground state solution of rotating spin-1 BEC. By properly choosing the constraint conditions associated with the components of the parameter variable, the proposed algorithm can effectively compute the ground states of spin-1 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{87}Rb$$\end{document}87Rb and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{23}Na$$\end{document}23Na under rapid rotation. Extensive numerical results demonstrate the efficiency of the proposed algorithm. In particular, the affect of the magnetization on the CGPEs is investigated.

" * " over the components ϕ i of the wave functions denotes the complex conjugate. Notice that the interaction is either repulsive or attractive depending on the constant c 0 is positive or negative. Furthermore, the spin-exchange interaction can be antiferromagnetic or ferromagnetic depending on the constant c 2 is positive or negative. The GPE based on the mean-field theory is treated as the model for describing the physical phenomena of BEC at zero temperature. However, there are still several models which have been proposed to take into account finite temperatures effects in a quantum fluid. A well-known example is the Zaremba-Nikuni-Griffin (ZNG) method 15,16 , in which a dissipative GPE for the condensate wave function is coupled with a heat base described by the Boltzmann equation. Another simpler model is the so-called stochastic projected Gross-Pitaevskii equation (SPGPE), in which thermal fluctuations of the bosonic field are taken into account at stochastic forcing [17][18][19] .
For completeness we define the state of lowest energy of a BEC system with fixed number of particles as the ground state. The states with energies greater than the ground-state energy are called excited states. The linear Zeeman (LZ) energy and the quadratic Zeeman (QZ) energy are given by [20][21][22] . and respectively. Both the parameters p 0 and q 0 play important roles in the ground state phase diagram as well as the dynamics of spin-1 condensates. From Eqs. (2) and (3) we have Various numerical methods have been proposed for computing the ground state solution of both one-and two-component rotating BEC [23][24][25][26][27][28][29][30][31][32][33] . In particular, the preconditioned imaginary time evolution method (PITEM), or the so-called continuous normalized gradient flow (CNGF) was widely used. See 34,35 . Recently, the performance of PITEM and continuation methods on some test problems in boson-fermion mixtures was compared in 36 . Published articles on numerical study of spin-1 BEC is also abundant. See e.g., [37][38][39][40] . Research papers on numerical investigation of rotating spin-1 BEC can be found e.g., in 18,19,[41][42][43][44][45][46][47] .
In this paper, we investigate the existence of nontrivial solution curves of the CGPEs using the Lyapunov-Schmidt reduction. By performing a small perturbation of the cubic nonlinearity, we show how the nontrivial solution curves branching from bifurcation points on the trivial one. Our result is a 2D and three-component generalization of that in 59 . Next, we describe a novel multi-level continuation algorithm to compute the ground states of Eq. (10) for various values of the parameters, where the Fourier sine functions are used as the basis functions to discretize the CGPEs. In the first two levels of the algorithm we use the chemical potential µ , and then add the magnetic potential as the first and the second continuation parameters, respectively. In our numerical computations we consider the cases with magnetization M = 0 and M = 0 . For convenience we omit the quadratic Zeeman energy q in Eq. (10). Note that the numerical computations of the ground states for spin-1 BEC with quadratic Zeeman energy has been widely investigated in 40 . Instead of using M = 0 as the constraint condition in Eq. (7), we impose a more reflexible one in the second level, where M ∈ [0, 1] is fixed, and the L 2 -norm � · � is defined by Note that as we start to switch from the trivial solution curve to the nontrivial one of Eq. (10), the two-norm of the components ψ 1 and ψ −1 are relatively small compared to the value M = 0 . It is impossible that Eq. (7) will hold. Thus we multiply the magnetization M in Eq. (11) by 2 in the continuation process in order to keep Newton's method from divergence in the corrector step of the continuation algorithm. In the third level of the algorithm we intend to use the number of particles N as the third component of the parameter variable. Since the scales of µ, and the number of the particles N are quite different, we impose an artificial parameter ν as the third component of the parameter variable to control the increment of N . We will also apply the proposed algorithm to study how the wave function of Eq. (10) changes with respect to the angular velocity when ω > 1 , where we impose a harmonic plus quartic trap on the system. Note that the numerical computations for the ground states of fast rotating spin-1 BEC become difficult when the angular velocity ω > 1 , and it is getting more challenging as ω increases. To our knowledge, the physical phenomena of the ground states of rotating spin-1 BEC with M > 0 and rapidly rotating spin-1 BEC have never been reported in the literature.
The organization of this paper is as follows. In section "Existence of nontrivial solution curves" we present the existence of nontrivial solution curves branching from bifurcation points of the CGPEs. A multi-parameter continuation algorithm is proposed in section "A multilevel continuation algorithm" for computing the ground state of (rapidly) rotating spin-1 BEC. In section "Numerical results" we investigate numerically how the magnetization www.nature.com/scientificreports/ may affect the behavior of the CGPEs. Our numerical results demonstrate that various vortex lattices of 87 Rb and 23 Na can be observed. Finally, some concluding remarks are given in section "Conclusions".

Existence of nontrivial solution curves
In this section, we will show the existence of nontrivial solution branches of Eq. (10), where the wave functions near the bifurcation point satisfy �ψ 1 � 2 + �ψ 0 � 2 + �ψ −1 � 2 = O(ε) . We consider the scaling ψ l (x) = ε 1/2 φ l (x) , l = 1, 0, −1 , Then Eq. (10) becomes . In order to be consistent with the continuation algorithm we describe in Section 3, we set = 0 . The linear eigenvalue problem associated with Eq. (13) is given by For simplicity we let ω = 0 . The eigenpairs of Eq. (14) are as follows: where H k is the kth degree Hermite polynomial. The first few eigenfunctions are Note that the set of eigenfunctions u m,n (x, y) | m, n = 0, 1, 2, . . . forms an orthonormal basis for L 2 (R 2 ) under the inner product f , g = R 2 f (x)g(x) dx. Using the Lyapunov-Schmidt reduction it was shown in 58 that the bifurcations of a single NLS are pitchfork. For the case of BEC the coefficient of the cubic nonlinear term is positive. Thus the pitchfork bifurcations are supercritical where the solution curves turn to right. It is straightforward to prove that the bifurcations of Eq. (10) have the same properties mentioned above. The stability analysis for the CGPES was studied in 61 . It is expected that the stability analysis for spin-1 BEC can be treated in a similar way.
Various types of continuation algorithms have been proposed for computing the ground state and excited states of (rotating) BEC 60,[62][63][64] . In this section, we describe a multi-parameter continuation algorithm for computing the ground state solution of Eq. (10). It suffices to trace the solution curve branching from the minimum eigenvalue of the linearized Schrödinger equation (LSE) associated with Eq. (10). Starting with � � ≈ 0 near the trivial solution curve, we will follow this primary solution curve by the proposed continuation algorithm described below until the target point is reached, where the normalization � � 2 = 1 is satisfied. See Eq. (6). The target point we obtained is indeed the ground state solution of Eq. (10). A detailed comparison between the performance of the PITEM/CNGF and continuation methods was reported in 36 . We also refer to 60,61 for further discussions. Furthermore, we will obtain all solutions of rotating spin-1 BEC for any values of the particle number N (or the angular velocity ω ) on certain interval, say, N ∈ [ N 0 , N * ] for some positive number N * (or ω ∈ [ω 0 , ω * ] ). Note that in some cases the ground state solution of the NLS does not necessarily lie on the nontrivial solution curve branching from the minimum eigenvalue of the associated LSE. See e.g., 65,66 .
Theoretically we can use both the chemical potential µ and the parameter N as the two continuation parameters simultaneously. However, the continuation increment, namely, the step size for curve-tracing is relatively small, say, from 10 −1 to 10 −2 , depending on the curvature of the solution curve, compared to the scale of N . Therefore, it requires large number of continuation steps to trace the solution curve which can be very expensive.  To begin with, we consider the first level continuation algorithm with � = µ ∈ R as the continuation parameter, and set = ν = 0 . The discrete analogue of Eq. (23) is a nonlinear system of equations involving the parameter µ , and is given as where the tangent vector ẏ(s) is normalized so that and the Jacobian matrix DH(y(s)) ∈ R (6N 2 )×(6N 2 +1) is of full rank. It follows from Eq. (28) that the augmented Jacobian matrix is nonsingular for all s ∈ I except that at the primary bifurcation points on the trivial solution curve {(0, 0, 0, 0, 0, 0, µ) | µ ∈ R} , where the Jacobian matrix DH(y(s)) has rank deficiency. To switch from the trivial solution curve to the primary solution branch near the bifurcation point, we solve the perturbed nonlinear system for some perturbation vector d ∈ R 6N 2 . In general, the vector d in Eq. (29) is chosen so that it has the same mode as the eigenfunction of the associated linear eigenvalue problem. We refer to 68 and the further reference cited therein for details.
Right after we switch from the trivial solution curve to the primary solution curve, we perform the second level continuation algorithm by adding the magnetic potential as the second component of the parameter variable defined in Eq. (27). That is, we set � 2 : = (µ, ) ∈ R 2 , where the value of the wave function � � 2 ≈ k 0 is small enough for some positive constant k 0 . Note that if k 0 is too large, which means that we implement the first level continuation algorithm to trace the solution curve by neglecting the affect of the magnetic potential . The price is that the algorithm can not mimic the physical systems of rotating spin-1 BEC precisely. This would make the algorithm either diverge or fail to trace the solution curve we wish to follow. We refer to 39 for detailed discussions. Now we rewrite the magnetization (7) as which is added as the last equation to the nonlinear system of equations H(x, �) = 0 defined in Eq. (27). More precisely, we update H(x, �) = 0 by setting Now the second level continuation is exploited to trace the ground state solution curve of Eq. (31). We stop the implementation of the second level algorithm when the normalization condition � � 2 = 1 is satisfied.
Finally, we update the parameter variable � 2 = (µ, ) to a three-component variable � 3 = (µ, , ν) ∈ R 3 , where the last component of 3 is defined in Eq. (23). We express the normalization condition (6) as   www.nature.com/scientificreports/ Equation (32) will be added as the last equation to the nonlinear system of equations H (x, � 2 ) = 0 defined in Eq. (31). In other words, the nonlinear system H (x, � 2 ) = 0 will be updated to a new one, namely, Denote the Jacobian matrix of H by D H ∈ R (6N 2 +2)×(6N 2 +3) . We implement the third level algorithm to trace the ground state solution of Eq. (33). To compute the unit tangent vector

we solve the linear system
A new approximating point is predicted by the Euler predictor where δ (k) > 0 is the step length. Next, we correct the predicted point z (k+1,1) by performing Newton's iteration. We solve the linear system where z (k+1,j+1) = z (k+1,j) + w (j) , j = 1, 2, . . . . If the corrector increment w (j 0 ) and � H(z (k+1,j 0 +1) )� are small enough for some j 0 ∈ N , we obtain the next approximating point y (k+1) = z (k+1,j 0 +1) . We stop the curve-tracing when the target point is reached. Now the ground state solutions for various values of the coefficient N are available on the solution curve. Note that the state variables of Eqs. (27), (31), and (33) are the same. More precisely, the solution curves connect consecutively except that we would gain more information for the ground states as the number of components of the parameter variable increases.
The algorithm described above may be briefly summarized as follows.

Algorithm 3.2 A multi-level continuation method for rapidly rotating spin-1 BEC.
Initialization: k 0 := a given small positive number for implementing Level 2.
N 0 := initial particle number used in Levels 1 and 2.
Remark For tracing the ground state solution curve of ultrarapidly rotating spin-1 BEC, we use the angular velocity ω as the third components of the parameter variable by letting ω = ω 0 + ν σ in Eq. (24) for some constants ω 0 and σ . That is, ω ∈ [ω 0 , ω * ]. Table 1. Comparison of the parameters and the associated constraints between the algorithms for spin-1 BEC and rotating spin-1 BEC.

Numerical results
In this section we report the implementation results of Algorithm 3.2. The ratio of spin-independent and spindependent interactions is ∼ 0.48% 69 for 87 Rb , ∼ 1.5% 70 for 23 Na , and strongly ferromagnetic effect ∼ 45% 5,71 for 17 Li . For typical parameters the character length of trap a m ∼ 1µ m and s-wave scattering length ∼ 10nm. It has been satisfied for most experiments. By using a m as the length unit, for the ferromagnetic case 87 Rb we chose g n = 0.0885 , and g s = −0.00041 , and for the antiferromagnetic case 23 Na we chose g n = 0.0241 , and g s = 0.00075 . From 72,73 , we have −0.67 < Na 0 ( or a 2 )/a m < ∞ . The particle number of condensate is typically between 10 4 and 10 6 . In these cases, the GPE model is valid under the dilute condition. Specifically, in Example 4.1 we studied the convergence behavior of Algorithm 3.   16) . We traced the ground state solution curve of spin-1 BEC using 87 Rb and 23 Na for different number of basis functions until the normalization condition � � 2 = 1 was satisfied. Denote the corresponding chemical potentials and exact chemical potential by µ (N) , and by µ * , respectively, and the convergence rate and convergence order for 87 Rb or 23 Na by e −mN and N −Order , respectively 74 , which are given by and respectively. Tables 2-3 list the chemical potentials, the corresponding absolute errors, and the values of m and Order associated with the convergence behavior of 87 Rb and 23 Na , respectively. Figure 1 displays the convergence behavior of the chemical potentials for 87 Rb and 23 Na . The results given above show that the convergence rate of Algorithm 3.2 combined with the Fourier sine functions for the CGPEs, a nonlinear elliptic eigenvalue problem, is indeed exponential. Figure 2a,b show the graphs of the ground state solution of 87 Rb and 23 Na , respectively, with N = 1024.  1,19,46,47,75,76 . In addition, the two-norm of the third component ψ −1 almost equals to zero for all values of N because of the magnetization M = 0.5 . It is expected that if we increase the value of magnetization from M = 0.5 gradually, the two-norm of the third component ψ −1 will be zero, and the three-coupled GPEs reduce to rotating two-component BEC. Our result is similar to that of rotating two-component BEC shown in 23 . The result verifies the prediction numerically shown in Lemma 3.1.
That is, as the magnetization M increases, the governing Eq. (10) will change gradually from the twocoupled GPEs, and then to the single GPE.   Figure 5a shows the relationship between the chemical potential µ and the particle numbers N on the two-norm solution curve of the wave function obtained in implementing Levels 1-3. Figure 5b displays how the vortex lattice of the components evolve with respect to the particle number N . More precisely, when N ≈ 40,720 , two vortices of the components ψ 1 and ψ −1 start to be pinned together. When N ≈ 48,199 , vortices in the three component were pinned together to form a vortex-pair lattice, where the vortices of each pair had the same circulation. When the particle number is large enough, say, N ≥ 50,069 , the vortices of the three components exhibit a square lattice where two vortices remained to be pinned together. Similar phenomenon has been observed in published literature on rotating spin-1 BEC. See e.g., 19 . (ii) M = 0.3 : We set N 0 = 20,000 , and N * = 100,000 . Algorithm 3.2 was implemented to compute the ground state solution curve of the wave function . The result was depicted in Fig. 6a, which showed www.nature.com/scientificreports/ that the projection on the second component ψ 0 equals zero in Levels 1 and 2, and then increases slowly in Level 3. On the other hand, the projections on the components ψ 1 and ψ −1 increase in Levels 1 and 2, and decrease in Level 3, which separate from each other owing to the magnetization. Figure 6b displayed how the contour plots of the three components varied with respect to the value of N . When N ≥ 78,022 , we found that the vortices of the three components exhibited a square lattice, which are similar to the case in (i) with M = 0 . However, the two-norm of ψ −1 is smaller than the counterpart of (i).    www.nature.com/scientificreports/ which was also surrounded by a ring with yellow color. When ω = 3.30 , the contours of the three components were similar to those of the components when ω = 3.10.

Conclusions
We have applied the Lyapunov-Schmidt reduction to show the existence of nontrivial solution curves branching from eigenvalues of the linearized CGPEs. Based on the existence theory a multilevel pseudo-arclength continuation algorithm has been proposed which can efficiently compute the ground state solutions of rapidly rotating spin-1 BEC for both 87 Rb and 23 Na . Our numerical results have demonstrated that various types of vortex lattices could be obtained for both 87 Rb and 23 Na.
We remark the phenomenon exhibited in Example 4.2(ii). Owing to the repulsive interspecies interaction, it is intuitional that vacancies like vortices in one component are filled by droplets in another in order to lowering the energy of the system. However, the spinor degrees of freedom can provide a platform to study topological quantum phenomena in such multi-component system. This kind of BEC are called spinor BEC. For a spinor F = 1 BEC, the individual topological defects, half-quantum vortices in the polar phase, polar-core vortices, skyrmions (in 3D), and baby-skyrmions (in 2D) in the ferromagnetic phase have been discussed in 75,76 .
From Lemma 3.1 we may conclude that the magnetization M plays a key factor which makes Eq. (10) reduce to a single GPE when M = ±1 . The contours displayed in Figs. 4, 7, and 8 verify numerically when M ≥ 0.5 , Eq. (10) almost decays to a two-coupled GPEs. As we increase the magnetization from M = 0.5 gradually to M = 1 , it is expected that the two-coupled GPEs will decay to a single GPE. Finally, it would be of interest yet challenging to propose numerical methods for the ground state solution of the SPGPE for future studies.   www.nature.com/scientificreports/