{'fraction_non_alphanumeric': 0.05682749576525873, 'fraction_numerical': 0.02431561116878859, 'mean_word_length': 4.278915488895299, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Computing Ground States of Spin-1 Bose-Einstein Condensates by the Normalized Gradient Flow

In this paper, we propose an efficient and accurate numerical method for computing the ground state of spin-1 Bose-Einstein condensates (BEC) by using the normalized gradient flow or imaginary time method. The key idea is to find a third projection or normalization condition based on the relation between the chemical potentials so that the three projection parameters used in the projection step of the normalized gradient flow are uniquely determined by this condition as well as the other two physical conditions given by the conservation of total mass and total magnetization. This allows us to successfully extend the most popular and powerful normalized gradient flow or imaginary time method for computing the ground state of single component BEC to compute the ground state of spin-1 BEC. An efficient and accurate discretization scheme, the backward-forward Euler sine-pseudospectral method (BFSP), is proposed to discretize the normalized gradient flow. Extensive numerical results on ground states of spin-1 BEC with ferromagnetic/antiferromagnetic interaction and harmonic/optical lattice potential in one/three dimensions are reported to demonstrate the efficiency of our new numerical method.

1. Introduction. Research in low temperature dilute atomic quantum gases remains active for more than ten years after the experimental realizations of Bose-Einstein condensation (BEC) in alkali atomic gases in 1995 [2,12,17]. Extensive theoretical and experimental studies have been carried out to investigate various novel phenomena of the condensates. In earlier BEC experiments, the atoms were confined in magnetic trap [2,12,17], in which the spin degrees of freedom is frozen. The particles are described by a scalar model and the wavefunction of the particles is governed by the Gross-Pitaevskii equation (GPE) within the mean-field approximation [16,11,23]. In recent years, experimental achievement of spin-1 and spin-2 condensates [10,18,21,26,28] offers new regimes to study various quantum phenomena that are generally absent in a single component condensate. The spinor condensate is achieved experimentally when an optical trap, instead of a magnetic trap, is used to provide equal confinement for all hyperfine states.
The theoretical studies of spinor condensate have been carried out in several papers since the achievement of it in experiments [19,20,22,27]. In contrast to single component condensate, a spin-F (F ∈ N) condensate is described by a generalized coupled GPEs which consists of 2F + 1 equations, each governing one of the 2F + 1 hyperfine states (m F = −F, −F +1, ..., F −1, F ) within the mean-field approximation. For a spin-1 condensate, at temperature much lower than the critical temperature T c , the three-components wavefunction Ψ(x, t) = (ψ 1 (x, t), ψ 0 (x, t), ψ −1 (x, t) T are well described by the following coupled GPEs [27,29,30,31], Here, x = (x, y, z) T is the Cartesian coordinate vector, t is time, is the Planck constant, m is the atomic mass, and V (x) is the external trapping potential. When a harmonic trap potential is considered, with ω x , ω y and ω z being the trap frequencies in the x-, y-and z-direction, respectively.f and Re(f ) denote the conjugate and real part of the function f , respectively. There are two atomic collision terms, c 0 = 4π 2 3m (a 0 + 2a 2 ) and c 2 = 4π 2 3m (a 2 − a 0 ), expressed in terms of the s-wave scattering lengths, a 0 and a 2 , for scattering channel of total hyperfine spin 0 (anti-parallel spin collision) and spin 2 (parallel spin collision), respectively. The usual mean-field interaction, c 0 , is positive for repulsive interaction and negative for attractive interaction. The spin-exchange interaction, c 2 , is positive for antiferromagnetic interaction and negative for ferromagnetic interaction. The wave function is normalized according to where N is the total number of particles in the condensate. By introducing the dimensionless variables: t → t/ω m with ω m = min{ω x , ω y , ω z }, x → x a s with a s = mωm , ψ l → √ N ψ l /a 3/2 s (l = −1, 0, 1), we get the dimensionless coupled GPEs from (1.1)-(1.3) as [30,32,9]: and V (x) = 1 2 (γ 2 x x 2 +γ 2 y y 2 + γ 2 z z 2 ) with γ x = ωx ωm , γ y = ωy ωm and γ z = ωz ωm . Similar as those in single component BEC [24,8,3,6], in a disk-shaped condensation, i.e. ω x ≈ ω y and ω z ≫ ω x (⇐⇒ γ x = 1, γ y ≈ 1 and γ z ≫ 1 with ω m = ω x ), the 3D coupled GPEs (1.6)-(1.8) can be reduced to a 2D coupled GPEs; and in a cigar-shaped condensation, i.e. ω y ≫ ω x and ω z ≫ ω x (⇐⇒ γ x = 1, γ y ≫ 1 and γ z ≫ 1 with ω m = ω x ), the 3D coupled GPEs (1.6)-(1.8) can be reduced to a 1D coupled GPEs. Thus here we consider the dimensionless coupled GPEs in d-dimensions (d = 1, 2, 3): In the equations above, V (x) is a real-valued potential whose shape is determined by the type of system under investigation, β n ∝ N and β s ∝ N correspond to the dimensionless mean-field (spin-independent) and spin-exchange interaction, respectively. Three important invariants of (1.9)-(1.11) are the mass (or normalization) of the wave function (1.12) N (Ψ(·, t)) := Ψ(·, t) 2 := 0)) = M and the energy per particle A fundamental problem in studying BEC is to find the condensate stationary states Φ(x), in particular the ground state which is the lowest energy stationary state. In other words, the ground state, Φ g (x), is obtained from the minimization of the energy functional subject to the conservation of total mass and magnetization: Find (Φ g ∈ S) such that where the nonconvex set S is defined as (1.16) This is a nonconvex minimization problem. When β n ≥ 0 and β n ≥ |β s | and lim |x|→∞ V (x) = ∞, the existence of a minimizer of the nonconvex minimization problem (1.15) follows from the standard theory [25]. For understanding the uniqueness question note that E(α · Φ g ) = E(Φ g ) for all α = e iθ1 , e iθ0 , e iθ−1 T with θ 1 + θ −1 = 2θ 0 . Thus additional constraints have to be introduced to show the uniqueness. As derived in [9], by defining the Lagrangian we get the Euler-Lagrange equations associated to the minimization problem (1.15): Here µ and λ are the Lagrange multipliers (or chemical potentials) of the coupled GPEs (1.9)-(1.11). In addition, (1.18)-(1.20) is also a nonlinear eigenvalue problem with two constraints In fact, the nonlinear eigenvalue problem (1.18)-(1.20) can also be obtained from the coupled GPEs (1.9)-(1.11) by plugging ψ l (x, t) = e −iµ l t φ l (x) (l = 1, 0, −1) with Thus it is also called as time-independent coupled GPEs. In physics literatures, any eigenfucntion Φ of the nonlinear eigenvalue problem (1.18)-(1.20) under constraints (1.21) and (1.22), whose energy is larger than the energy of the ground state is called as an excited state of the coupled GPEs (1.9)-(1.11).
A widely used numerical method for computing the ground state of a single component condensate is the imaginary time method followed by an appropriate discretization scheme [15,5,3] to evolve the resulted gradient flow equation under normalization of the wavefunction, which is mathematically justified by using the normalized gradient flow [5,3]. However, it is not obvious that this most popular and powerful normalized gradient flow (or imaginary time method) could be directly extended to compute the ground state of spin-1 BEC. The reason is that we only have two normalization conditions (i.e. the two constraints: conservation of total mass and magnetization) which are insufficient to determine the three projection constants for the three components of the wavefunction used in the normalization step. In physics literatures, the imaginary time method is still applied to compute the ground state of spin-1 BEC through the introduction of a random variable to choose the three projection parameters in the projection step [30,32]. Of course, this is not a determinate and efficient way to compute the ground state of spin-1 BEC due to the choice of the random variable. Recently, Bao and Wang [9] have proposed a continuous normalized gradient flow (CNGF) for computing the ground state of spin-1 BEC. The CNGF is discretized by Crank-Nicolson finite difference method with a proper and very special way to deal with the nonlinear terms and thus the discretization scheme can be proved to be mass and magnetization conservative and energy diminishing in the discretized level [9]. However, at each time step, a fully nonlinear system must be solved which is a little tedious from computational point of view since the CNGF is an integral-differential equations (see details in (A.1)-(A.9)) which involves implicitly the Lagrange multipliers in the normalized gradient flow evolution [9]. The aim of this paper is to introduce a third normalization condition based on the relation between the chemical potentials of spin-1 BEC, in addition to the two existing normalization conditions given by the conservation of the total mass and magnetization. Thus we can completely determine the three projection constants used in the normalization step for the normalized gradient flow. This allows us to develop the most popular and powerful normalized gradient flow or imaginary time method to compute the ground state of spin-1 BEC.
The paper is organized as follows. In section 2, the normalized gradient flow is constructed by introducing the third projection or normalization condition for computing the ground state of spin-1 BEC. In section 3, the backward-forward Euler sine-pseudospectral method (BESP) is presented to discretize the normalized gradient flow. In section 4, ground states of spin-1 BEC are reported with ferromagnetic/antiferromagnetic interaction and harmonic/optical lattice potential in one/three dimensions, respectively. Finally, some conclusions are drawn in section 5.
2. The normalized gradient flow. In this section, we will construct the normalized gradient flow for computing the ground state of spin-1 BEC by introducing the third normalization condition.
Various algorithms for computing the minimizer of the nonconvex minimization problem (1.15) have been studied in literature. For instance, a continuous normalized gradient flow (CNGF) and its discretization that preserve the total mass and magnetization conservation and energy diminishing properties were presented in [9]. Perhaps one of the more popular and efficient techniques for dealing with the normalization constraints in (1.16) is through the following construction: choose a time step k = ∆t > 0 and denote time steps as t n = n k for n = 0, 1, 2, · · · . To adapt an efficient algorithm for the solution of the usual gradient flow to the minimization problem under constraints, it is natural to consider the following splitting (or projection) scheme, which was widely used in the physics literature for computing the ground state of BECs: followed by a projection step as and σ n l (l = −1, 0, 1) are projection constants and they are chosen such that In fact, the gradient flow (2.1)-(2.3) can be viewed as applying the steepest decent method to the energy functional E(Φ) in (1.14) without constraints, and (2.4)-(2.6) project the solution back to the unit sphere S in order to satisfy the constraints in (1.16). In addition, (2.1)-(2.3) can also be obtained from the coupled GPEs (1.9)-(1.11) by the change of variable t → −i t, that is why the algorithm is usually called as the imaginary time method in the physics literatures [15,5,3]. Plugging (2.4)-(2.6) into (2.7), we obtain There are three unknowns and only two equations in the above nonlinear system, so the solution is undetermined! In order to determine the projection constants σ n l (l = −1, 0, 1), we need to find an additional equation. Based on the relation between the chemical potentials in (1.23) and the continuous normalized gradient flow proposed in [9] for computing the ground state of spin-1 BEC, see details in Appendix A, we propose to use the following equation as the third normalization condition (2.10) σ n 1 σ n −1 = (σ n 0 ) 2 . Solving the nonlinear system (2.8), (2.9) and (2.10), see details in Appendix B, we get explicitly the projection constants as (2.11) .
From the numerical point of view, the gradient flow (2.1)-(2.3) can be solved via traditional techniques, and the normalization of the gradient flow is simply achieved by a projection at the end of each time step.
Due to the trapping potential V (x) given by (1.4), the solution Φ(x, t) decays to zero exponentially fast when |x| → ∞. Thus in practical computation, we truncate the problem into a bounded computational domain Ω x (chosen as an interval (a, b) in 1D, a rectangle (a, b) × (c, d) in 2D, and a box (a, b) × (c, d) × (e, f ) in 3D, with |a|, |c|, |e|, b, d and f sufficiently large) with homogeneous Dirichlet boundary conditions.
For simplicity of notation we introduce the method for the case of one spatial dimension (d = 1) defined over the interval (a, b) with homogeneous Dirichlet boundary conditions. Generalization to higher dimension are straightforward for tensor product grids, and the results remain valid without modifications. For d = 1, we choose the spatial mesh size h = ∆x > 0 with h = (b − a)/M for M an even positive integer, and let the grid points be In the discretization, we use sine-pseudospectral method for spatial derivatives and backward/forward Euler scheme for linear/nonlinear terms in time discretization. The gradient flow (2.1)-(2.3) is discretized, for j = 1, 2, . . . , M − 1 and n ≥ 1, as Here, D s xx , a pseudospectral differential operator approximation of ∂ xx , is defined as and α l (l = −1, 0, 1) are the stabilization parameters which are chosen in the 'optimal' form (such that the time step can be chosen as large as possible) as [4] (3.7) The homogeneous Dirichlet boundary conditions are discretized as The projection step (2.4)-(2.4) is discretized, for 0 ≤ j ≤ M and n ≥ 1, as The initial data (A.10) is discretized as φ 0 l,j = φ l (x j , 0), j = 0, 1, 2, · · · , M, l = −1, 0, 1.
The linear system (3.1)-(3.3) can be solved very efficiently by using the fast sine transform. In fact, take discrete sine transform at both sides, we get Solve the above system in the phase space, we obtain Remark 3.1. The gradient flow (2.1)-(2.3) can also be discretized by using the backward Euler finite difference method proposed in [5] or the backward Euler sine-pseudospectral method proposed in [4] for computing the ground state of onecomponent BEC. . We first test that the converged solution is independent of different choices of the initial data in (A.10) and energy diminishing property of the normalized gradient flow. In order to do so, we choose the initial data in (A.10) as • Gaussian profiles satisfying the constraints in (1.16) initially, i.e.
where κ is a constant satisfying 0 < κ < 1 − |M |. • Unnormalized Gaussian profiles, i.e.  time step ∆t is small. Based on our numerical experiments, for 0 ≤ M ≤ 1, we suggest the initial data in (A.10) be chosen as: i) with ferromagnetic interaction, i.e. β s ≤ 0 and ii) with antiferromagnetic interaction, i.e. β s > 0 where φ ap g (x) can be chosen as the approximate ground state solution of single component BEC, e.g. the harmonic oscillator approximation when β n small and the Thomas-Fermi approximation when β n ≫ 1 [5,8,7]. Based on these choices of initial data, we report the ground states computed by our numerical method.  i.e. β s ≤ 0, for fixed number of particles N in the condensate, when the magnetization M increases from 0 to 1, the mass N 1 increases from 0.25 to 1, the mass N −1 decreases from 0.25 to 0 and the mass N 0 decreases from 0.5 to 0 (c.f. Fig. 4.9a); while for antiferromagnetic interaction, i.e. β s ≥ 0, N 1 increases from 0.5 to 1, N −1 decreases from 0.5 to 0 and N 0 = 0 (c.f. Fig. 4.9b)   chemical potentials increase (c.f. Fig. 4.8). These observations agree with those obtained in [9] and [30] by different numerical methods.

Application in 1D
with optical lattice potential. In this subsection, our method is applied to compute the ground state of spin-1 BEC in one dimension (1D) with an optical lattice potential. Again, two different interaction are considered: •    to that, similar conclusions as those in the end of previous subsection can also be observed in this case.

Applications in 3D
with optical lattice potential. In this subsection, our method is applied to compute the ground state of spin-1 BEC in three dimensions (3D) with an optical lattice potential. Again, two different interaction are considered: • Case I. For 87 Rb with dimensionless quantities in (1.9)-(1.11) used as: d = 3, V (x) = 1 2 x 2 + y 2 + z 2 + 100 sin 2 πx   From Fig. 4.12, we can see that our method can be used to compute the ground state of spin-1 BEC in 3D with general trapping potential.

Conclusions.
We have proposed an efficient and accurate normalized gradient flow or imaginary time method to compute the ground state of spin-1 Bose-Einstein condensates by introducing a third normalization condition, in addition to the conservation of total particle number and the conservation of total magnetization. The condition is derived from the relationships between the chemical potentials of the three spinor components together with a splitting scheme applied to the continuous normalized gradient flows proposed to compute the ground state of spin-1 BEC. The backward-forward sine-pseudospectral method is applied to discretize the normalized gradient flow for practical computation. The ground state solutions and fraction of each component are reported for both ferromagnetic and antiferromagnetic interaction cases. The energy and chemical potentials of the condensate are also reported. In addition, the method may be further extended to other spinor condensate with higher degree of freedom as well as spinor condensate in the presence of external magnetic field, which will be our future study.
µ Φ (t) and λ Φ (t) are chosen such that the above CNGF is mass (or normalization) and magnetization conservative and they are given as [9]   with For the above CNGF, for any given initial data it was proven that the total mass and magnetization are conservative and the energy is diminishing [9], i.e.