# Computing ground states of spin-2 Bose-Einstein condensates by the
normalized gradient flow††thanks: Funding: The work was partially supported by
the Ministry of Education of Singapore under its Academic Research Fund
MOE-T2EP20122-0002 (A-8000962-00-00) (W. Bao), the National Natural Science
Foundation of China Grant No. 11971335 (Q. Tang) and No. 11971007 (Y. Yuan).

Weizhu Bao Department of Mathematics, National University of Singapore,
Singapore 119076<EMAIL_ADDRESS>URL: https://blog.nus.edu.sg/matbwz/).
Qinglin Tang  School of Mathematics, Sichuan University, Chengdu 610064, P. R.
China<EMAIL_ADDRESS>URL: http://www.qinglin-tang-scu.com). Yongjun
Yuan LCSM (MOE), School of Mathematics and Statistics, Hunan Normal
University, Changsha, Hunan 410081, P. R. China (yyj1983@hunnu.edu.cn).

###### Abstract

We propose and analyze an efficient and accurate numerical method for
computing ground states of spin-2 Bose-Einstein condensates (BECs) by using
the normalized gradient flow (NGF). In order to successfully extend the NGF to
spin-2 BECs which has five components in the vector wave function but with
only two physical constraints on total mass conservation and magnetization
conservation, two important techniques are introduced for designing the
proposed numerical method. The first one is to systematically investigate the
ground state structure and property of spin-2 BECs within a spatially uniform
system, which can be used on how to properly choose initial data in the NGF
for computing ground states of spin-2 BECs. The second one is to introduce
three additional projection conditions based on the relations between the
chemical potentials, together with the two existing physical constraints, such
that the five projection parameters used in the projection step of the NGF can
be uniquely determined. Then a backward-forward Euler finite difference method
is adapted to discretize the NGF. We prove rigorously that there exists a
unique solution of the nonlinear system for determining the five projection
parameters in the full discretization of the NGF under a mild condition on the
time step size. Extensive numerical results on ground states of spin-2 BECs
with ferromagnetic/nematic/cyclic phase and harmonic/optical lattice/box
potential in one/two dimensions are reported to show the efficiency and
accuracy of the proposed numerical method and to demonstrate several
interesting physical phenomena on ground states of spin-2 BECs.

###### keywords:

spin-2 Bose-Einstein condensate, coupled Gross-Pitaevskii equation, ground
state, normalized gradient flow, backward-forward Euler finite difference
method.

###### AMS:

35Q55, 49J45, 65N06, 65N12, 65Z05, 81-08

##  1 Introduction

Since its experimental realization in 1995 [2, 17, 27], the Bose-Einstein
condensate (BEC) has stimulated great excitement in the physical community and
regains vast interests in atomic and molecular as well as condense matter
physics. At early stage, atoms were magnetically trapped in BEC experiments
and hence their spin degrees of freedom were frozen [2, 17, 27]. Nevertheless,
recently developed optical trapping techniques [40] have enabled to release
the spin internal degrees of freedom, opening up a new research arena of
quantum many-body systems named spinor BEC [28, 21, 36]. Extensive theoretical
and experimental studies have been carried out to reveal numerous new quantum
phenomena which are generally absent in a spin-frozen condensate [16, 24, 29,
32, 33, 31, 30, 38, 44, 46, 47]. Within the mean-field approximation, in
contrast with a spin-frozen BEC whose order parameter can be well described by
a single-component Gross-Pitaevskii equation (GPE), a spin-$F$
$(F\in\mathbb{N})$ BEC is described by a coupled Gross-Pitaevskii equation
(CGPE) consisting of $2F+1$ equations, each of which governs one of the $2F+1$
hyperfine states $(m_{F}=-F,\cdots,F)$. In this paper, we consider spin-2
BECs, i.e. $F=2$.

At temperatures $T$ much smaller than the critical temperature $T_{c}$, a
spin-2 BEC can be well described by the macroscopic complex-valued vector wave
function
$\Psi:=\Psi(\bm{x},t)=\big{(}\psi_{2}(\bm{x},t),\psi_{1}(\bm{x},t),\psi_{0}(\bm{x},t),\psi_{-1}(\bm{x},t),\psi_{-2}(\bm{x},t)\big{)}^{T}:=\big{(}\psi_{2},\psi_{1},\psi_{0},\psi_{-1},\psi_{-2})^{T}$
that satisfies the following dimensionless CGPE [7, 42, 43, 30]:

(1) $\begin{split}&i\partial_{t}\psi_{\pm 2}=\left(H+\beta_{0}\rho\pm
2\beta_{1}F_{z}\right)\psi_{\pm 2}+\beta_{1}F_{\mp}\psi_{\pm
1}+\frac{\beta_{2}}{\sqrt{5}}A_{00}\bar{\psi}_{\mp 2},\\\
&i\partial_{t}\psi_{\pm 1}=\left(H+\beta_{0}\rho\pm\beta_{1}F_{z}\
\right)\psi_{\pm
1}+\beta_{1}\big{(}\frac{\sqrt{6}}{2}F_{\mp}\psi_{0}+F_{\pm}\psi_{\pm
2}\big{)}-\frac{\beta_{2}}{\sqrt{5}}A_{00}\bar{\psi}_{\mp 1},\\\
&i\partial_{t}\psi_{0}=\left(H+\beta_{0}\rho\right)\psi_{0}+\frac{\sqrt{6}\beta_{1}}{2}\big{(}F_{+}\psi_{1}+F_{-}\psi_{-1}\big{)}+\frac{\beta_{2}}{\sqrt{5}}A_{00}\bar{\psi}_{0};\end{split}$

where

(2) $\begin{split}&H=-\frac{1}{2}\nabla^{2}+V(\bm{x}),\quad
F_{z}:=F_{z}(\Psi)=2\big{(}|\psi_{2}|^{2}-|\psi_{-2}|^{2}\big{)}+|\psi_{1}|^{2}-|\psi_{-1}|^{2},\\\
&\rho(\bm{x},t)=\sum_{\ell=-2}^{2}|\psi_{\ell}(\bm{x},t)|^{2},\quad
A_{00}:=A_{00}(\Psi)=\frac{1}{\sqrt{5}}\left[2\psi_{2}\psi_{-2}-2\psi_{1}\psi_{-1}+\psi_{0}^{2}\right],\\\
&F_{+}:=F_{+}(\Psi)=\bar{F}_{-}:=\bar{F}_{-}(\Psi)=2\big{(}\bar{\psi}_{2}\psi_{1}+\bar{\psi}_{-1}\psi_{-2}\big{)}+\sqrt{6}\big{(}\bar{\psi}_{1}\psi_{0}+\bar{\psi}_{0}\psi_{-1}\big{)}.\end{split}$

Here, $\bm{x}\in\mathbb{R}^{d}$ ($d=1,2,3$) is the Cartesian coordinate
vector, $t$ is time, $\rho:=\rho(\bm{x},t)$ is the total density, $\beta_{0}$,
$\beta_{1}$, $\beta_{2}$ are real constants characterizing the spin-
independent interaction, spin-exchange interaction and spin-singlet
interaction, respectively, and $\bar{f}$ denotes the complex conjugate of $f$.
The real function $V(\bm{x})$ represents the external trapping potential.

There are three important invariants of the CGPE (1), i.e. the mass (or
normalization) of the wave function

(3)
$\mathcal{N}(\Psi(\cdot,t)):=\|\Psi(\cdot,t)\|^{2}=\sum_{\ell=-2}^{2}\int_{\mathbb{R}^{d}}|\psi_{\ell}(\bm{x},t)|^{2}d\bm{x}\equiv\mathcal{N}(\Psi(\cdot,0))=1,\quad
t\geq 0,$

the magnetization

(4)
$\mathcal{M}(\Psi(\cdot,t)):=\sum_{\ell=-2}^{2}\int_{\mathbb{R}^{d}}\ell|\psi_{\ell}(\bm{x},t)|^{2}d\bm{x}\equiv\mathcal{M}(\Psi(\cdot,0))=M,\quad
t\geq 0,$

where $M\in[-2,2]$ is a given constant, and the energy per particle

(5) $\begin{split}\mathcal{E}(\Psi(\cdot,t)):=&\int_{\mathbb{R}^{d}}\bigg{[}\
\sum_{\ell=-2}^{2}\left(\frac{1}{2}|\nabla\psi_{\ell}|^{2}+V(\bm{x})|\psi_{\ell}|^{2}\right)+\frac{\beta_{0}}{2}\rho^{2}+\frac{\beta_{2}}{2}|A_{00}(\Psi)|^{2}\\\
&+\frac{\beta_{1}}{2}\bigg{(}|F_{+}(\Psi)|^{2}+|F_{z}(\Psi)|^{2}\bigg{)}\bigg{]}\;d\bm{x}\equiv\mathcal{E}(\Psi(\cdot,0)),\quad
t\geq 0.\end{split}$

One important problem in the theoretical study of a spin-2 BEC is to find its
ground state so as to initialize its dynamics and to predict new important
phase categories of the ground state which can be later compared/confirmed
with/by those physical experimental observations. In [7], the authors carried
out a theoretical study on the existence and uniqueness of the ground state.
Meanwhile, validity of the so called single mode approximation (SMA) of the
ground state (which simplifies the ground state computation) is partially
investigated. Different numerical methods have been proposed to compute the
ground state of a scalar BEC [1, 3, 4, 5, 8, 9, 11, 14, 15, 19, 20, 22, 23,
25, 26, 45] and a spin-1 BEC [10, 13, 15, 18, 35, 48, 49, 50]. Among them, a
simple and most popular method is the normalized gradient flow (NGF) (or
imaginary time method) incorporated with a proper discretization scheme to
evolve the resulted gradient flow under normalization of the wave function [5,
11, 13, 15]. However, to extend the NGF from single-component BEC and spin-1
BEC to spin-2 BEC, due to the fact that we only have two conservation
conditions (3) and (4) and there are five projection constants to be
determined in the projection step, it is unclear that the NGF method could be
easily and straightforwardly extended to compute ground states of spin-2 BECs.
A projection gradient method [42, 43] was proposed to compute ground states of
spin-2 BECs, where a continuous normalized gradient flow (CNGF) was
discretized by the Crank-Nicolson finite difference method with a proper and
special way to deal with the nonlinear terms. This scheme was proved to be
energy-diminishing and conserve both the total mass and magnetization in the
discrete level. However, a fully nonlinear coupled system need to be solved at
each time step, which introduces much computational cost, especially in three
dimensions. Recently, numerical methods were presented for computing ground
states of spin-2 BEC based on the NGF with inaccurate projections [18, 41].
The main objective of this paper is to present and analyze a numerical method
for computing ground states of spin-2 BECs via the NGF. In order to do so, two
main techniques are presented, which are: (i) to carry out a systematic study
on the ground state structure and property of a spin-2 BEC within a spatially
uniform system, which is then used to choose simple and proper initial data in
the NGF for computing ground states of spin-2 BECs, and (ii) to introduce
three additional projection conditions based on the relations between the
chemical potentials of a spin-2 BEC for overcoming the fact that there are
five components in the vector wave function but with only two constraints on
total mass conservation and total magnetization conservation. In fact, the
proposed three additional projection conditions, together with total mass
conservation and magnetization conservation, can completely determine the five
projection constants used in projection step of the NGF. This enables us to
extend the simple and powerful NGF method to compute ground states of spin-2
BECs.

The rest of the paper is organized as follows. In section 2, ground state
structure and property of spin-2 BEC in a spatially uniform system are
investigated systematically. In section 3, a NGF is constructed by introducing
three additional projection conditions and then a backward-forward Euler
finite difference method is presented to discretize the NGF. In section 4, we
report extensive numerical results on ground states of spin-2 BECs with
ferromagnetic/nematic/cyclic phase and under harmonic/optical lattice/box
potential in one/two dimensions. Finally, some conclusions are drawn in
section 5.

##  2 Ground states and their properties

In this section, we mainly investigate the ground state structure and property
of a spin-2 BEC.

###  2.1 Euler-Lagrange equations and classification of ground states

The ground state of a spin-2 BEC is defined as the minimizer of the following
nonconvex minimization problem:

(6) $\Phi_{g}:=\arg\min_{\Phi\in\mathbb{S}}\mathcal{E}(\Phi),$

where the nonconvex set $\mathbb{S}$ is defined as

(7)
$\mathbb{S}=\left\\{\Phi=(\phi_{2},\phi_{1},\phi_{0},\phi_{-1},\phi_{-2})^{T}\
|\ \mathcal{N}(\Phi)=1,\ \mathcal{M}(\Phi)=M,\
\mathcal{E}(\Phi)<\infty\right\\}.$

Define the Lagrangian

(8)
$\mathcal{L}(\Phi,\mu,\lambda):=\mathcal{E}(\Phi)-\mu\left(\sum_{\ell=-2}^{2}\|\phi_{\ell}\|^{2}-1\right)-\lambda\left(\sum_{\ell=-2}^{2}\ell\|\phi_{\ell}\|^{2}-M\right)$

with $\|\phi\|^{2}:=\int_{\mathbb{R}^{d}}|\phi(\bm{x})|^{2}d\bm{x}$, then the
Euler-Lagrange equation associated to the constrained minimization problem (6)
is given by

(9) $\begin{split}&(\mu\pm 2\lambda)\phi_{\pm 2}=\left(H+\beta_{0}\rho\pm
2\beta_{1}F_{z}\right)\phi_{\pm 2}+\beta_{1}F_{\mp}\phi_{\pm
1}+\frac{\beta_{2}}{\sqrt{5}}A_{00}\bar{\phi}_{\mp 2},\\\
&(\mu\pm\lambda)\phi_{\pm
1}=\left(H+\beta_{0}\rho\pm\beta_{1}F_{z}\right)\phi_{\pm
1}+\beta_{1}\big{(}\frac{\sqrt{6}}{2}F_{\mp}\phi_{0}+F_{\pm}\phi_{\pm
2}\big{)}-\frac{\beta_{2}}{\sqrt{5}}A_{00}\bar{\phi}_{\mp 1},\\\
&\mu\phi_{0}=\left(H+\beta_{0}\rho\right)\phi_{0}+\frac{\sqrt{6}}{2}\beta_{1}\big{(}F_{+}\phi_{1}+F_{-}\phi_{-1}\big{)}+\frac{\beta_{2}}{\sqrt{5}}A_{00}\bar{\phi}_{0};\end{split}$

under the following two constraints

(10)
$\mathcal{N}(\Phi):=\|\Phi\|^{2}:=\sum_{\ell=-2}^{2}\|\phi_{\ell}\|^{2}=1,\quad\mathcal{M}(\Phi):=\sum_{\ell=-2}^{2}\ell\|\phi_{\ell}\|^{2}=M,\
M\in[-2,2].$

Here $\mu$ and $\lambda$ are the Lagrange multipliers (or chemical potentials)
of the CGPE (1). In fact, the problem of (9) with (10) is also a constrained
nonlinear eigenvalue problem, or time-independent CGPE. Thus the ground state
of (6) can also be viewed as the eigenfunction of the constrained nonlinear
eigenvalue problem (9) with (10), which has the lowest energy among all
eigenfunctions. Other eigenfunctions with higher energy are usually called as
excited states of the CGPE (1).

Existence and uniqueness of the ground state of (6) was carried out
mathematically in [7, 30, 31], and the uniqueness is understood in the sense
that
$\Phi_{g}^{1}:=(\phi_{2}^{1,g},\phi_{1}^{1,g},\phi_{0}^{1,g},\phi_{-1}^{1,g},\phi_{-2}^{1,g})^{T}$
and
$\Phi_{g}^{2}:=(\phi_{2}^{2,g},\phi_{1}^{2,g},\phi_{0}^{2,g},\phi_{-1}^{2,g},\phi_{-2}^{2,g})^{T}$
are regarded as the same ground state if there exists a
$\bm{\alpha}:=(\alpha_{2},\alpha_{1},\alpha_{0},\alpha_{-1},\alpha_{-2})^{T}=\left(e^{i(2\theta_{1}-\theta_{0})},\,e^{i\theta_{1}},\,e^{i\theta_{0}}\right.$,
$\left.e^{i(2\theta_{0}-\theta_{1})},\,e^{i(3\theta_{0}-2\theta_{1})}\right)^{T}$
with $\theta_{0},\theta_{1}\in\mathbb{R}$ such that
$\phi_{\ell}^{1,g}=\alpha_{\ell}\phi_{\ell}^{2,g}$ for $\ell=-2,-1,0,1,2$ [7,
30, 31]. According to [24, 30], the ground state $\Phi_{g}$ of (6) can be
classified into three categories based on the values of $|F_{+}(\Phi_{g})|$
and $|A_{00}(\Phi_{g})|$: (i) ferromagnetic if $|F_{+}(\Phi_{g})|>0$ and
$A_{00}(\Phi_{g})=0$, (ii) nematic if $F_{+}(\Phi_{g})=0$ and
$|A_{00}(\Phi_{g})|>0$, and (iii) cyclic if
$F_{+}(\Phi_{g})=A_{00}(\Phi_{g})=0$. We remark here that, when $M=\pm 2$, the
CGPE (1) reduces to a single component GPE with $\psi_{\pm 2}$, respectively,
and thus the ground state computation is reduced to a single component BEC,
which has been widely studied in the literatures [9, 11, 12, 14]. In addition,
if one replaces the magnetization $\mathcal{M}(\Phi)=M$ in (7) by
$\mathcal{M}(\Phi)=-M$, it is easy to see that the ground state
$\Phi_{g}:=(\phi_{2}^{g},\phi_{1}^{g},\phi_{0}^{g},\phi_{-1}^{g},\phi_{-2}^{g})^{T}$
of (6) can be simply replaced by
$\tilde{\Phi}_{g}:=(\phi_{-2}^{g},\phi_{-1}^{g},\phi_{0}^{g},\phi_{1}^{g},\phi_{2}^{g})^{T}$.
Thus for the simplicity of notations and presentation, from now on, we only
consider the magnetization $M\in[0,2)$.

######  Remark 2.1.

In the literature [7, 30], instead of the ground state defined as the
minimizer of the energy function $\mathcal{E}(\Phi)$ under two constraints of
the total mass conservation $\mathcal{N}(\Phi)=1$ and total magnetization
$\mathcal{M}(\Phi)=M$ with $M\in[-2,2]$, i.e. (6), another type of ground
state has also been studied, which is defined as the minimizer of the energy
function $\mathcal{E}(\Phi)$ under only one constraint of the total mass
conservation $\mathcal{N}(\Phi)=1$, i.e.

(11) $\tilde{\Phi}_{g}:=\arg\min_{\|\Phi\|=1}\mathcal{E}(\Phi).$

In fact, the above minimization problem can be obtained via the minimization
problem (6) by further minimizing for $M\in[-2,2]$, i.e.

(12) $\tilde{\Phi}_{g}:=\arg\min_{M\in[-2,2]}\Phi_{g}^{M},\quad\hbox{\rm
with}\quad\Phi_{g}^{M}:=\arg\min_{\Phi\in\mathbb{S}}\mathcal{E}(\Phi).$

The Euler-Lagrange equation associated to the minimization problem (11) is
given by

(13) $\begin{split}&\mu\phi_{\pm 2}=\left(H+\beta_{0}\rho\pm
2\beta_{1}F_{z}\right)\phi_{\pm 2}+\beta_{1}F_{\mp}\phi_{\pm
1}+\frac{\beta_{2}}{\sqrt{5}}A_{00}\bar{\phi}_{\mp 2},\\\ &\mu\phi_{\pm
1}=\left(H+\beta_{0}\rho\pm\beta_{1}F_{z}\right)\phi_{\pm
1}+\beta_{1}\big{(}\frac{\sqrt{6}}{2}F_{\mp}\phi_{0}+F_{\pm}\phi_{\pm
2}\big{)}-\frac{\beta_{2}}{\sqrt{5}}A_{00}\bar{\phi}_{\mp 1},\\\
&\mu\phi_{0}=\left(H+\beta_{0}\rho\right)\phi_{0}+\frac{\sqrt{6}}{2}\beta_{1}\big{(}F_{+}\phi_{1}+F_{-}\phi_{-1}\big{)}+\frac{\beta_{2}}{\sqrt{5}}A_{00}\bar{\phi}_{0};\end{split}$

under only one constraint $\|\Phi\|=1$. Here $\mu$ is the Lagrange multiplier
(or chemical potential) of the CGPE (1) under the total mass conservation (3).

###  2.2 Ground states in a spatially uniform system

Consider a spin-2 BEC in a spatially uniform system, i.e. the CGPE (1) without
potential on a bounded domain $\mathcal{D}$ with measure $|\mathcal{D}|=1$ and
periodic boundary condition. i.e., $V(\bm{x})\equiv 0$. In this special case,
all possible ground states are constants in the form of [30]

(14)
$\Phi(\bm{x})\equiv\bm{\xi}=:(\xi_{2},\xi_{1},\xi_{0},\xi_{-1},\xi_{-2})^{T},\quad\bm{x}\in\overline{\mathcal{D}},$

where $\xi_{j}\in\mathbb{C}$ for $j=-2,-1,0,1,2$. Plugging (14) into (5) with
$V(\bm{x})\equiv 0$ and $\Psi(\cdot,t)=\bm{\xi}$ and replacing
$\mathbb{R}^{d}$ by $\mathcal{D}$, after a detailed computation, we obtain

(15)
$\mathcal{E}(\Phi)=\mathcal{E}_{U}(\bm{\xi}):=\frac{1}{2}\left[\left(\beta_{1}|\tau|^{2}+\frac{\beta_{2}}{5}|\delta|^{2}\right)+\bigg{(}\beta_{0}+\beta_{1}M^{2}\bigg{)}\right]=:E\left(\tau,\delta\right),$

where

(16)
$\left\\{\begin{array}[]{l}\tau=:\tau(\bm{\xi})=F_{+}(\bm{\xi})=2\big{(}\bar{\xi}_{2}\xi_{1}+\bar{\xi}_{-1}\xi_{-2}\big{)}+\sqrt{6}\big{(}\bar{\xi}_{1}\xi_{0}+\bar{\xi}_{0}\xi_{-1}\big{)},\\\\[5.0pt]
\delta=:\delta(\bm{\xi})=A_{00}(\bm{\xi})=2\xi_{2}\xi_{-2}-2\xi_{1}\xi_{-1}+\xi_{0}^{2}.\end{array}\right.$

Meanwhile, the conservation of mass (3) and magnetization  (4) leads to:

(17)
$\left\\{\begin{array}[]{lr}|\xi_{2}|^{2}+|\xi_{-2}|^{2}+|\xi_{0}|^{2}+|\xi_{1}|^{2}+|\xi_{-1}|^{2}=1,&\\\\[2.5pt]
2(|\xi_{2}|^{2}-|\xi_{-2}|^{2})+|\xi_{1}|^{2}-|\xi_{-1}|^{2}=M,&\end{array}\right.\quad\forall\,M\in[0,2).$

Solving (6)-(7) for the ground state $\Phi_{g}(\bm{x})$ is now equivalent to
solve the following minimization problem for minimizer
$\bm{\xi}_{g}=:(\xi^{g}_{2},\xi^{g}_{1},\xi^{g}_{0},\xi^{g}_{-1},\xi^{g}_{-2})^{T}$
as

(18)
$\quad\bm{\xi}_{g}:=\arg\min_{\bm{\xi}\in\mathbb{S}_{C}}\mathcal{E}_{U}(\bm{\xi}),\quad\mathbb{S}_{C}=\left\\{\bm{\xi}\in\mathbb{C}^{5}\
\bigg{|}\ \sum_{\ell=-2}^{2}|\xi_{l}|^{2}=1,\
\sum_{\ell=-2}^{2}\ell|\xi_{\ell}|^{2}=M\right\\}.$

In the following, we show that the minimization problem (18) on the complex
manifold $\mathbb{S}_{C}$ can be reduced to a minimization problem on the real
manifold $\mathbb{S}_{R}=\mathbb{R}^{5}\cap\mathbb{S}_{C}$.

######  Lemma 2.1.

If $\bm{\xi}\in\mathbb{R}^{5}$, then the system (17) has a real solution if
and only if

(19) $\tau^{2}(\bm{\xi})+4\,\delta^{2}(\bm{\xi})\leq 4-M^{2}.$

###### Proof.

Firstly, we prove the necessary condition, i.e. we assume
$\bm{\xi}\in\mathbb{R}^{5}$ and the system (17) has a real solution. By (16)
and (17), we have

(20) $(\xi_{2}-\xi_{-2})^{2}+(\xi_{1}+\xi_{-1})^{2}=1-\delta.$

Noticing $|\delta|\leq\sum_{\ell=-2}^{2}\xi_{\ell}^{2}=1$ and denoting
$p=\sqrt{1-\delta}\geq 0$, then there exists a constant $\theta\in[0,2\pi)$
such that

(21) $\left\\{\begin{array}[]{l}\xi_{-2}-\xi_{2}=p\cos{\theta},\\\\[5.0pt]
\xi_{1}+\xi_{-1}=p\sin{\theta},\end{array}\right.\quad\Longleftrightarrow\quad\left\\{\begin{array}[]{l}\xi_{-2}=p\cos{\theta}+\xi_{2},\\\\[5.0pt]
\xi_{-1}=p\sin{\theta}-\xi_{1}.\end{array}\right.$

We then prove (19) is valid by considering four different parameter cases. To
simplify the presentation, in the following, we only state the formula for
$\xi_{0}$, $\xi_{1}$ and $\xi_{2}$, and the expressions of $\xi_{-1}$ and
$\xi_{-2}$ can be obtained directly from (21).

Case (i). If $\delta=1$, then $p=0$. From (21), we obtain $\xi_{2}=\xi_{-2}$
and $\xi_{1}=-\xi_{-1}$. By (17) and (16), we have

(22) $M=\tau=0,\quad\Longrightarrow\quad\tau^{2}+4\delta^{2}\leq 4-M^{2}.$

Case (ii). If $\delta<1$ (and thus $p>0$) and $\cos\theta=0$. Plugging (21)
into (17), we have

(23)
$\displaystyle\xi_{1}=\frac{1}{2p}(p^{2}+M)\sin\theta,\qquad\quad\xi_{0}=\frac{1}{p\sqrt{6}}(\tau\sin\theta-2p\xi_{2}),$
(24) $\displaystyle
16p^{2}\xi_{2}^{2}-4p\,\sin\theta\,\tau\,\xi_{2}+3p^{4}-6p^{2}+3M^{2}+\tau^{2}=0.$

Note that (24) has a real solution if and only if the corresponding
discriminant

(25) $\Delta:=-48p^{2}(4\delta^{2}+\tau^{2}+4M^{2}-4)\geq 0,$

which immediately implies

(26) $\tau^{2}+4\delta^{2}\leq 4-4M^{2}\leq 4-M^{2}.$

Case (iii). If $\delta<1$ (and thus $p>0$) and $\sin\theta=0$. Similarly, we
obtain

(27)
$\xi_{2}=\frac{(2p^{2}+M)\cos\theta}{4p},\quad\xi_{1}=\frac{\tau\cos\theta}{2p},\quad\xi_{0}^{2}=\frac{4-M^{2}-4\delta^{2}-4\tau^{2}}{8p^{2}}.$

Thus $\xi_{0}$ is real if and only if $4-M^{2}-4\delta^{2}-4\tau^{2}\geq 0$,
which again implies

(28) $\tau^{2}+4\delta^{2}\leq 4-M^{2}-3\tau^{2}\leq 4-M^{2}.$

Case (iv). If $\delta<1$ (and thus $p>0$) and $\sin(2\theta)\neq 0$. Plugging
(21) into (17), we have

(29)
$\displaystyle\qquad\xi_{2}=-(M+2p^{2}-p^{2}\sin^{2}\theta-2p\sin\theta\xi_{1})/(4p\cos\theta),$
(30)
$\displaystyle\qquad\xi_{0}=\frac{4-6\sin^{2}\theta}{\sqrt{6}\sin(2\theta)}\,\xi_{1}+\frac{M-2p^{2}+3p^{2}\sin^{2}\theta}{2\sqrt{6}p\cos\theta}+\frac{\tau}{\sqrt{6}p\sin\theta},$
(31) $\displaystyle\qquad A\xi_{1}^{2}+B\xi_{1}+C=0,$

with

$\begin{split}A=&8\csc^{2}(2\theta),\\\
B=&\frac{\csc^{2}\theta\sec^{2}\theta}{4p}\Big{[}5\tau\cos\theta+3\tau\cos(3\theta)-5M\sin\theta-8p^{2}\sin\theta+3M\sin(3\theta)\Big{]},\\\
C=&\frac{1}{2p^{2}}\Big{[}\tau^{2}\csc^{2}\theta+(M-2p^{2})\tau\csc\theta\sec\theta-\frac{1}{2}\sec^{2}\theta\big{(}-2M^{2}+6\delta
p^{2}-Mp^{2}+p^{4}\\\
&+3p^{2}(2\delta+M+p^{2})\cos(2\theta)-3p^{2}\tau\sin(2\theta)\big{)}\Big{]}.\end{split}$

Similarly, (31) has a real solution if and only if the corresponding
discriminant

(32)
$\Delta=-\frac{3\csc^{2}\theta\sec^{2}\theta}{p^{2}}\;\Big{[}3(M\sin\theta+\tau\cos\theta)^{2}+4\delta^{2}+\tau^{2}+M^{2}-4\Big{]}\geq
0,$

which implies

(33) $4-M^{2}-(4\delta^{2}+\tau^{2})\geq 3(M\sin\theta+\tau\cos\theta)^{2}\geq
0\quad\Longrightarrow\quad\tau^{2}+4\delta^{2}\leq 4-M^{2}.$

Secondly, we prove the sufficient condition, i.e. we assume (19) is valid.
When $\delta=1$, then $M=\tau=0$, and thus $\xi=(0,0,1,0,0)^{T}$ is a real
solution of (17); when $\delta<1$ and $M=0$, then (25) is satisfied, thus
(23)-(24) is a real solution; and when $\delta<1$ and $M\neq 0$, by choosing
$\theta=\pi-\arctan(\tau/M)$ such that $M\sin\theta=-\tau\cos\theta$ and (32)
is fulfilled, then (29)-(31) is a real solution. ∎

######  Lemma 2.2.

For $\forall\,\bm{\xi}\in\mathbb{S}_{C}$, we have

(34) $|\tau(\bm{\xi})|^{2}+4\,|\delta(\bm{\xi})|^{2}\leq 4-M^{2}.$

Hence there exists a $\bm{\zeta}_{R}\in\mathbb{S}_{R}$ such that
$\mathcal{E}_{U}(\bm{\zeta}_{R})=\mathcal{E}_{U}(\bm{\xi})$. Therefore the
minimization problem (18) has at least a real ground state.

###### Proof.

We prove this lemma by considering two different cases.

Case (i). If $\delta(\bm{\xi})=0$, noticing that

(35)
$\bm{\xi}\in\mathbb{S}_{C}\qquad\Longrightarrow\qquad\bm{\zeta}:=\left(|\xi_{2}|,|\xi_{1}|,|\xi_{0}|,|\xi_{-1}|,|\xi_{-2}|\right)^{T}\in\mathbb{S}_{R},$

hence the system (17) is fulfilled for $\bm{\zeta}$. By lemma 2.1, we have

(36)
$|\tau(\bm{\xi})|^{2}+4|\delta(\bm{\xi})|^{2}=|\tau(\bm{\xi})|^{2}\leq\tau^{2}(\bm{\zeta})\leq\tau^{2}(\bm{\zeta})+4\delta^{2}(\bm{\zeta})\leq
4-M^{2}.$

Case (ii). If $\delta(\bm{\xi})\neq 0$, noticing that

(37)
$\begin{split}\max_{\bm{\xi}\in\mathbb{S}_{C}}\left\\{|\tau(\bm{\xi})|^{2}+4|\delta(\bm{\xi})|^{2}+M^{2}\right\\}=&\max_{\bm{\xi}\in\mathbb{S}_{C}}\left\\{|\tau(\bm{\xi})|^{2}+4|\delta(\bm{\xi})|^{2}+F_{z}^{2}(\bm{\xi})\right\\}\\\
\leq&\max_{\bm{\zeta}\in\mathbb{S}_{1}}\left\\{|\tau(\bm{\zeta})|^{2}+4|\delta(\bm{\zeta})|^{2}+F_{z}^{2}(\bm{\zeta})\right\\},\end{split}$

with $\mathbb{S}_{1}=\left\\{\bm{\zeta}\in\mathbb{C}^{5}\ |\
\sum_{\ell=-2}^{2}|\zeta_{l}|^{2}=1\right\\}$, thus in order to prove the
inequality (34), we only need to show

(38)
$\max_{\bm{\zeta}\in\mathbb{S}_{1}}\left\\{|\tau(\bm{\zeta})|^{2}+4|\delta(\bm{\zeta})|^{2}+F_{z}^{2}(\bm{\zeta})\right\\}\leq
4.$

Consider an auxiliary minimization problem

(39)
$\bm{\zeta}_{g}:=(\zeta_{2}^{g},\zeta_{1}^{g},\zeta^{g}_{0},\zeta_{-1}^{g},\zeta_{-2}^{g})^{T}=\arg\min_{\bm{\zeta}\in\mathbb{S}_{1}}\mathcal{F}(\bm{\zeta}),$

where the auxiliary functional $\mathcal{F}(\bm{\zeta})$ is defined as

(40)
$\mathcal{F}(\bm{\zeta})=-\big{[}|\tau(\bm{\zeta})|^{2}+4|\delta(\bm{\zeta})|^{2}+F_{z}^{2}(\bm{\zeta})\big{]}.$

It is clear that

(41)
$\mathcal{F}(\bm{\zeta}_{g})=-\max_{\bm{\zeta}\in\mathbb{S}_{1}}\left\\{|\tau(\bm{\zeta})|^{2}+4|\delta(\bm{\zeta})|^{2}+F_{z}^{2}(\bm{\zeta})\right\\},$

and $\bm{\zeta}_{g}$ satisfies the following Euler-Lagrange equation
associated with problem (39) as
$\nabla_{\bar{\bm{\zeta}}}\,\mathcal{F}(\bm{\zeta})=\lambda_{\bm{\zeta}}\,\bm{\zeta}$
with $\lambda_{\bm{\zeta}}\in\mathbb{R}$ being the Lagrange multiplier and
$\bar{\bm{\zeta}}$ being the complex conjugate of $\bm{\zeta}$, i.e.
$\nabla_{\bar{\bm{\zeta}}}\,\mathcal{F}(\bm{\zeta}_{g})=\lambda_{\bm{\zeta}_{g}}\,\bm{\zeta}_{g}$.
Hence, by denoting
$\bm{\eta}_{g}=(\zeta^{g}_{-2},-\zeta^{g}_{-1},\zeta^{g}_{0},-\zeta^{g}_{1},\zeta^{g}_{2})^{T}$,
we have

(42)
$\left\\{\begin{array}[]{l}\bar{\bm{\zeta}_{g}}\cdot\nabla_{\bar{\bm{\zeta}}}\,\mathcal{F}(\bm{\zeta}_{g})=\lambda_{\bm{\zeta}_{g}},\\\\[5.0pt]
\bm{\eta}_{g}\cdot\nabla_{\bar{\bm{\zeta}}}\,\mathcal{F}(\bm{\zeta}_{g})=\lambda_{\bm{\zeta}_{g}}\,\bm{\eta}_{g}\cdot\bm{\zeta}_{g},\end{array}\right.\quad\Longrightarrow\;\quad\left\\{\begin{array}[]{l}-2(|\tau|^{2}+4|\delta|^{2}+F_{z}^{2})=\lambda_{\bm{\zeta}_{g}},\\\\[5.0pt]
(\lambda_{\bm{\zeta}_{g}}+8)\delta=0,\end{array}\right.$

which leads to

(43) $\left\\{\begin{array}[]{l}\lambda_{\bm{\zeta}_{g}}=-8,\\\\[5.0pt]
|\tau(\bm{\zeta}_{g})|^{2}+4|\delta(\bm{\zeta}_{g})|^{2}+F_{z}^{2}(\bm{\zeta}_{g})=4,\end{array}\right.\quad\Longrightarrow\;\quad\mathcal{F}(\bm{\zeta}_{g})=-4.$

Noticing (41), one gets (38). This, together with (37), concludes the desired
inequality (34).

Therefore, combining (34) and Lemma 2.1, one obtains that for
$\forall\,\bm{\xi}\in\mathbb{S}_{C}$, there exists a
$\bm{\zeta}_{R}\in\mathbb{S}_{R}$ such that

$\tau(\bm{\zeta}_{R})=|\tau(\bm{\xi})|,\quad\delta(\bm{\zeta}_{R})=|\delta(\bm{\xi})|\quad\mbox{and}\quad\mathcal{E}_{U}(\bm{\zeta}_{R})=\mathcal{E}_{U}(\bm{\xi}),$

which implies

(44)
$\min_{\bm{\xi}\in\mathbb{S}_{C}}\mathcal{E}_{U}(\bm{\xi})=\min_{\bm{\zeta}_{R}\in\mathbb{S}_{R}}\mathcal{E}_{U}(\bm{\zeta}_{R})=\min_{S_{\tau,\delta}}\;E(\tau,\delta)=\frac{1}{2}(\beta_{0}+\beta_{1}M^{2})+\frac{1}{2}\;\min_{S_{\tau,\delta}}\;f(\tau,\delta),$

where

(45)
$f(\tau,\delta)=\beta_{1}\tau^{2}+\frac{\beta_{2}}{5}\delta^{2}\quad\mbox{and}\quad
S_{\tau,\delta}:=\left\\{\tau\in{\mathbb{R}},\;\delta\in{\mathbb{R}},\;\tau^{2}+4\delta^{2}\leq
4-M^{2}\right\\}.$

Obviously, for any given $\beta_{1}$ and $\beta_{2}$, the minimization problem

(46) $(\tau_{g},\delta_{g}):=\arg\;\min_{S_{\tau,\delta}}\;f(\tau,\delta)$

possesses at least a solution. As a result, according to Lemma 2.1, the
minimization problem (18) has at least a real ground state. ∎

Thanks to Lemmas 2.1 & 2.2, noticing (15), to find the ground state
$\bm{\xi}_{g}$ of (18) is then reduced to find minimizers of the minimization
problem (46). Actually, it can be solved analytically since it is to find
minimizers of a quadratic function over an elliptic region. To illustrate
this, Fig. 1 shows contour plots of the energy $E(\tau,\delta)$ for
$\beta_{0}=M=0$ with different $\beta_{1}$ & $\beta_{2}$ on the elliptic
region (19). From the figure, one can obtain the minimizers
$(\tau_{g},\delta_{g})$. Then from the proof of Lemma 2.1, one can obtain
$p_{g}=\sqrt{1-\delta_{g}}$ and choose $\theta_{g}\in[0,2\pi)$ satisfying
$M\sin\theta_{g}=-\tau_{g}\cos\theta_{g}$, and finally one can obtain all real
ground states $\bm{\xi}_{g}$ of (18) based on the four different parameter
cases as shown in the proof of Lemma 2.1.

Fig. 1: Contour plots of the energy $E(\tau,\delta)$ in (46) for
$\beta_{0}=M=0$ with different $\beta_{1}=-1$ and $\beta_{2}=20$ (left),
$\beta_{1}=1$ and $\beta_{2}=-20$ (middle), and $\beta_{1}=1$ and
$\beta_{2}=20$ (right). ${\color[rgb]{1,0,0}`+^{\prime}}$ denote those points
where the minimum value of $E(\tau,\delta)$ are achieved. According to the
value (and position on the graph) of those points
${\color[rgb]{1,0,0}`+^{\prime}}$, one can immediately conclude the phase of
the corresponding ground state are in ferromagnetic (left), nematic (middle)
and cyclic (right).

Noticing that if $\bm{\xi}$ solves (16)–(17) with $\tau$, then
$\tilde{\bm{\xi}}=(\xi_{2},-\xi_{1},\xi_{0},-\xi_{-1},\xi_{-2})^{T}$ solves
(16)–(17) with the parameter $-\tau$. Thus, it suffices to assume $\tau\geq 0$
for simplicity hereafter. Combining Lemmas 2.1 & 2.2 and noticing (46), we
have the following theorem for the ground state of a spin-2 BEC in a spatially
uniform system, i.e. (18). Denote

(47) $m_{1}=\sqrt{2+M},\ m_{2}=\sqrt{2-M},\ m_{3}=\sqrt{1+M},\
m_{4}=\sqrt{1-M},\quad M\in[0,2).$

######  Theorem 2.1.

Assume $\bm{\xi}_{g}\in\mathbb{R}^{5}$ be the real ground state of (18), we
have

(i). When $\beta_{1}<0$ and $\beta_{2}>20\beta_{1}$, then $E(\tau,\delta)$
attains the minimum at $(\tau,\delta)=(\sqrt{4-M^{2}},0)$, i.e., the ground
state is ferromagnetic. Moreover, for $\forall\,M\in[0,2)$, $\bm{\xi}_{g}$ can
be taken as

(48)
$\bm{\xi}_{g}=\left(m_{1}^{4}/16,\;m_{1}^{3}m_{2}/8,\;\sqrt{6}m_{1}^{2}m_{2}^{2}/16,\;m_{1}m_{2}^{3}/8,\;m_{2}^{4}/16\right)^{T}.$

(ii). When $\beta_{2}<0$ and $\beta_{2}<20\beta_{1}$, then $E(\tau,\delta)$
attains the minimum at $(\tau,\delta)=(0,\pm\frac{\sqrt{4-M^{2}}}{2})$, i.e.,
the ground state is nematic. Moreover, if $0<M<2$, then $\bm{\xi}_{g}$ can be
taken as

(49) $\bm{\xi}_{g}=\left(m_{1}/2,\;0,\;0,\;0,\;m_{2}/2\right)^{T};$

and if $M=0$, $\bm{\xi}_{g}$ can be taken in two different types as

(50)
$\begin{split}&\bm{\xi}_{g}=\left(\gamma_{1}\cos\theta,\;\gamma_{1}\sin\theta,\;\gamma,\;-\gamma_{1}\sin\theta,\;\gamma_{1}\cos\theta\right)^{T},\quad\hbox{\rm
or}\\\
&\bm{\xi}_{g}=\left(\cos\theta/\sqrt{2},\;\sin\theta/\sqrt{2},\;0,\;\sin\theta/\sqrt{2},\;-\cos\theta/\sqrt{2}\right)^{T};\end{split}$

where $\gamma_{1}=\sqrt{\frac{1-\gamma}{2}}$ for $|\gamma|\leq 1$ and
$\theta\in[0,2\pi)$.

(iii). When $\beta_{1}>0$ and $\beta_{2}>0$, then $E(\tau,\delta)$ attains the
minimum at $(\tau,\delta)=(0,0)$, i.e., the ground state is cyclic. Moreover,
if $M\in[0,1]$, then $\bm{\xi}_{g}$ can be taken in three different types as

(51)
$\begin{split}&\bm{\xi}_{g}=\left(m_{1}^{2}/4,\;0,\;\sqrt{2}\,m_{1}m_{2}/4,\;0,\;m_{2}^{2}/4\right)^{T},\quad\hbox{\rm
or}\\\
&\bm{\xi}_{g}=\left(\sqrt{3}\,m_{3}m_{4}/4,\;m_{3}^{2}/2,\;-\sqrt{2}\,m_{3}m_{4}/4,\;m_{4}^{2}/2,\;\sqrt{3}\,m_{3}m_{4}/4\right)^{T},\quad\hbox{\rm
or}\end{split}$ (52)
$\bm{\xi}_{g}=(\xi^{g}_{2},\xi^{g}_{1},\xi^{g}_{0},\xi^{g}_{-1},\xi^{g}_{-2})^{T},$

where

(53)
$\left\\{\begin{array}[]{l}\xi^{g}_{0}=-\frac{3\sqrt{6}}{8}M\sin^{2}\theta\cos\theta\mp-\frac{\sqrt{2}}{8}\big{(}2\cot(2\theta)+\cot\theta\big{)}g(\theta),\\\\[5.0pt]
\xi^{g}_{1}=\frac{3}{4}M\sin^{3}\theta+\frac{1}{4}m_{2}^{2}\sin\theta\pm\frac{\sqrt{3}}{4}g(\theta),\qquad\quad\xi^{g}_{-1}=\sin\theta-\xi^{g}_{1},\\\\[5.0pt]
\xi^{g}_{2}=\frac{1}{8}\big{(}3M\sin^{2}\theta+2m_{1}^{2}\big{)}\cos\theta\mp\frac{\sqrt{3}}{8}g(\theta)\tan\theta,\qquad\xi^{g}_{-2}=\xi^{g}_{2}-\cos\theta,\end{array}\right.$

with
$g(\theta)=\sqrt{\big{(}m_{1}m_{2}-3M^{2}\sin^{2}\theta\big{)}\sin^{2}\theta\cos^{2}\theta}$
for $\theta\in(0,2\pi)$ satisfying $m_{1}m_{2}-3M^{2}\sin^{2}\theta\geq 0$ and
$|\sin\theta|\neq 1$; and if $M\in(1,2)$, then $\bm{\xi}_{g}$ can only be
taken as (52) with (53).

###### Proof.

Combining Lemmas 2.1 & 2.2, noticing (15), to find the ground state
$\bm{\xi}_{g}$ of (18) is then reduced to find minimizers of the minimization
problem (46). For $(\tau,\delta)\in S_{\tau,\delta}$, we adapt the following
elliptic-polar coordinates as

(54)
$\tau=\rho\cos\theta,\quad\delta=\frac{1}{2}\rho\sin\theta,\quad\theta\in[0,2\pi),\quad\rho\in\Big{[}0,\sqrt{4-M^{2}}\Big{]}.$

Plugging (54) into (45), we get

(55)
$\begin{split}f(\tau,\delta)&=\beta_{1}\tau^{2}+\frac{\beta_{2}}{5}\delta^{2}=\beta_{1}\rho^{2}\cos^{2}\theta+\frac{\beta_{2}}{20}\rho^{2}\sin^{2}\theta\\\
&=\beta_{1}\rho^{2}+\left(\frac{\beta_{2}}{20}-\beta_{1}\right)\rho^{2}\sin^{2}\theta\\\
&=\frac{\beta_{2}}{20}\rho^{2}+\left(\beta_{1}-\frac{\beta_{2}}{20}\right)\rho^{2}\cos^{2}\theta,\quad\theta\in[0,2\pi),\quad\rho\in\Big{[}0,\sqrt{4-M^{2}}\Big{]}.\end{split}$

Combining (55) and (46), we immediately obtain

(i). When $\beta_{1}<0$ and $\beta_{2}>20\beta_{1}$, then $f(\tau,\delta)$
(and thus $E(\tau,\delta)$) attains the minimum at
$(\tau_{g},\delta_{g})=(\sqrt{4-M^{2}},0)$. This, together with (16), we get
$|F_{+}(\bm{\xi}_{g})|=|\tau_{g}|=\sqrt{4-M^{2}}>0$ and
$A_{00}(\bm{\xi}_{g})=\delta_{g}=0$, i.e., the ground state $\bm{\xi}_{g}$ of
(18) is ferromagnetic. Then the ground state (48) can be obtained directly
from the results in Lemma 2.1.

(ii). When $\beta_{2}<0$ and $\beta_{2}<20\beta_{1}$, then $f(\tau,\delta)$
(and thus $E(\tau,\delta)$) attains the minimum at
$(\tau_{g},\delta_{g})=(0,\pm\frac{\sqrt{4-M^{2}}}{2})$. This, together with
(16), we get $F_{+}(\bm{\xi}_{g})=\tau_{g}=0$ and
$|A_{00}(\bm{\xi}_{g})|=|\delta_{g}|=\frac{\sqrt{4-M^{2}}}{2}>0$, i.e., the
ground state $\bm{\xi}_{g}$ of (18) is nematic. Then the ground states
(49)-(50) can be obtained straightforwardly from the results in Lemma 2.1.

(iii). When $\beta_{1}>0$ and $\beta_{2}>0$, then $f(\tau,\delta)$ (and thus
$E(\tau,\delta)$) attains the minimum at $(\tau_{g},\delta_{g})=(0,0)$. This,
together with (16), we get $F_{+}(\bm{\xi}_{g})=\tau_{g}=0$ and
$A_{00}(\bm{\xi}_{g})=\delta_{g}=0$, i.e., the ground state $\bm{\xi}_{g}$ of
(18) is cyclic. Then the ground states (51)-(52) can be obtained directly from
the results in Lemma 2.1. ∎

######  Remark 2.2.

For a spin-2 BEC in the spatially uniform system, i.e. under the ansatz (14),
the minimization problem (11) collapses to

(56)
$\tilde{\bm{\xi}}_{g}:=\arg\min_{|\bm{\xi}|=1}\mathcal{E}_{U}(\bm{\xi})=\arg\min_{M\in[-2,2]}\mathcal{E}_{U}(\bm{\xi}_{g}^{M}),\quad\hbox{\rm
with}\quad\bm{\xi}_{g}^{M}:=\arg\;\min_{\bm{\xi}\in\mathbb{S}_{C}}\mathcal{E}_{U}(\bm{\xi}).$

Define

(57) $\beta(M)=:\mathcal{E}_{U}(\bm{\xi}_{g}^{M}),\qquad M\in[-2,2].$

By the results in Theorem 2.1, we have for $M\in[-2,2]$

(58)
$\beta(M)=:\mathcal{E}_{U}(\bm{\xi}_{g}^{M})=\frac{1}{2}\left\\{\begin{array}[]{ll}\beta_{0}+4\beta_{1},&\beta_{1}<0\;\;\&\;\;\beta_{2}>20\beta_{1},\\\\[2.5pt]
\beta_{0}+\frac{\beta_{2}}{5}+\frac{\left(20\beta_{1}-\beta_{2}\right)M^{2}}{20},&\beta_{2}<0\;\;\&\;\;\beta_{2}<20\beta_{1},\\\\[2.5pt]
\beta_{0}+\beta_{1}M^{2},&\beta_{1}>0\;\;\&\;\;\beta_{2}>0.\end{array}\right.$

Therefore, the energy of the ground state $\tilde{\bm{\xi}}_{g}$ defined in
the minimization problem (56) is given as

(59)
$\mathcal{E}_{U}(\tilde{\bm{\xi}}_{g})=\min_{M\in[-2,2]}\beta(M)=\frac{1}{2}\left\\{\begin{array}[]{ll}\beta_{0}+4\beta_{1},&\beta_{1}<0\;\;\&\;\;\beta_{2}>20\beta_{1},\\\
\beta_{0}+\beta_{2}/5,&\beta_{2}<0\;\;\&\;\;\beta_{2}<20\beta_{1},\\\
\beta_{0},&\beta_{1}>0\;\;\&\;\;\beta_{2}>0.\end{array}\right.$

This immediately suggests that, for the ground state of a spin-2 BEC defined
as the minimizer of the energy functional under the total mass conservation,
the ground state energy in a spatially uniform system is achieved in nematic
and cyclic phases when $M=0$, and respectively, in the ferromagnetic phase for
any $M\in[-2,2]$. The ground state is not unique. In fact, in the
ferromagnetic phase, the ground state $\tilde{\bm{\xi}}_{g}$ can be taken as
(48) for any $M\in[0,2]$, and in nematic and cyclic phases,
$\tilde{\bm{\xi}}_{g}$ can be taken as (50) and (51)-(52) with $M=0$,
respectively.

###  2.3 Single mode approximation of ground states

In the literature [7, 30], the single mode approximation (SMA) is an
interesting and useful tool for obtaining approximate ground states of spinor
BEC. It can reduce solving the ground state of a spin-2 BEC to solving the
ground state of a single component BEC. In fact, in the SMA for a spin-2 BEC,
one assumes an ansatz for $\Phi\in{\mathbb{S}}$ in (6) as

(60)
$\Phi(\bm{x})=\phi(\bm{x})(\xi_{2},\;\xi_{1},\;\xi_{0},\;\xi_{-1},\;\xi_{-2})^{T}=\phi(\bm{x})\,\bm{\xi}:=\Phi_{\rm
sma}(\bm{x}),$

where
$\bm{\xi}=(\xi_{2},\;\xi_{1},\;\xi_{0},\;\xi_{-1},\;\xi_{-2})^{T}\in\mathbb{S}_{C}$
and
$\phi:=\phi(\bm{x})\in\widetilde{\mathbb{S}}_{1}=:\\{\varphi|\int_{\mathbb{R}^{d}}|\varphi(\bm{x})|^{2}d\bm{x}=1\\}.$
Plugging (60) into (5) with $\Psi=\Phi$, noticing (15), we obtain

(61) $\mathcal{E}(\Phi)=\mathcal{E}(\Phi_{\rm
sma})=\int_{\mathbb{R}^{d}}\left[\frac{1}{2}|\nabla\phi|^{2}+V|\phi|^{2}+\mathcal{E}_{U}(\bm{\xi})|\phi|^{4}\right]d\bm{x}:=\mathcal{E}_{\rm
sma}(\phi,\bm{\xi}).$

Then the single mode approximation of the ground state is to find
$\Phi^{g}_{\rm sma}=\phi_{g}\,\bm{\xi}_{g}$ with
$\bm{\xi}_{g}\in\mathbb{S}_{C}$ and $\phi_{g}\in\widetilde{\mathbb{S}}_{1}$
such that

(62) $\Phi^{g}_{\rm sma}:=\arg\;\min_{\Phi_{\rm
sma}\in\mathbb{S}}\mathcal{E}(\Phi_{\rm sma}).$

Combining (62), (61), (56) and (57), we get

(63) $\begin{split}\mathcal{E}(\Phi^{g}_{\rm sma})=&\min_{\Phi_{\rm
sma}\in\mathbb{S}}\mathcal{E}(\Phi_{\rm
sma})=\min_{\phi\in\widetilde{\mathbb{S}}_{1}}\min_{\bm{\xi}\in\mathbb{S}_{C}}\mathcal{E}_{\rm
sma}(\phi,\bm{\xi})\\\
=&\min_{\phi\in\widetilde{\mathbb{S}}_{1}}\left\\{\int_{\mathbb{R}^{d}}\left[\frac{1}{2}|\nabla\phi|^{2}+V|\phi|^{2}+\left(\min_{\bm{\xi}\in\mathbb{S}_{C}}\mathcal{E}_{U}(\bm{\xi})\right)\,|\phi|^{4}\right]d\bm{x}\right\\}\\\
=&\min_{\phi\in\widetilde{\mathbb{S}}_{1}}\;E_{\rm sma}(\phi),\end{split}$

where

(64) $E_{\rm
sma}(\phi)=\int_{\mathbb{R}^{d}}\left[\frac{1}{2}|\nabla\phi|^{2}+V|\phi|^{2}+\beta(M)\,|\phi|^{4}\right]d\bm{x}.$

Thus $\bm{\xi}_{g}$ is the ground state of a spin-2 BEC in a spatially uniform
system, i.e. (18), and $\phi_{g}$ is the ground state of a singe-component
BEC, i.e.

(65) $\phi_{g}:=\arg\;\min_{\phi\in\widetilde{\mathbb{S}}_{1}}\;E_{\rm
sma}(\phi).$

As it has been observed numerically and proved mathematically in the
literature [7, 30], the above single mode approximation of ground state indeed
gives the ground state of the spin-2 BEC in the following cases: (i) when
$M=\pm 2$, (ii) when $M=0$, and (iii) in the ferromagnetic phase when
$\beta_{1}<0$ and $\beta_{2}>20\beta_{1}$ for $M\in[-2,2]$. Thus, in these
cases, the computation of ground state of a spin-2 BEC can be reduced to the
computation of the ground state of a single component BEC, i.e. (65).
Certainly, for all the other cases, one has to solve the original minimization
problem (6).

##  3 An efficient and accurate numerical method

In this section, we first present a normalized gradient flow (NGF) to compute
the ground state of the spin-2 BEC, then introduce three additional equations
for determining the five projection constants in the projection step for semi-
discretization of the NGF in time, and finally a full discretization of the
NGF is proposed.

###  3.1 A continuous normalized gradient flow (CNGF)

In order to compute the ground state of spin-2 BEC (6), similar as for the
single-component BEC [6, 11] and spin-1 BEC [6, 13], here we first present a
continuous normalized gradient flow (CNGF) for
$\Phi:=\Phi(\bm{x},t)=(\phi_{2},\phi_{1},\phi_{0},\phi_{-1},\phi_{-2})^{T}:=(\phi_{2}(\bm{x},t),\phi_{1}(\bm{x},t),\phi_{0}(\bm{x},t)$,
$\phi_{-1}(\bm{x},t),\phi_{-2}(\bm{x},t))^{T}$ as [43]:

(66)
$\qquad\partial_{t}\phi_{\ell}=-[H+a_{\ell}(\Phi)]\phi_{\ell}-f_{\ell}(\Phi)+[\mu_{\Phi}(t)+\ell\lambda_{\Phi}(t)]\phi_{\ell}:=\big{(}{\bf
H}\Phi\big{)}_{\ell},\ -2\leq\ell\leq 2,$

where $a_{\ell}:=a_{\ell}(\Phi)$ and $f_{\ell}:=f_{\ell}(\Phi)$
($\ell=2,\cdots,-2$) are given as

$\begin{split}a_{0}=&\beta_{0}\rho+3\beta_{1}(|\phi_{1}|^{2}+|\phi_{-1}|^{2})+0.2\beta_{2}|\phi_{0}|^{2},\\\\[5.0pt]
f_{0}=&\beta_{1}\bigg{[}\bigg{(}\frac{\sqrt{6}}{2}F_{+}-3\bar{\phi}_{1}\phi_{0}\bigg{)}\phi_{1}+\bigg{(}\frac{\sqrt{6}}{2}F_{-}-3\bar{\phi}_{-1}\phi_{0}\bigg{)}\phi_{-1}\bigg{]}+\beta_{2}\bigg{(}\frac{A_{00}}{\sqrt{5}}+0.2\beta_{2}\phi_{0}^{2}\bigg{)}\bar{\phi}_{0},\\\\[5.0pt]
a_{\pm\ell}=&\beta_{0}\rho+\beta_{1}\big{[}2\,|\phi_{\pm(3-\ell)}|^{2}\pm\ell\,F_{z}(\Phi)+(6-3\ell)\,|\phi_{0}|^{2}\big{]}+0.4\beta_{2}|\phi_{\mp\ell}|^{2},\qquad\ell=2,1\\\\[5.0pt]
f_{\pm\ell}=&\beta_{1}\bigg{\\{}\Big{[}(\ell-1)F_{\mp}+(2-\ell)F_{\pm}-2\,\bar{\phi}_{\pm(3-\ell)}\phi_{\pm\ell}\Big{]}\phi_{\pm(3-\ell)}\\\\[0.50003pt]
&+(2-\ell)\,\bigg{[}\frac{\sqrt{6}}{2}F_{\mp}-3\bar{\phi}_{0}\phi_{\pm\ell}\Big{]}\phi_{0}\bigg{\\}}+\beta_{2}\bigg{[}\frac{(-1)^{\ell}\,A_{00}}{\sqrt{5}}-0.4\phi_{\pm\ell}\phi_{\mp\ell}\bigg{]}\bar{\phi}_{\mp\ell},\quad\ell=2,1.\end{split}$

Here $\mu_{\Phi}(t)$ and $\lambda_{\Phi}(t)$ are the Lagrangian multipliers
such that both the mass (3) and magnetization (4) are conserved during the
dynamics and they can be taken as [43]:

(67)
$\begin{split}&\mu_{\Phi}(t)=\frac{\mathcal{R}(\Phi(\cdot,t))\,\mathcal{K}(\Phi(\cdot,t))-\mathcal{M}(\Phi(\cdot,t))\,\mathcal{P}(\Phi(\cdot,t))}{\mathcal{R}(\Phi(\cdot,t))\,\mathcal{N}(\Phi(\cdot,t))-\mathcal{M}^{2}(\Phi(\cdot,t))},\\\
&\lambda_{\Phi}(t)=\frac{\mathcal{N}(\Phi(\cdot,t))\,\mathcal{P}(\Phi(\cdot,t))-\mathcal{M}(\Phi(\cdot,t))\,\mathcal{K}(\Phi(\cdot,t))}{\mathcal{R}(\Phi(\cdot,t))\,\mathcal{N}(\Phi(\cdot,t))-\mathcal{M}^{2}(\Phi(\cdot,t))},\end{split}$

with $\mathcal{N}(\Phi(\cdot,t))$ and $\mathcal{M}(\Phi(\cdot,t))$ given in
(3) and (4) with $\Psi=\Phi$, respectively, and

$\begin{split}&\mathcal{R}(\Phi(\cdot,t))=\sum_{\ell=-2}^{2}\ell^{2}\|\phi_{\ell}(\cdot,t)\|^{2},\quad\mathcal{K}(\Phi(\cdot,t))=\sum_{\ell=-2}^{2}\int_{\mathbb{R}^{d}}\bar{\phi}_{\ell}(\bm{x},t)\,\big{(}{\bf
H}\Phi\big{)}_{-\ell}(\bm{x},t)\,d\bm{x},\\\
&\mathcal{P}(\Phi(\cdot,t))=\sum_{\ell=-2}^{2}\int_{\mathbb{R}^{d}}\,\ell\,\bar{\phi}_{\ell}(\bm{x},t)\,\big{(}{\bf
H}\Phi\big{)}_{-\ell}(\bm{x},t)\,d\bm{x},\qquad t\geq 0.\end{split}$

For any given initial data $\Phi(\bm{x},0):=\Phi_{0}(\bm{x})$ satisfying

(68)
$\mathcal{N}(\Phi_{0}(\bm{x}))=1,\qquad\qquad\mathcal{M}(\Phi_{0}(\bm{x}))=M,$

it is easy to show that the CNGF (66) conserves the total mass and
magnetization as well as diminishes the total energy [43], i.e.,

(69)
$\begin{split}&\mathcal{N}(\Phi(\cdot,t))\equiv\mathcal{N}(\Phi_{0})=1,\qquad\mathcal{M}(\Phi(\cdot,t))\equiv\mathcal{M}(\Phi_{0})=M,\qquad
t\geq 0,\\\
&\mathcal{E}(\Phi(\cdot,t))\leq\mathcal{E}(\Phi(\cdot,s))\leq\mathcal{E}(\Phi_{0}),\qquad{\rm
for}\quad\forall\;t\geq s\geq 0.\end{split}$

Thus the ground state of spin-2 BEC (6) can be obtained as the steady state of
the CNGF (66) with proper choice of the initial data $\Phi_{0}$ in
(LABEL:eq:CNGF-Ini).

###  3.2 A gradient flow with discrete normalization (GFDN)

Choose a time step size $\tau>0$ and denote time steps as $t_{n}=n\tau$ for
$n\geq 0$. Then a gradient flow with discrete normalization (GFDN) for
computing the ground state of the spin-2 BEC (6) can be constructed by
applying the first-order time-splitting semi-discretization of the CNGF (66)
(or as in the imaginary time method by setting $t\to it$ in the CGPE (1) and
then projecting the solution back to the manifold $\mathbb{S}$ in (7) at the
end of each time interval) as

(70)
$\partial_{t}\phi_{\ell}=-\left(H+a_{\ell}(\Phi)\right)\phi_{\ell}-f_{\ell}(\Phi),\quad
t\in[t_{n-1},t_{n}),\quad\ell=2,1,\cdots,-2,$

followed by a projection step as

(71)
$\phi_{\ell}(\bm{x},t_{n}):=\phi_{\ell}(\bm{x},t_{n}^{+})=\sigma_{\ell}^{n}\phi_{\ell}(\bm{x},t_{n}^{-}),\qquad\ell=2,1,\cdots,-2.$

Here $\phi_{\ell}(\bm{x},t_{n}^{\pm})=\lim_{t\rightarrow
t_{n}^{\pm}}\phi_{\ell}(\bm{x},t)$. Moreover, the projection constants
$\sigma_{\ell}^{n}\geq 0$ for $\ell=2,1,\cdots,-2$ are to be chosen such that

(72)
$\|\Phi(\cdot,t_{n})\|^{2}=\sum_{\ell=-2}^{2}\|\phi_{\ell}(\cdot,t_{n})\|^{2}=1,\quad\quad\sum_{\ell=-2}^{2}\ell\|\phi_{\ell}(\cdot,t_{n})\|^{2}=M.$

Plugging (71) into (72), we have

(73)
$\sum_{\ell=-2}^{2}(\sigma_{\ell}^{n})^{2}\|\phi_{\ell}(\cdot,t_{n}^{-})\|^{2}=1,\qquad\qquad\sum_{\ell=-2}^{2}\ell(\sigma_{\ell}^{n})^{2}\|\phi_{\ell}(\cdot,t_{n}^{-})\|^{2}=M.$

In fact, in the projection step, we have to determine the five projection
constants $\sigma_{\ell}^{n}$ for $\ell=2,1,\cdots,-2$ in (71). However, we
only have two equations in (73). In order to find additional proper
constraints for determining the five projection constants in the projection
step (71), we can view the GFDN (70)-(71) as a first-order time-splitting
semi-discretization of the CNGF (66). In this regard, the projection step (71)
is equivalent to solving the following nonlinear ordinary differential
equations (ODEs):

(74)
$\partial_{t}\phi_{\ell}(\bm{x},t)=[\mu_{\Phi}(t)+\ell\lambda_{\Phi}(t)]\phi_{\ell},\quad
t_{n-1}\leq t\leq t_{n},\quad\ell=2,1,\cdots,-2.$

Solving the above ODEs, one obtains for $\ell=2,1,\cdots,-2$ as

(75)
$\phi_{\ell}(\bm{x},t_{n})=\phi_{\ell}(\bm{x},t_{n-1})\exp\left(\int_{t_{n-1}}^{t_{n}}[\mu_{\Phi}(s)+\ell\lambda_{\Phi}(s)]\
ds\right):=\tilde{\sigma}_{\ell}^{n}\phi_{\ell}(\bm{x},t_{n-1}),$

which immediately suggests the following three relationships between the
constants $\tilde{\sigma}_{\ell}^{n}$ ($\ell=2,1,\cdots,-2$) as

(76)
$\tilde{\sigma}_{2}^{n}\tilde{\sigma}_{-2}^{n}=(\tilde{\sigma}_{0}^{n})^{2},\qquad\tilde{\sigma}_{1}^{n}\tilde{\sigma}_{-1}^{n}=(\tilde{\sigma}_{0}^{n})^{2},\qquad\tilde{\sigma}_{2}^{n}\tilde{\sigma}_{0}^{n}=(\tilde{\sigma}_{1}^{n})^{2}.$

Based on the above observation, we propose and adapt the following three
additional constrains for determining the five projection constants in the
project step (71) as

(77)
$\sigma_{2}^{n}\sigma_{-2}^{n}=(\sigma_{0}^{n})^{2},\qquad\quad\sigma_{1}^{n}\sigma_{-1}^{n}=(\sigma_{0}^{n})^{2},\qquad\quad\sigma_{2}^{n}\sigma_{0}^{n}=(\sigma_{1}^{n})^{2}.$

For the existence and uniqueness (in most cases) of the five project constants
$\sigma_{\ell}^{n}\geq 0$ for $\ell=2,1,\cdots,-2$ governed by (73) and (77),
we have the following result.

######  Theorem 3.1.

For sufficiently small time step size $\tau>0$ and for any $M\in[0,2)$, the
solution $\sigma_{\ell}^{n}$ ($\ell=2,\cdots,-2$) governed by (73) and (77)
can be obtained as:

(i). when $M=1$, $\|\phi_{1}(\bm{x},t_{n}^{-})\|>0$ and
$\|\phi_{\ell}(\bm{x},t_{n}^{-})\|=0\ (\ell=2,0,-1,-2)$, then

(78)
$\sigma_{1}^{n}=\frac{1}{\|\phi_{1}(\bm{x},t_{n}^{-})\|},\qquad\sigma_{0}^{n}=1,\qquad\sigma_{\ell}^{n}=(\sigma_{1}^{n})^{\ell},\quad\ell=2,-1,-2;$

(ii). when $M=0$, $\|\phi_{0}(\bm{x},t_{n}^{-})\|>0$ and
$\|\phi_{\ell}(\bm{x},t_{n}^{-})\|=0\ (\ell=2,1,-1,-2)$, then

(79)
$\sigma_{0}^{n}=\frac{1}{\|\phi_{0}(\bm{x},t_{n}^{-})\|},\qquad\sigma_{1}^{n}=1,\qquad\sigma_{\ell}^{n}=(\sigma_{0}^{n})^{1-\ell},\quad\ell=2,-1,-2;or$

(iii). for all other cases, then

(80)
$\sigma_{0}^{n}=\frac{1}{\sqrt{\sum_{\ell=-2}^{2}\lambda_{*}^{\ell}\|\phi_{\ell}(\bm{x},t_{n}^{-})\|^{2}}},\qquad\sigma_{\ell}^{n}=\sigma_{0}^{n}(\lambda_{*})^{\ell/2},\quad\ell=-2,-1,1,2,$

where $\lambda_{*}$ is the unique positive solution of the following fourth-
order algebraic equation with respect to the unknown $\lambda$ as

(81)
$\sum_{\ell=-2}^{2}(\ell-M)\|\phi_{\ell}(\bm{x},t_{n}^{-})\|^{2}\lambda^{\ell+2}=0.$

###### Proof.

Combining (73) and (77), it is straightforward to check that (78) is a
solution (not unique) of (73) and (77) in case (i), and respectively, (79) is
a solution (not unique) of (73) and (77) in case (ii).

Now we prove (80) in case (iii). Noticing (72) is also valid when $t_{n}$ is
replaced by $t_{n-1}$ for $n\geq 1$, i.e.

(82)
$\sum_{\ell=-2}^{2}\|\phi_{\ell}(\bm{x},t_{n-1})\|^{2}=1,\qquad\qquad\sum_{\ell=-2}^{2}\ell\|\phi_{\ell}(\bm{x},t_{n-1})\|^{2}=M,$

we have

(83) $\sum_{\ell=-2}^{2}(\ell-M)\|\phi_{\ell}(\bm{x},t_{n-1})\|^{2}=0.$

In addition, for sufficiently small time step size $\tau>0$,
$\|\phi_{\ell}(\cdot,t)\|\in\mathcal{C}([t_{n-1},t_{n}))$ $(\ell=2\cdots,-2)$
implies

(84)
$\|\phi_{\ell}(\bm{x},t_{n-1})\|>0\qquad\Longleftrightarrow\qquad\|\phi_{\ell}(\bm{x},t_{n}^{-})\|>0.$

Combining (84) and (82), we get
$\sum_{\ell=-2}^{2}\|\phi_{\ell}(\bm{x},t_{n}^{-})\|^{2}>0$. Therefore, case
$(iii)$ can be divided into three sub-cases:

(a) $M=1$ and
$\|\phi_{2}(\bm{x},t_{n}^{-})\|+\|\phi_{0}(\bm{x},t_{n}^{-})\|+\|\phi_{-1}(\bm{x},t_{n}^{-})\|+\|\phi_{-2}(\bm{x},t_{n}^{-})\|>0$,

(b) $M=0$ and
$\|\phi_{2}(\bm{x},t_{n}^{-})\|+\|\phi_{1}(\bm{x},t_{n}^{-})\|+\|\phi_{-1}(\bm{x},t_{n}^{-})\|+\|\phi_{-2}(\bm{x},t_{n}^{-})\|>0$,
or

(c) $M\neq 1$ and $M\neq 0$, thus $\ell-M\neq 0\ (\ell=2,\cdots,-2)$ by
noticing $M\in[0,2)$.

For all the above three sub-cases (a)-(c), noting (83)-(84), it holds

(85)
$\sum_{\ell<M}\|\phi_{\ell}(\bm{x},t_{n}^{-})\|^{2}>0,\qquad\qquad\sum_{\ell>M}\|\phi_{\ell}(\bm{x},t_{n}^{-})\|^{2}>0.$

Denote $\lambda=\sigma_{2}^{n}/\sigma_{0}^{n}$, by (73) and (77), for any
$M\in[0,2)$, we have

(86)
$\begin{split}&\sigma_{0}^{n}=\frac{1}{\sqrt{\sum_{\ell=-2}^{2}\|\phi_{\ell}(\bm{x},t_{n}^{-})\|^{2}\lambda^{\ell}}},\qquad\sigma_{j}^{n}=\sigma_{0}^{n}\lambda^{j/2},\quad
j=2,1,-1,-2,\\\
&g(\lambda):=\sum_{\ell=-2}^{2}(\ell-M)\|\phi_{\ell}(\bm{x},t_{n}^{-})\|^{2}\lambda^{\ell+2}=0.\end{split}$

Then we need only to show that $g(\lambda)$ has a unique positive root to
finish the proof.

Define $\ell_{1}^{n}$, $\ell_{2}^{n}\in\\{2,1,0,-1,-2\\}$ as

$\ell_{1}^{n}:=\min\\{\ell\;|\;\|\phi_{\ell}(\bm{x},t_{n}^{-})\|\neq
0\\}<\ell_{2}^{n}:=\max\\{\ell\;|\;\|\phi_{\ell}(\bm{x},t_{n}^{-})\|\neq
0\\}.$

Then (85) indicates that $\ell_{1}^{n}<M$ and $\ell_{2}^{n}>M$. Hence
$g(\lambda)$ can be reformulated as

$\begin{split}g(\lambda)&=\sum_{\ell_{1}^{n}\leq\ell<M}(\ell-M)\|\phi_{\ell}(\bm{x},t_{n}^{-})\|^{2}\lambda^{\ell+2}+\sum_{M<\ell\leq\ell_{2}^{n}}(\ell-M)\|\phi_{\ell}(\bm{x},t_{n}^{-})\|^{2}\lambda^{\ell+2}\\\
&=:h(\lambda)\lambda^{\ell_{1}^{n}+2},\end{split}$

where

$\begin{split}h(\lambda)&=\sum_{\ell_{1}^{n}\leq\ell<M}(\ell-M)\|\phi_{\ell}(\bm{x},t_{n}^{-})\|^{2}\lambda^{\ell-\ell_{1}^{n}}+\sum_{M<\ell\leq\ell_{2}^{n}}(\ell-M)\|\phi_{\ell}(\bm{x},t_{n}^{-})\|^{2}\lambda^{\ell-\ell_{1}^{n}}\\\
&=:h_{1}(\lambda)+h_{2}(\lambda).\end{split}$

A simple calculation shows

$\lim_{\lambda\rightarrow
0^{+}}h(\lambda)=(\ell_{1}^{n}-M)\|\phi_{\ell_{1}^{n}}(\bm{x},t_{n}^{-})\|^{2}<0,\qquad\qquad\lim_{\lambda\rightarrow+\infty}h(\lambda)=+\infty,$

which immediately implies that $h(\lambda)$ has at least one positive root
$\lambda_{*}>0$. In addition, at any positive root $\lambda=\lambda_{*}>0$ of
$h(\lambda)$, noticing (85), we have

$\begin{split}h^{\prime}(\lambda_{*})=&h_{1}^{\prime}(\lambda_{*})+h_{2}^{\prime}(\lambda_{*})\\\
=&\sum_{\ell_{1}^{n}\leq\ell<M}\left(\ell-\ell_{1}^{n}\right)(\ell-M)\|\phi_{\ell}(\bm{x},t_{n}^{-})\|^{2}(\lambda_{*})^{\ell-\ell_{1}^{n}-1}+h_{2}^{\prime}(\lambda_{*})\\\
=&\lambda_{*}^{-1}\;(M-\ell_{1}^{n})\sum_{\ell_{1}^{n}\leq\ell<M}(\ell-M)\|\phi_{\ell}(\bm{x},t_{n}^{-})\|^{2}(\lambda_{*})^{\ell-\ell_{1}^{n}}+h_{2}^{\prime}(\lambda_{*})\\\
=&\lambda_{*}^{-1}\;(M-\ell_{1}^{n})\;h_{1}(\lambda_{*})+h_{2}^{\prime}(\lambda_{*})\\\
=&\lambda_{*}^{-1}\;(M-\ell_{1}^{n})\;\big{[}h(\lambda_{n})-h_{2}(\lambda_{*})\big{]}+h_{2}^{\prime}(\lambda_{*})=\lambda_{*}^{-1}\;(\ell_{1}^{n}-M)\;h_{2}(\lambda_{*})+h_{2}^{\prime}(\lambda_{*})\\\
=&\sum_{M<\ell\leq\ell_{2}^{n}}\left(\ell_{1}^{n}-M\right)(\ell-M)\|\phi_{\ell}(\bm{x},t_{n}^{-})\|^{2}(\lambda_{*})^{\ell-\ell_{1}^{n}-1}+h_{2}^{\prime}(\lambda_{*})\\\
=&\sum_{M<\ell\leq\ell_{2}^{n}}\big{[}(\ell-M)^{2}-\left(\ell-\ell_{1}^{n}\right)(\ell-M)\big{]}\;\|\phi_{\ell}(\bm{x},t_{n}^{-})\|^{2}(\lambda_{*})^{\ell-\ell_{1}^{n}-1}+h_{2}^{\prime}(\lambda_{*})\\\
=&\sum_{M<\ell\leq\ell_{2}^{n}}(\ell-M)^{2}\|\phi_{\ell}(\bm{x},t_{n}^{-})\|^{2}(\lambda_{*})^{\ell-\ell_{1}^{n}-1}>0.\end{split}$

Therefore $h(\lambda)$ (and thus $g(\lambda)$) has exactly one positive root
$\lambda_{*}$. Substituting $\lambda=\lambda_{*}$ into (86) leads to the
formulas for the projection constants in (80). ∎

###  3.3 A backward-forward Euler finite difference discretization

Due to the trapping potential $V(\bm{x})$, the solution $\Phi(\bm{x},t)$ of
the CNGF (66) (or the GFDN (70)-(71)) decays exponentially fast to zero as
$|\bm{x}|\rightarrow\infty$. Hence, one can truncate the problem into a
bounded domain $\mathcal{D}$ with the homogeneous Dirichelet boundary
condition in practical computation. Various methods such as the backward
(-forward) Euler finite difference/sine-spectral method [9, 11, 13] can be
applied to discretize the GFDN (70)-(71). In this section, we adapt a
backward-forward Euler finite difference method (BEFD) to discretize the GFDN
(70)-(71).

To simplify the presentation, we introduce the scheme for the case of one
spatial dimension, i.e. $d=1$, defined on an interval
$\mathcal{D}=\left(a,b\right)$ with the homogeneous Dirichelet boundary
condition. Generalization to higher dimension is straightforward by tensor
product. For $d=1$, the spatial mesh size is chosen as $h=\left(b-a\right)/N$
with $N$ an even positive integer. Let

$x_{j}:=a+jh,\quad j=0,\cdots,N$

be the grid points and denote
$\Phi_{j}^{n}=(\phi_{2,j}^{n},\phi_{1,j}^{n},\phi_{0,j}^{n},\phi_{-1,j}^{n},\phi_{-2,j}^{n})^{T}$
$(n=0,1,\cdots)$ as the approximation of $\Phi(x_{j},t_{n})$. Moreover, we
denote $\Phi^{n}$ as the solution vector with component $\Phi_{j}^{n}$. Then
the GFDN (70)-(71) is discretized as

(87)
$\begin{split}&\frac{\phi^{\ast}_{\ell,j}-\phi_{\ell,j}^{n}}{\tau}=\left(\frac{1}{2}\delta_{h}^{2}-V_{j}-a_{\ell}(\Phi_{j}^{n})\right)\phi_{\ell,j}^{\ast}-f_{\ell}(\Phi_{j}^{n}),\quad
1\leq j\leq N-1,\\\
&\phi_{\ell,j}^{n+1}=\sigma_{\ell}^{n}\,\phi_{\ell,j}^{\ast},\qquad\qquad\ell=2,\cdots,-2,\qquad
j=0,1,\cdots,N.\end{split}$

Here, $V_{j}=V(x_{j})$, $\delta_{h}^{2}$ is the second-order central finite
difference operator and $\sigma_{\ell}^{n}(\ell=2,\cdots,-2)$ are the
projection constants chosen as (78)-(80). In addition, the homogeneous
Dirichelet boundary condition and initial data are discretized as

(88)
$\phi_{\ell,0}^{\ast}=\phi_{\ell,N}^{\ast}=0,\qquad\phi_{\ell,j}^{0}=\phi_{\ell}(x_{j},0),\qquad\ell=2,\cdots,-2,\quad
j=0,1,\cdots,N.$

##  4 Numerical results

In this section, we first study how to choose proper initial data for
computing numerically the ground states of spin-2 BECs, then apply the
numerical method to compute the ground states under different interaction
parameters $\beta_{1}$ and $\beta_{2}$ as well as the magnetization $M$ in one
and two dimensions. Uniqueness and non-uniqueness of the ground state are
tested and discussed based on our extensive numerical results. In our
numerical computations, the ground state $\Phi_{g}:=\lim_{n\to\infty}\Phi^{n}$
is reached numerically when
$\frac{\|\Phi^{n+1}-\Phi^{n}\|_{\infty}}{\tau}\leq\varepsilon:=10^{-7}$.

In practice, unless stated, we fix $\beta_{0}=100,$ $\tau=0.005$ and
$\mathcal{D}=[-10,10]^{d}$ for $d=1,2$. The mesh size is taken as $h_{x}=1/64$
when $d=1$, and respectively, $h_{x}=h_{y}=1/16$ when $d=2$. Moreover,
$V(\bm{x})$ is chosen either as the harmonic plus optical lattice potential

(89)
$V(\bm{x})=\sum_{j=1}^{d}\left[\frac{1}{2}\nu^{2}_{j}+\eta\,(d-1)\sin^{2}\left(q_{j}\nu_{j}\right)\right],\quad\bm{x}\in\mathcal{D},\quad
d=1,2,$

with $\nu_{1}=x$, $\nu_{2}=y$, $\eta$ and $q_{j}\;(j=1,2)$ given constants, or
the box potential

(90) $V_{\rm
box}(\bm{x})=\left\\{\begin{array}[]{ll}0,&\bm{x}\in\mathcal{D},\\\\[2.5pt]
+\infty,&\hbox{\rm otherwise}.\end{array}\right.$

###  4.1 Choice of initial data and uniqueness of the ground state

A proper choice of the initial data $\Phi_{0}(\bm{x})$ usually improves
significantly the efficiency and accuracy of the GFDN (70)-(71). For cases
where SMA is valid, e.g. the nematic and/or cyclic phase with $M=0$ and the
ferromagnetic phase [7], see also Fig. 5), one can either simply construct the
ground state via (65) by solving the ground state of the single component BEC
(65), or directly solve the ground state via the GFDN (87)-(71) with initial
data (it is a reasonable and natural choice by noticing section 2.3) taken as

(91) $\bm{\Phi}_{0}(\bm{x})=\bm{\xi}_{g}\,\phi(\bm{x}),$

where $\bm{\xi}_{g}$ is given in Theorem 2.1 and $\phi(\bm{x})$ is an
approximation of the ground state of the single component BEC (65), e.g., the
harmonic oscillator approximation $\phi_{g}^{\rm hos}(\bm{x})$ for small
$\beta(M)$ and/or the Thomas-Fermi approximation $\phi_{g}^{\rm TF}(\bm{x})$
for large $\beta(M)$ [11, 13]. For other cases where the SMA is invalid, the
initial data $\Phi_{0}$ can be chosen either as (91) or as a more general
initial set-up

(92) $\bm{\Phi}^{0}(\bm{x})=\phi(\bm{x})\bm{\xi},$

where

(93)
$\bm{\xi}=\frac{1}{2}\Big{(}\sqrt{2+M-2\sigma},\,\sqrt{\sigma},\,\sqrt{2\sigma},\,\sqrt{\sigma},\,\sqrt{2-M-2\sigma}\Big{)}^{T},\quad\sigma\in[0,1-M/2].$

Extensive numerical comparison (not shown here for brevity) for different
initial data show that the GFDN would usually converge faster with initial
data (91) than with other types of initial data. Based on those comparison, we
would conclude and suggest the choice of initial data as follows:

a). For ferromagnetic phase, $\bm{\xi}_{g}$ is suggested to be chosen as (48)
for $\forall M\in[0,2)$. Meanwhile, the ground state is found to be unique.

b). For nematic phase, if $M\in(0,2)$, $\bm{\xi}_{g}$ is suggested to be
chosen as (49), the ground state is found to be unique. However, if $M=0$, the
ground state is not unique. Hence, any initial data chosen as (49) works and
probably converge to different ground states.

c). For cyclic phase, if $M\in(0,2)$, then $\bm{\xi}_{g}$ is suggested to be
chosen as (52) with $\theta=\arctan\sqrt{(2-M)/(1+M)}$, i.e.,

(94) $\bm{\xi}_{g}=\left(\sqrt{(M+1)/3},\ 0,\ 0,\ \sqrt{(2-M)/3},\
0\right)^{T}.$

The ground state is found to be unique, which is essentially different from
the spatially uniform system where the ground state is not unique. While if
$M=0$, similar as the nematic phase, the ground state is not unique, thus any
initial data works and probably converge to different ground states.

Fig. 2: Time evolution of the total energy $\mathcal{E}(t)$ in Example 4.1
with different time step size $\tau$.

######  Example 4.1.

Here we show the energy-diminishing property of our numerical method. To this
end, we let $M=0.5$, $\beta_{1}=1$ and $\beta_{2}=-2$ (i.e. nematic phase).
The ground state of a spin-2 BEC (1) is computed by the BEFD (87) with time
step $\tau=0.1/0.01/0.005$ and the initial data (92)-(93) with $\sigma=0$,
i.e., $\bm{\xi}$ is chosen as the ground state in a spatially uniform system
(49). Fig. 2 shows the evolution of the energy
$\mathcal{E}(t):=\mathcal{E}(\Phi(\cdot,t))$ with different time step $\tau$.

From Fig. 2 and additional experiments not shown here for brevity, we can see
that: (i). the energy is diminishing for different time step size $\tau$, even
for the relatively large step size $\tau=0.1$ (cf. Fig. 2), and (ii). the GFDN
with different initial data converge to the same ground state. In addition,
when $\sigma=0$, i.e., $\bm{\xi}$ is chosen as the ground state in a spatially
uniform system (49), the GFDN usually converges in the fastest way.

Fig. 3: Plots of the wave function of the ground states $\phi_{\ell}$
($\ell=2,1,0,-1,-2$) in Cases 1-4 (from top to bottom) that computed by
different initial data in the Example 4.2.

###  4.2 Applications

In this subsection, we apply our numerical method to compute the ground state
of spin-2 BEC with different parameter regimes.

######  Example 4.2.

Here we further study the non-uniqueness of the ground state under some
parameter regimes. To this end, we take $d=1$, and carry out the following
cases:

(i) Case 1. $\beta_{1}=-1$, $\beta_{2}=-20$, $M=0.5$, $\bm{\xi}_{g}$ is taken
as (48)) and (49), respectively;

(ii) Case 2. $\beta_{1}=0$, $\beta_{2}=1$, $M=0.8$, $\bm{\xi}_{g}$ is taken as
(48) and (94), respectively;

(iii) Case 3. $\beta_{1}=10$, $\beta_{2}=2$, $M=0$, $\bm{\xi}_{g}$ is taken as
(51) and (52) with $\theta=\arcsin(1/5)$, respectively;

(iv) Case 4. $\beta_{1}=1$, $\beta_{2}=-2$, $M=0$, $\bm{\xi}_{g}$ is taken as
(50) with $(\gamma_{1}=0.4,\theta=\arcsin\sqrt{2/15})$ and
$\theta=\arcsin(2\sqrt{2}/5)$, respectively.

Fig. 3 depicts the wave functions of different ground states computed by
different initial data in Cases 1-4. For each case, different initial data
converges to different ground states with the same energy. From this figure
and additional results not shown here for brevity, we can see that the ground
states are not unique for the following four cases (cf. Fig. 3): (a).
$\forall\,M\in[0,2)$ and $\beta_{1}<0$, $\beta_{2}=20\beta_{1}$. (b).
$\forall\,M\in[0,2)$ and $\beta_{2}>0$, $\beta_{1}=0$. (c). $M=0$ and
$\beta_{1}>0$, $\beta_{2}>0$. (d). $M=0$ and $\beta_{2}<0$,
$\beta_{2}<20\beta_{1}$.

Fig. 4: Wave functions of ground states, i.e., $\phi_{\ell}^{g}$
($\ell=2,1,0,-1,-2$) with different magnetizations $M=0,\;0.5,\;1.5$ (left to
right) for Cases 5-7 (top to bottom) in the Example 4.3.

Fig. 5: Plots of $\phi^{g}_{1}/\phi^{g}_{2}$ (dashed line),
$\phi^{g}_{0}/\phi^{g}_{2}$ (solid line), $\phi^{g}_{-1}/\phi^{g}_{2}$ (dotted
line), and $\phi^{g}_{-2}/\phi^{g}_{2}$ (dashed-dotted line) for Cases 5-7
(top to bottom) in the Example 4.3, respectively, to analyze the SMA property
for different parameters $(\beta_{1},\beta_{2},M)$.

Fig. 6: The component masses $\mathcal{N}_{\ell}$ ($\ell=2,1,0,-1,-2$) of the
ground states for different $\beta_{1}$ & $\beta_{2}$ with $M=0$ (upper) and
$M=0.5$ (lower) in Case 8 in the Example 4.3.

######  Example 4.3.

We study the wave functions and SMA property of the ground states in different
parameter regimes. To this end, we take $d=1$ and consider following three
cases:

(i) Case 5. ferromagnetic phase, we take $\beta_{1}=-1$ and $\beta_{2}=2$;

(ii) Case 6. nematic phase, we choose $\beta_{1}=1$ and $\beta_{2}=-2$; and

(iii) Case 7. cyclic phase, we let $\beta_{1}=10$ and $\beta_{2}=2$.

(iv) Case 8. let $M=0$ or $0.5$, take initial as (91) with $\bm{\xi}_{g}$
reading as (48), (49) and (94), vary both $\beta_{1}$ and $\beta_{2}$.

Fig. 4 depicts the wave functions of the ground states in Cases 5-7 for
different magnetizations $M$. Fig. 5 shows the SMA property for different
parameters $(\beta_{1},\beta_{2},M)$, while Fig. 6 shows the changes of mass
$\mathcal{N}_{\ell}$ ($\ell=2,1,0,-1,-2$) with respect to parameters
$\beta_{1}$ & $\beta_{2}$ for fixed $M=0$ and $M=0.5$, respectively.

From Figs. 4-6 and extensive numerical experiments not shown here for brevity,
we observe that: (i). When $M\in(0,2)$, for ferromagnetic phase,
$\phi_{\ell}^{g}>0$ for all $\ell=-2,-1,0,1,2$ (cf. Fig. 4 (top row)). For
nematic phase, $\phi_{2}^{g}>0$ & $\phi_{-2}^{g}>0$, while
$\phi_{1}^{g}=\phi_{-1}^{g}=\phi_{0}^{g}\equiv 0$ (cf. Fig. 4 (middle row)).
For cyclic phase, $\phi_{2}^{g}>0$ & $\phi_{-1}^{g}>0$, while
$\phi_{1}^{g}=\phi_{-2}^{g}=\phi_{0}^{g}\equiv 0$ (cf. Fig. 4( bottom row)).
(ii). The SMA is valid for the ground states with ferromagnetic phase and
those with the nematic or cyclic phases as well as zero magnetization,
however, it is invalid for the other cases (cf. Fig. 5). (iii). For fixed $M$,
when vary parameters $(\beta_{1},\beta_{2})$ within the same parameters regime
(i.e., the regime of ferromagnetic, nematic and cyclic phases), the mass
$\mathcal{N}_{\ell}$ $(\ell=2,1,0,-1,-2)$ are independent of
$(\beta_{1},\beta_{2})$ (cf. Fig. 6).

Fig. 7: $M=0.5$. Contour plots of the components of the ground states
$\phi_{\ell}^{g}$ (from left to right, $\ell=2,1,0,-1,-2$) in Cases 9-11 (top
to bottom) with the harmonic potential in the Example 4.4.  Fig. 8: $M=0.5$.
Contour plots of the components of the ground states $\phi_{\ell}^{g}$ (from
left to right, $\ell=2,1,0,-1,-2$) in Cases 9-11 (top to bottom) with the box
potential in the Example 4.4.  Fig. 9: $M=0.5$. Contour plots of the
components of the ground states $\phi_{\ell}^{g}$ (from left to right,
$\ell=2,1,0,-1,-2$) in Cases 9-11 (top to bottom) with the optical lattice
potential in the Example 4.4.

######  Example 4.4.

Here, we study the ground state of a two-dimensional spin-2 BEC with
harmonic/box/ optical lattice potentials. To this end, we take $d=2$,
$q_{1}=q_{2}=\pi/2$, $\eta=0$/$\eta=10$ in (89) for the harmonic/optical
lattice potential, and choose a box potential $V_{box}(\bm{x})$ as in (90). We
consider the following three cases:

(i) Case 9. ferromagnetic phase, let $\beta_{1}=-1$ & $\beta_{2}=-5$;

(ii) Case 10. nematic phase, choose $\beta_{1}=-1$ & $\beta_{2}=-25$;

(iii) Case 11. cyclic phase, take $\beta_{1}=10$ & $\beta_{2}=2$.

Figs. 7-9 show the plots of the wave functions of the ground states in Cases
9-11 with the harmonic, box and optical lattice potential, respectively. From
these results and additional numerical experiments not shown here for brevity,
one finds that our method can be applied to compute the ground state of spin-2
BEC with general potentials. Additionally, similar as the 1-d case, the
uniqueness, validity of SMA and phenomena of vanishing-component of the ground
state can also be concluded in the two-dimensional case.

##  5 Conclusion

We proposed an efficient and accurate normalized gradient flow method for
computing the ground states of spin-2 BEC by introducing three additional
projection constraints, in addition to the conservation of the total mass and
magnetization. A backward-forward finite difference method was applied to
fully discretize the gradient flow with discrete normalization. Moreover, the
ground states in spatially uniform system, i.e $V(\bm{x})=0$, were solved
analytically, which give important hints to understand the properties of
ground states in spatially non-uniform systems as well as to the choice of
initial data for numerical calculations. The numerical method was then applied
to study the properties of the ground state of spin-2 BECs with the
ferromagnetic/nematic/cyclic phases and harmonic/box/optical lattice trapping
potentials in one/two dimensions. Various numerical experiments on the ground
states of spin-2 BEC with different parameters and potentials were carried
out, which suggest some interesting properties about the ground states. For
example, the parameter regimes for the uniqueness, validity of SMA and
phenomena of vanishing-component were numerically partially found, and the
property that the mass of each component of the ground state in the same phase
is independent of the interaction parameters $(\beta_{1},\beta_{2})$ was
observed. Rigorous mathematical justifications for these observations are on-
going.

## References

  * [1] A. Aftalion and Q. Du, Vortices in a rotating Bose-Einstein condensate: Critical angular velocities and energy diagrams in the Thomas-Fermi regime, Phys. Rev. A, 64 (2001), 063603.
  * [2] M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman and E. A. Cornell, Observation of Bose-Einstein condensation in a dilute atomic vapor, Science, 269 (1995), 198–201.
  * [3] X. Antoine and R. Duboscq, Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose-Einstein condensates, J. Comput. Phys., 258 (2014), 509–523.
  * [4] X. Antoine, A. Levitt and Q. Tang, Efficient spectral computation of the stationary states of rotating Bose-Einstein condensates by preconditioned nonlinear conjugate gradient methods, J. Comput. Phys., 343 (2017), 92–109.
  * [5] W. Bao, Ground states and dynamics of multicomponent Bose-Einstein condensates, Multiscale Model. Simul., 2 (2004), 210–236.
  * [6] W. Bao and Y. Cai, Mathematical theory and numerical methods for Bose-Einstein condensation, Kinet. Relat. Mod., 6 (2013), 1-135.
  * [7] W. Bao and Y. Cai, Mathematical models and numerical methods for spinor Bose-Einstein condensates, Commun. Comput. Phys., 24 (2018), 899–965.
  * [8] W. Bao and M. Chai, A uniformly convergent numerical method for singularly perturbed nonlinear eigenvalue problems, Commun. Comput. Phys., 4 (2008), 135–160.
  * [9] W. Bao, I-L Chern and F. Lim, Efficient and spectrally accurate numerical methods for computing ground and first excited states in Bose-Einstein condensates, J. Comput. Phys., 219 (2006), 836–854.
  * [10] W. Bao, I-L. Chern and Y. Zhang, Efficient methods for computing ground states of spin-1 Bose-Einstein condensates based on their characterizations, J. Comput. Phys., 253 (2013), 189–208.
  * [11] W. Bao and Q. Du, Computing the ground state solution of Bose-Einstein condensates by a normalized gradient flow, SIAM J. Sci. Comput., 25 (2004), 1674–1697.
  * [12] W. Bao, D. Jaksch and P. A. Markowich, Numerical solution of the Gross-Pitaevskii equation for Bose-Einstein condensation, J. Comput. Phys., 187 (2003), 318–342.
  * [13] W. Bao and F. Lim, Computing ground states of spin-1 Bose-Einstein condensates by the normalized gradient flow, SIAM J. Sci. Comput., 30 (2008), 1925–1948.
  * [14] W. Bao and W. Tang, Ground state solution of trapped interacting Bose-Einstein condensate by directly minimizing the energy functional, J. Comput. Phys., 187 (2003), 230–254.
  * [15] W. Bao and H. Wang, A mass and magnetization conservative and energy-diminishing numerical method for computing ground state of spin-1 Bose-Einstein condensates, SIAM J. Numer. Anal., 45 (2007), 2177–2200.
  * [16] R. Barnett, S. Mukerjee and J. E. Moore, Vortex lattice transitions in cyclic spinor condensates, Phys. Rev. Lett., 100 (2008), 240405.
  * [17] C. C. Bradley, C. A. Sackett, J. J. Tollett and R. G. Hulet, Evidence of Bose-Einstein condensation in an atomic gas with attractive interaction, Phys. Rev. Lett., 75 (1995), 1687–1690.
  * [18] Y. Cai and W. Liu, Efficient and accurate gradient flow methods for computing ground states of spinor Bose-Einstein condensates, J. Comput. Phys., 433 (2021), 110183.
  * [19] S.-L. Chang, C.-S. Chien and B. Jeng, Liapunov-Schmidt reduction and continuation for nonlinear Schrödinger equations, SIAM J. Sci. Comput., 29 (2007), 729–755.
  * [20] S.-L. Chang, C.-S. Chien and B. Jeng, Computing wave functions of nonlinear Schrödinger equations: A time independent approach, J. Comput. Phys., 226 (2007) 104–130.
  * [21] M. Chang, C. Hamley, M. D. Barrett, J. A. Sauer, K. M. Fortier, W. Zhang, L. You and M. S. Chapman, Observation of spinor dynamics in optically trapped 87Rb Bose-Einstein condensates, Phys. Rev. Lett., 92 (2004), 140403.
  * [22] S. Chang, W. Lin and S. F. Shieh, Gauss-Seidel-type methods for energy states of a multi-component Bose-Einstein condensate, J. Comput. Phys., 202 (2005), 367–390.
  * [23] M. L. Chiofalo, S. Succi and M. P. Tosi, Ground state of trapped interacting Bose-Einstein condensates by an explicit imaginary-time algorithm, Phys. Rev. E, 62 (2000), 007438.
  * [24] C. V. Ciobanu, S. K. Yip and T. L. Ho, Phase diagrams of F=2 spinor Bose condensates, Phys. Rev. A, 61 (2000), 033607.
  * [25] F. Dalfovo, S. Giorgini, L. P. Pitaevskii and S. Stringari, Theory of Bose-Einstein condensation in trapped gases, Rev. Mod. Phys., 71 (1999), 463–512.
  * [26] I. Danaila and P. Kazemi, A new Sobolev gradient method for direct minimization of the Gross-Pitaevskii energy with rotation, SIAM J. Sci. Comput., 32 (2010), 2447–2467.
  * [27] K. B. Davis, M. O. Mewes, M. R. Andrews, N. J. van Druten, D. S. Durfee, D. M. Kurn and W. Ketterle, Bose-Einstein condensation in a gas of sodium atoms, Phys. Rev. Lett., 75 (1995), 3969–3973.
  * [28] A. Görlitz, T. L. Gustavson, A. E. Leanhardt, R. Löw, A. P. Chikkatur, S. Gupta, S. Inouye, D. E. Pritchard and W. Ketterle, Sodium Bose-Einstein condensates in the F=2 state in a large-volume optical trap, Phys. Rev. Lett., 90 (2003), 090401.
  * [29] T. Isoshima, M. Okano, H. Yasuda, K. Kasa, J. A. M. Huhtamäki, M. Kumakura and Y. Takahashi, Spontaneous splitting of a quadruply charged vortex, Phys. Rev. Lett., 99 (2007), 200403.
  * [30] Y. Kawaguchi and M. Ueda, Spinor Bose-Einstein condensates, Phys. Rep., 520 (2012), 253–381.
  * [31] M. Koashi and M. Ueda, Exact eigenstates and magnetic response of spin-1 and spin-2 Bose-Einstein condensates, Phys. Rev. Lett., 84 (2000), 1066.
  * [32] T. Kuwamoto, K. Araki, T. Eno and T. Hirano, Magnetic field dependence of the dynamics of 87Rb spin-2 Bose-Einstein condensates, Phys. Rev. A, 69 (2004), 063604.
  * [33] T. Kuwamoto, H. Usuda, S. Tojo and T. Hirano, Dynamics of quadruply quantized vortices in 87Rb Bose-Einstein condensates confined in magnetic and optical traps, J. Phys. Soc. Japan, 79 (2010), 034004.
  * [34] L. Lin and I-Liang Chern, A kinetic energy reduction technique and characterizations of the ground states of spin-1 Bose-Einstein condensates, Discrete Cont. Dyn-B, 19 (2014), 1119-1128.
  * [35] W. Liu and Y. Cai, Normalized gradient flow with Lagrange multiplier for computing ground states of Bose–Einstein condensates, SIAM J. Sci. Comput., 43 (2021), B219–B242.
  * [36] B. Pasquiou, E. Maréchal, G. Bismut, P. Pedri, L. Vernac, O. Gorceix, B. Laburthe-Tolra, Spontaneous demagnetization of a dipolar spinor Bose gas in an ultralow magnetic field, Phys. Rev. Lett., 106 (2011), 255303.
  * [37] L. P. Pitaevskii and S. Stringari, Bose-Einstein Condensation, Clarendon Press, Oxford, 2003.
  * [38] H. Schmaljohann, M. Erhard, J. Kronjäger, M. Kottke, S. van Staa, L. Cacciapuoti, J. J. Arlt, K. Bongs and K. Sengstock, Dynamics of F=2 spinor Bose-Einstein condensates, Phys. Rev. Lett., 92 (2004), 040402.
  * [39] L. Simon, Asymptotics for a class of nonlinear evolution equations, with applications to geometric problems, Ann. of Math., 118 (1983), 525–571.
  * [40] D. M. Stamper-Kurn, M. R. Andrews, A. P. Chikkatur, S. Inouye, H. J. Miesner, J. Stenger and W. Ketterle, Optical confinement of a Bose-Einstein condensate, Phys. Rev. Lett., 80 (1998), 2027–2030.
  * [41] T. Tian, Y. Cai, X. Wu and Z. Wen, Ground states of spin-F Bose-Einstein condensates, SIAM J. Sci. Comput., 42 (2020), B983–B1013.
  * [42] H. Wang, An efficient numerical method for computing dynamics of spin F=2 Bose-Einstein condensates, J. Comput. Phys., 230 (2011), 6155–6168.
  * [43] H. Wang, A projection gradient method for computing ground state of spin-2 Bose-Einstein condensates, J. Comput. Phys., 274 (2014), 473–488.
  * [44] A. Widera, F. Gerbier, S. Folling, T. Gericke, O. Mandel and I. Bloch, Precision measurement of spin-dependent interaction strengths for spin-1 and spin-2 87Rb atoms, New J. Phy., 8 (2006), 152.
  * [45] X. Wu, Z. Wen and W. Bao, A regularized Newton method for computing ground states of Bose-Einstein condensates, J. Sci. Comput. 73 (2017), 303–329.
  * [46] L. Xia, X. Xu, F. Yang, W. Xiong, J. Li, Q. Ma, X. Zhou, H. Guo and X. Chen, Manipulation of the quantum state by Majorana transition in spinor Bose-Einstein condensates, Phys. Rev. A, 77 (2008), 043622.
  * [47] Z. Yu and Z. Jiao, Tunnelling effect of the spin-2 Bose condensate driven by external magnetic fields, Phys. Lett. A, 320 (2004), 318–321.
  * [48] Y. Yuan, Z. Xu, Q. Tang and H. Wang, The numerical study of the ground states of spin-1 Bose-Einstein condensates with spin-orbit-coupling, E. Asian J. Appl. Math., 8 (2018), 598–610.
  * [49] W. Zhang and L. You, An effective quasi-one-dimensional description of a spin-1 atomic condensate, Phys. Rev. A, 71 (2005), 025603.
  * [50] W. Zhang, S. Yi and L. You, Mean field ground state of a spin-1 condensate in a magnetic field, New J. Phys., 5 (2003), 77–89.