text:
# Theory of spin-2 Bose-Einstein condensates: spin correlations, magnetic response, and excitation spectra

## I Introduction

Bose-Einstein condensates (BECs) of trapped alkali atoms have internal degrees of freedom due to the hyperfine spin of the atoms. When a BEC is trapped in a magnetic potential, these degrees of freedom are frozen and the state of the BEC is described at a mean-field level by a scalar order parameter. When a BEC is trapped in an optical potential, however, the spin degrees of freedom are liberated, giving rise to a rich variety of phenomena such as spin domains  and textures . Here the order parameter has $`2f+1`$ components that transform under spatial rotation as the spherical tensor of rank $`f`$, where $`\mathrm{}f`$ is the hyperfine spin of bosons.

Mean field theories (MFTs) of spinor BECs were put forth for both spin-1  and spin-2  cases. According to them, the $`m=0`$ magnetic sublevel of an antiferromagnetic BEC is not populated in the presence of magnetic field for both spin-1 and spin-2 cases. However, Law et al.  used many-body theory to show that the $`m=0`$ sublevel of a spin-1 BEC is macroscopically populated due to the formation of spin-singlet “pairs” of bosons. It was subsequently shown  that the $`m=0`$ sublevel of a spin-2 BEC is also macroscopically occupied due to the formation of spin-singlet “trios” of bosons. The physics common to both cases is that the spin-singlet state is isotropic and therefore each magnetic sublevel shares the equal population.

A question then arises as to what extent and under what conditions MFTs are applicable. It is now understood  that the validity of MFTs is quickly restored with increasing an applied magnetic field. Thus for the many-body spin correlations to manifest themselves, the external magnetic field have to be very low.

The spin-singlet pairs of bosons should be distinguished from Cooper pairs of electrons or those of 3He due to the statistical difference of constituent particles. The Cooper pairs consist of fermions, so that the state is symmetric only under the permutations that do not break any pair. For the case of spin-singlet pairs of bosons, the state is symmetric under any permutation of the constituents. The Bose-Einstein statistics leads to a constructive interference among permuted terms, giving rise to a highly nonlinear magnetic response to be discussed later. In contrast with usual antiferromagnets, where antiparallel spins are alternately aligned (Neel order), “antiferromagnetic” BECs do not possess such a long-range spatial order because the system lacks crystal order. The antiferromagnetic phase of BECs is also called polar .

In Refs. , only spin degrees of freedom are considered by assuming that the spatial degrees of freedom are frozen. In this paper, we relax this restriction and develop a theory of spin-2 BECs that enables us to study many-body ground states and the excitation spectrum thereof on an equal footing. For spin-1 BECs, this program has been carried out in Ref. . Many-body spin correlations and magnetic response of BECs, the results of which were briefly reported in Ref. , are expounded. The role of symmetry of the ground state in determining the character of the excitation spectrum is also elucidated.

This paper is organized as follows. Section II derives an effective Hamiltonian that enables us to study many-body spin correlations and low-lying excitation spectrum of spin-2 BECs on an equal footing. Section III reviews mean-field properties of each phase of BECs. Section IV studies many-body spin correlations and magnetic response of BEC. The energy eigenstate is explicitly constructed using the creation operators of boson pairs and trios. The degeneracy of the eigenstate is examined and some novel magnetic response such as a huge jump in magnetization and the robustness of the minimum-magnetization state against an applied magnetic field are discussed. Section V derives excitation spectra of BECs using the Bogoliubov approximation. All excitation spectra are obtained analytically and the relations of their characters to the symmetry of the ground state are discussed. Section VI summarizes the main results of the present paper. Appendix A recapitulates the parametrization of the order parameter of spin-2 BECs, and appendix B describes a method of calculating Zeeman-level populations.

## II Formulation of the problem

### A Interaction Hamiltonian

Consider a system of identical bosons with hyperfine spin $`f`$ and let $`\widehat{\mathrm{\Psi }}_m(𝒓)`$ ($`m=f,f1,\mathrm{},f`$) be the field operator that annihilates at position $`𝒓`$ a boson with magnetic quantum number $`m`$. The field operators are assumed to obey the canonical commutation relations

$`[\widehat{\mathrm{\Psi }}_m(𝒓),\widehat{\mathrm{\Psi }}_n^{}(𝒓^{})]=\delta _{mn}\delta (𝒓𝒓^{}),[\widehat{\mathrm{\Psi }}_m(𝒓),\widehat{\mathrm{\Psi }}_n(𝒓^{})]=0,[\widehat{\mathrm{\Psi }}_m^{}(𝒓),\widehat{\mathrm{\Psi }}_n^{}(𝒓^{})]=0,`$ (1)
where the Kronecker’s delta $`\delta _{mn}`$ takes on the value of 1 if $`m=n`$ and 0 otherwise. The Bose-Einstein statistics requires that the total spin of any two bosons whose relative orbital angular momentum is zero be restricted to $`F=2f,2f2,\mathrm{},0`$. We may therefore use $`F`$ as an index for classifying binary interactions between identical bosons:

$`\widehat{V}={\displaystyle \underset{F=0,2,\mathrm{},2f}{}}\widehat{V}^{(F)},`$ (2)
where $`\widehat{V}^{(F)}`$ describes an interaction between two bosons whose total spin is $`F`$. To construct $`\widehat{V}^{(F)}`$, consider the operator $`\widehat{A}_{FM}(𝒓,𝒓^{})`$ that annihilates at positions $`𝒓`$ and $`𝒓^{}`$ two bosons with total spin $`F`$ and total magnetic quantum number $`M`$:

$`\widehat{A}_{FM}(𝒓,𝒓^{})={\displaystyle \underset{m_1,m_2=f}{\overset{f}{}}}F,M|f,m_1;f,m_2\widehat{\mathrm{\Psi }}_{m_1}(𝒓)\widehat{\mathrm{\Psi }}_{m_2}(𝒓^{}),`$ (3)
where $`F,M|f,m_1;f,m_2`$ is the Clebsch-Gordan coefficient. We may use $`\widehat{A}_{FM}`$ to construct $`\widehat{V}^{(F)}`$ as

$`\widehat{V}^{(F)}={\displaystyle \frac{1}{2}}{\displaystyle 𝑑𝒓𝑑𝒓^{}v^{(F)}(𝒓,𝒓^{})\underset{M=F}{\overset{F}{}}\widehat{A}_{FM}^{}(𝒓,𝒓^{})\widehat{A}_{FM}(𝒓,𝒓^{})},`$ (4)
where $`v^{(F)}(𝒓,𝒓^{})`$ describes the dependence of the interaction on the positions of the particles. Because of the completeness relation $`_{F,M}|F,MF,M|=\widehat{1}`$, where $`\widehat{1}`$ is the identity operator, we find that

$`{\displaystyle \underset{F=0,2,\mathrm{},2f}{}}{\displaystyle \underset{M=F}{\overset{F}{}}}\widehat{A}_{FM}^{}(𝒓,𝒓^{})\widehat{A}_{FM}(𝒓,𝒓^{})`$ $`=`$ $`:\widehat{n}(𝒓)\widehat{n}(𝒓^{}):,`$ (5)
where

$`\widehat{n}(𝒓){\displaystyle \underset{m=f}{\overset{f}{}}}\widehat{\mathrm{\Psi }}_m^{}(𝒓)\widehat{\mathrm{\Psi }}_m(𝒓)`$ (6)
is the total density operator and :: denotes normal ordering, that is, annihilation operators are placed to the right of creation operators. Integrating Eq. (5) over $`𝒓,𝒓^{}`$ yields

$`{\displaystyle 𝑑𝒓𝑑𝒓^{}\underset{F,M}{}\widehat{A}_{FM}^{}(𝒓,𝒓^{})\widehat{A}_{FM}(𝒓,𝐫^{})}=\widehat{N}(\widehat{N}1),`$ (7)
where $`\widehat{N}\widehat{n}(𝒓)𝑑𝒓`$ is the total number of bosons.

In the case of a dilute Bose-Einstein condensate of neutral atoms, we may to a good approximation assume that $`v^{(F)}(𝒓,𝒓^{})=g_F\delta (𝒓𝒓^{})`$, where $`g_F`$ characterizes the strength of the interaction between two bosons whose total spin is $`F`$, and is related to the corresponding s-wave scattering length $`a_F`$ as

$`g_F={\displaystyle \frac{4\pi \mathrm{}^2}{M}}a_F.`$ (8)
Equation (4) then becomes

$`\widehat{V}^{(F)}={\displaystyle \frac{g_F}{2}}{\displaystyle 𝑑𝒓\underset{M=F}{\overset{F}{}}\widehat{A}_{FM}^{}(𝒓,𝒓)\widehat{A}_{FM}(𝒓,𝒓)}.`$ (9)
In the following discussions we shall focus on this case and therefore denote $`\widehat{A}_{FM}(𝒓,𝒓)`$ simply as $`\widehat{A}_{FM}(𝒓)`$.

When $`f=2`$, $`F`$ can take on values 0, 2, and 4. For $`F=0`$, we have

$`\widehat{V}^{(0)}={\displaystyle \frac{g_0}{2}}{\displaystyle 𝑑𝒓\widehat{A}_{00}^{}(𝒓)\widehat{A}_{00}(𝒓)},`$ (10)
where

$`\widehat{A}_{00}(𝒓)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{5}}}\left[2\widehat{\mathrm{\Psi }}_2(𝒓)\widehat{\mathrm{\Psi }}_2(𝒓)2\widehat{\mathrm{\Psi }}_1(𝒓)\widehat{\mathrm{\Psi }}_1(𝒓)+\widehat{\mathrm{\Psi }}_0^2(𝒓)\right].`$ (11)
For $`F=2`$, we have

$`\widehat{V}^{(2)}={\displaystyle \frac{g_2}{2}}{\displaystyle 𝑑𝒓\underset{M=2}{\overset{2}{}}\widehat{A}_{2M}^{}(𝒓)\widehat{A}_{2M}(𝒓)}.`$ (12)
For $`F=4`$ we have

$`\widehat{V}^{(4)}`$ $`=`$ $`{\displaystyle \frac{g_4}{2}}{\displaystyle 𝑑𝒓\underset{M=4}{\overset{4}{}}\widehat{A}_{4M}^{}(𝒓)\widehat{A}_{4M}(𝒓)}`$ (13)
$`=`$ $`{\displaystyle \frac{g_4}{2}}{\displaystyle }d𝒓[:\widehat{n}^2(𝒓):\widehat{A}_{00}^{}(𝒓)\widehat{A}_{00}(𝒓){\displaystyle \underset{M=2}{\overset{2}{}}}\widehat{A}_{2M}^{}(𝒓)\widehat{A}_{2M}(𝒓)],`$ (14)
where Eq. (5) was used in obtaining the second equality. Summing Eqs. (10), (12) and (14), we obtain the interaction Hamiltonian as

$`\widehat{V}`$ $`=`$ $`\widehat{V}^{(0)}+\widehat{V}^{(2)}+\widehat{V}^{(4)}`$ (15)
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle }d𝒓[g_4:\widehat{n}^2(𝒓):+(g_0g_4)\widehat{A}_{00}^{}(𝒓)\widehat{A}_{00}(𝒓)+(g_2g_4){\displaystyle \underset{M=2}{\overset{2}{}}}\widehat{A}_{2M}^{}(𝒓)\widehat{A}_{2M}(𝒓)].`$ (16)
To eliminate the last term in Eq. (16), we note the following operator identity:

$`{\displaystyle \frac{1}{7}}:\widehat{𝑭}^2(𝒓):+{\displaystyle \underset{M=2}{\overset{2}{}}}\widehat{A}_{2M}^{}(𝒓)\widehat{A}_{2M}(𝒓)+{\displaystyle \frac{10}{7}}\widehat{A}_{00}^{}(𝒓)\widehat{A}_{00}(𝒓)={\displaystyle \frac{4}{7}}:\widehat{n}^2(𝒓):,`$ (17)
where $`\widehat{𝑭}=(\widehat{F}^x,\widehat{F}^y,\widehat{F}^z)`$ represents the spin density operators defined by

$`\widehat{F}_i(𝒓)={\displaystyle \underset{m,n=1}{\overset{1}{}}}f_{mn}^i\widehat{\mathrm{\Psi }}_m^{}(𝒓)\widehat{\mathrm{\Psi }}_n(𝒓)(i=x,y,z)`$ (18)
with $`f_{mn}^i`$ $`(i=x,y,z)`$ being the $`(m,n)`$-components of spin-2 matrices $`f^i`$ given by

$`f^x={\displaystyle \frac{1}{2}}\left(\begin{array}{ccccc}0& 2& 0& 0& 0\\ 2& 0& \sqrt{6}& 0& 0\\ 0& \sqrt{6}& 0& \sqrt{6}& 0\\ 0& 0& \sqrt{6}& 0& 2\\ 0& 0& 0& 2& 0\end{array}\right),f^y={\displaystyle \frac{i}{2}}\left(\begin{array}{ccccc}0& 2& 0& 0& 0\\ 2& 0& \sqrt{6}& 0& 0\\ 0& \sqrt{6}& 0& \sqrt{6}& 0\\ 0& 0& \sqrt{6}& 0& 2\\ 0& 0& 0& 2& 0\end{array}\right),f^z=\left(\begin{array}{ccccc}2& 0& 0& 0& 0\\ 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 2\end{array}\right).`$ (34)
We may use Eq. (17) to eliminate the last term in Eq. (16), obtaining

$`\widehat{V}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle }d𝒓[c_0:\widehat{n}^2(𝒓):+c_1:\widehat{𝑭}^2(𝒓):+c_2\widehat{A}_{00}^{}(𝒓)\widehat{A}_{00}(𝒓)],`$ (35)
where $`c_0(4g_2+3g_4)/7`$, $`c_1(g_4g_2)/7`$, and $`c_2(7g_010g_2+3g_4)/7`$.

### B Total Hamiltonian

In the following discussions we shall assume that the external magnetic field is weak enough to ignore the quadratic Zeeman effect. Then the total Hamiltonian $`\widehat{H}`$ of the system consists of the kinetic energy term $`\widehat{H}_{\mathrm{KE}}`$, the trapping potential energy term $`\widehat{H}_{\mathrm{PE}}`$, the linear Zeeman term $`\widehat{H}_{\mathrm{LZ}}`$, and the interaction term $`\widehat{V}`$:

$`\widehat{H}_{\mathrm{KE}}`$ $`=`$ $`{\displaystyle 𝑑𝒓\underset{m=2}{\overset{2}{}}\widehat{\mathrm{\Psi }}_m^{}\left(\frac{\mathrm{}^2^2}{2M}\right)\widehat{\mathrm{\Psi }}_m},`$ (36)
$`\widehat{H}_{\mathrm{PE}}`$ $`=`$ $`{\displaystyle 𝑑𝒓\underset{m=2}{\overset{2}{}}\widehat{\mathrm{\Psi }}_m^{}U_{\mathrm{trap}}(𝒓)\widehat{\mathrm{\Psi }}_m},`$ (37)
$`\widehat{H}_{\mathrm{LZ}}`$ $`=`$ $`p{\displaystyle 𝑑𝒓\underset{m,n=2}{\overset{2}{}}\widehat{\mathrm{\Psi }}_mf_{mn}^z\widehat{\mathrm{\Psi }}_n}=p{\displaystyle 𝑑𝒓\underset{m=2}{\overset{2}{}}m\widehat{\mathrm{\Psi }}_m^{}\widehat{\mathrm{\Psi }}_m},`$ (38)
where $`p`$ $`(>0)`$ is the product of the gyromagnetic ratio and the external magnetic field which is assumed to be applied in the z-direction.

When the system is spatially uniform, i.e. $`U_{\mathrm{trap}}(𝒓)=0`$, it is convenient to expand the field operator in terms of plane waves:

$`\widehat{\mathrm{\Psi }}_m(𝒓)={\displaystyle \frac{1}{\sqrt{V}}}{\displaystyle \underset{𝒌}{}}\widehat{a}_{𝒌m}e^{i𝒌𝒓}.`$ (39)
Then the single-particle part of the Hamiltonian becomes

$`\widehat{H}_0`$ $`=`$ $`\widehat{H}_{\mathrm{KE}}+\widehat{H}_{\mathrm{LZ}}={\displaystyle \underset{𝒌,m}{}}(ϵ_𝒌pm)\widehat{n}_{𝒌,m},`$ (40)
where $`ϵ_𝒌\mathrm{}^2k^2/2M`$, and the interaction Hamiltonian becomes

$`\widehat{V}={\displaystyle \frac{1}{2V}}{\displaystyle \underset{𝒌}{}}:(c_0\widehat{\rho }_𝒌^{}\widehat{\rho }_𝒌+c_1\widehat{𝑭}_𝒌^{}\widehat{𝑭}_𝒌+c_2\widehat{A}_𝒌^{}\widehat{A}_𝒌):,`$ (41)
where

$`\widehat{\rho }_𝒌{\displaystyle \widehat{n}(𝒓)e^{i𝒌𝒓}𝑑𝒓}={\displaystyle \underset{𝒒m}{}}\widehat{a}_{𝒒,m}^{}\widehat{a}_{𝒒+𝒌,m},`$ (42)
$`\widehat{𝑭}_𝒌{\displaystyle \widehat{𝑭}(𝒓)e^{i𝒌𝒓}𝑑𝒓}={\displaystyle \underset{𝒒mn}{}}𝐟_{mn}\widehat{a}_{𝒒,m}^{}\widehat{a}_{𝒒+𝒌,n},`$ (43)
$`\widehat{A}_𝒌{\displaystyle \widehat{A}_{00}(𝒓)e^{i𝒌𝒓}𝑑𝒓}={\displaystyle \frac{1}{\sqrt{5}}}{\displaystyle \underset{𝒒,m}{}}(1)^m\widehat{a}_{𝒒,m}\widehat{a}_{𝒌𝒒,m}.`$ (44)
## III Mean-field theory

Mean field theory (MFT) of spin-2 BEC was discussed in Refs. . We present here a brief summary of as much of this theory as is relevant to later discussions. We shall also present some new results. When the system is uniform, BEC occurs in the $`𝒌=\mathrm{𝟎}`$ state. We therefore assume the following trial wave function

$`|𝜻={\displaystyle \frac{1}{\sqrt{N!}}}\left({\displaystyle \underset{m=2}{\overset{2}{}}}\zeta _m\widehat{a}_{\mathrm{𝟎},m}^{}\right)^N|\mathrm{vac},`$ (45)
where the complex amplitudes $`\zeta _m`$ are assumed to satisfy the normalization condition

$`{\displaystyle \underset{m=2}{\overset{2}{}}}|\zeta _m|^2=1.`$ (46)
The variational parameters $`\zeta _m`$ are determined so as to minimize the expectation value of $`\widehat{H}`$ over the state (45):

$`\widehat{H}={\displaystyle \frac{c_0}{2V}}N(N1)+{\displaystyle \frac{c_1}{2V}}N(N1)\widehat{𝒇}^2+{\displaystyle \frac{2c_2}{5V}}N(N1)|\widehat{s}_{}|^2pN\widehat{f}_z,`$ (47)
where

$`\widehat{𝒇}{\displaystyle \underset{mn}{}}𝒇_{mn}\zeta _m^{}\zeta _n,`$ (48)
$`\widehat{f}_z{\displaystyle \underset{m}{}}m|\zeta _m|^2,`$ (49)
$`\widehat{s}_{}{\displaystyle \frac{1}{2}}{\displaystyle \underset{m}{}}(1)^m\zeta _m\zeta _m.`$ (50)
When the external magnetic field is applied in the $`z`$-direction, only the $`z`$-component of $`\widehat{𝒇}`$ is nonzero. We thus obtain

$`\widehat{H}{\displaystyle \frac{c_0}{2V}}N(N1)+{\displaystyle \frac{N(N1)}{2V}}ϵ_0,`$ (51)
where

$`ϵ_0=c_1\widehat{f}_z^2+{\displaystyle \frac{4}{5}}c_2|\widehat{s}_{}|^2\stackrel{~}{p}\widehat{f}_z`$ (52)
with $`\stackrel{~}{p}2Vp/(N1)`$.

The mean-field solution should be determined so as to minimize $`ϵ_0`$ subject to the normalization condition (46):

$`{\displaystyle \frac{}{\zeta _m^{}}}\left(ϵ_0\lambda {\displaystyle \underset{m}{}}|\zeta _m|^2\right)=(2c_1\widehat{f}_z\stackrel{~}{p})m\zeta _m\lambda \zeta _m+{\displaystyle \frac{4}{5}}c_2(1)^m\widehat{s}_{}\zeta _m^{}=0,`$ (53)
where $`\lambda `$ is a Lagrange multiplier. Multiplying both sides by $`(1)^m\zeta _m`$ and summing over $`m`$ yield

$`\left(\lambda {\displaystyle \frac{2}{5}}c_2\right)\widehat{s}_{}=0.`$ (54)
On the other hand, multiplying both sides of Eq. (53) by $`\zeta _m^{}`$ and summing over $`m`$ yields

$`(2c_1\widehat{f}_z\stackrel{~}{p})\widehat{f}_z\lambda +{\displaystyle \frac{8}{5}}c_2|\widehat{s}_{}|^2=0.`$ (55)
### A Ferromangetic BEC

When $`\widehat{s}_{}=0`$, Eq. (53) gives

$`\left[(2c_1\widehat{f}_z\stackrel{~}{p})m\lambda \right]\zeta _m=0.`$ (56)
For nonzero components $`\zeta _m0`$, Eq. (56) gives

$`(2c_1\widehat{f}_z\stackrel{~}{p})m\lambda =0.`$ (57)
The case of $`\stackrel{~}{p}=2c_1\widehat{f}_z`$ will be discussed in Sec. III C. When $`\stackrel{~}{p}2c_1\widehat{f}_z`$, only one component can be nonzero. The mean-field solutions, magnetizations, and mean-field energies are therefore given by

$`𝜻=e^{i\varphi }(1,0,0,0,0),\widehat{f}_z=2,ϵ^\mathrm{F}=4c_12\stackrel{~}{p},`$ (58)
$`𝜻=e^{i\varphi }(0,1,0,0,0),\widehat{f}_z=1,ϵ^\mathrm{F}=c_1\stackrel{~}{p},`$ (59)
$`𝜻=e^{i\varphi }(0,0,0,1,0),\widehat{f}_z=1,ϵ^\mathrm{F}=c_1+\stackrel{~}{p},`$ (60)
$`𝜻=e^{i\varphi }(0,0,0,0,1),\widehat{f}_z=2,ϵ^\mathrm{F}=4c_1+2\stackrel{~}{p},`$ (61)
where $`\varphi `$ is an arbitrary global phase.

The ground state is degenerate with respect to the global phase $`\varphi `$. This represents the gauge invariance, i.e. conservation of the total number of particle and leads to a massless Goldstone mode, as will be shown in Sec. V B. The conservation of the spin angular momentum does not lead to a new Goldstone mode because in ferromagnets all spins are aligned in the same direction and therefore the total spin angular momentum has the same piece of information as the total number of particles.

### B Antiferromagnetic BEC

The antiferromagnetic (or polar) phase of a BEC is defined as the one having nonzero spin-singlet pair amplitude, $`\widehat{s}_{}0`$. When $`\widehat{s}_{}0`$, Eq. (54) gives $`\lambda =2c_2/5`$. Substituting this and Eq. (55) into Eq. (52) yields

$`ϵ_0^{\mathrm{AF}}={\displaystyle \frac{c_2}{5}}{\displaystyle \frac{\stackrel{~}{p}}{2}}\widehat{f}_z.`$ (62)
With $`\lambda =2c_2/5`$, Eq. (53) leads to

$`(2c_1\widehat{f}_z\stackrel{~}{p})\left[\left(2c_1m^2{\displaystyle \frac{2}{5}}c_2\right)\widehat{f}_z\stackrel{~}{p}m^2\right]\zeta _m=0.`$ (63)
When $`\stackrel{~}{p}2c_1\widehat{f}_z`$, the solutions of Eq. (63) is that only ($`\zeta _2`$, $`\zeta _2`$) or ($`\zeta _1`$, $`\zeta _1`$) or $`\zeta _0`$ is nonzero. Determining these values using conditions $`\widehat{f}_z=_mm|\zeta _m|^2`$ and Eq. (46), we obtain the mean-field solutions and the corresponding magnetizations as

$`𝜻`$ $`={\displaystyle \frac{1}{\sqrt{2}}}e^{i\varphi }(\sqrt{1+{\displaystyle \frac{\widehat{f}_z}{2}}},0,0,0,e^{i\chi }\sqrt{1{\displaystyle \frac{\widehat{f}_z}{2}}}),`$ $`\widehat{f}_z={\displaystyle \frac{2\stackrel{~}{p}}{4c_1c_2/5}},ϵ^{\mathrm{AF}}={\displaystyle \frac{c_2}{5}}{\displaystyle \frac{\stackrel{~}{p}^2}{4c_1c_2/5}}`$ (64)
$`𝜻`$ $`={\displaystyle \frac{1}{\sqrt{2}}}e^{i\varphi }(0,\sqrt{1+\widehat{f}_z},0,e^{i\chi _1}\sqrt{1\widehat{f}_z},0),`$ $`\widehat{f}_z={\displaystyle \frac{\stackrel{~}{p}}{2(c_1c_2/5)}},ϵ^{\mathrm{AF}}={\displaystyle \frac{c_2}{5}}{\displaystyle \frac{\stackrel{~}{p}^2}{4(c_1c_2/5)}}`$ (65)
$`𝜻`$ $`=e^{i\varphi }(0,0,1,0,0),`$ $`\widehat{f}_z=0,ϵ^{\mathrm{AF}}={\displaystyle \frac{c_2}{5}}.`$ (66)
The mean-field solutions (64) and (65) are degenerate with respect to two continuous phase variables, that is, the global phase $`\varphi `$ and the relative phase $`\chi _m=\varphi _m\varphi _m`$ ($`m=1,2`$) between the two nonvanishing amplitudes $`\zeta _{\pm m}`$. Corresponding to these two continuous degeneracies, we expect to have two Goldstone modes, as will be shown in Sec. V C.

When $`\stackrel{~}{p}=2c_1\widehat{f}_z`$, Eq. (55) with $`\lambda =2c_2/5`$ gives $`|\widehat{s}_{}|=1/2`$, which, together with (53), leads to $`\zeta _m=e^{2i\varphi _0}(1)^m\zeta _m^{}`$. Hence we have $`\widehat{f}_z=0`$. This is possible only when the external magnetic field is zero. The corresponding order parameter is given by

$`𝜻=e^{i\varphi _0}({\displaystyle \frac{e^{i(\varphi _2\varphi _0)}}{\sqrt{2}}}\mathrm{sin}\theta \mathrm{sin}\psi ,{\displaystyle \frac{e^{i(\varphi _1\varphi _0)}}{\sqrt{2}}}\mathrm{sin}\theta \mathrm{cos}\psi ,\mathrm{cos}\theta ,{\displaystyle \frac{e^{i(\varphi _1\varphi _0)}}{\sqrt{2}}}\mathrm{sin}\theta \mathrm{cos}\psi ,{\displaystyle \frac{e^{i(\varphi _2\varphi _0)}}{\sqrt{2}}}\mathrm{sin}\theta \mathrm{sin}\psi ).`$ (67)
This solution is degenerate with respect to five continuous variables: one global gauge ($`\varphi _0`$), two relative gauges ($`\varphi _2\varphi _0`$ and $`\varphi _1\varphi _0`$), and two variables $`\theta `$ and $`\psi `$ that specify the amplitudes of the order parameter. As a consequence of these degeneracies, we expect to have five (three density-like and two spin-like) Goldstone modes, as will be shown in Sec. V C.

### C Cyclic BEC

The remaining possibility is the case in which $`\widehat{s}_{}=0`$ and $`\widehat{f}_z=\stackrel{~}{p}/2c_1`$. This phase will be referred to as cyclic phase. The energy of this phase is given from Eq. (52) by

$`ϵ_0^\mathrm{C}=c_1\widehat{f}_z^2.`$ (68)
Let us now parameterize the order parameter of the cyclic phase as it will be needed to find the Bogoliubov spectrum. There are four equations (six real equations) that restrict the order parameter of this phase, that is,

$`{\displaystyle \underset{m}{}}|\zeta _m|^2=1,`$ (69)
$`\widehat{s}_{}={\displaystyle \frac{1}{2}}{\displaystyle \underset{m}{}}(1)^m\zeta _m\zeta _m=0,`$ (70)
$`\widehat{f}_z={\displaystyle \underset{m}{}}m|\zeta _m|^2={\displaystyle \frac{\stackrel{~}{p}}{2c_1}},`$ (71)
$`\widehat{f}_+=2(\zeta _2^{}\zeta _1+\zeta _1^{}\zeta _2)+\sqrt{6}(\zeta _1^{}\zeta _0+\zeta _0^{}\zeta _1)=0.`$ (72)
We use the representation (A16) of the order parameter derived in Appendix A to analyze the cyclic phase. This representation automatically satisfies the normalization condition (69). To meet the condition (70), we note that

$`\widehat{s}_{}={\displaystyle \frac{1}{2}}\mathrm{Tr}𝐌^2={\displaystyle \frac{1}{2}}\zeta _0^2\zeta _1\zeta _1+\zeta _2\zeta _2={\displaystyle \frac{i}{2}}\mathrm{sin}\chi \mathrm{cos}(\varphi \theta ).`$ (73)
The condition (70) therefore requires either

$`\chi =0\mathrm{or}\varphi \theta =\pi /2.`$ (74)
On the other hand, Eq. (72) becomes

$`\widehat{f}_+=2\mathrm{sin}\delta \mathrm{cos}\chi \left[\mathrm{cos}\psi \mathrm{sin}(\theta +\pi /6)i\mathrm{sin}\psi \mathrm{sin}(\theta \pi /6)\right]=0,`$ (75)
whence we obtain

$`(\mathrm{i})\chi =\pi /2,\mathrm{or}(\mathrm{ii})\delta =0,\mathrm{or}(\mathrm{iii})\psi =\pi /2\mathrm{and}\theta =\pi /6,\mathrm{or}(\mathrm{iv})\psi =0\mathrm{and}\theta =\pi /6.`$ (76)
From conditions (74) and (76), we find the following three solutions and the corresponding magnetization:

$`\zeta _{\pm 2}`$ $`={\displaystyle \frac{1}{2}}(\mathrm{cos}\theta \pm \mathrm{cos}\chi i\mathrm{sin}\theta \mathrm{sin}\chi )e^{i\varphi },\zeta _{\pm 1}=0,\zeta _0={\displaystyle \frac{1}{\sqrt{2}}}(\mathrm{sin}\theta +i\mathrm{cos}\theta \mathrm{sin}\chi )e^{i\varphi },`$ $`\widehat{f}_z=2\mathrm{cos}\chi \mathrm{cos}\theta ,`$ (77)
$`\zeta _{\pm 2}`$ $`={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\sqrt{3}}{2}}\pm \mathrm{cos}\delta \mathrm{cos}\chi {\displaystyle \frac{i}{2}}\mathrm{sin}\chi \right)e^{i\varphi },\zeta _{\pm 1}=\pm {\displaystyle \frac{i}{2}}e^{i\varphi }\mathrm{sin}\delta \mathrm{cos}\chi ,\zeta _0={\displaystyle \frac{1}{\sqrt{2}}}\left({\displaystyle \frac{1}{2}}+i{\displaystyle \frac{\sqrt{3}}{2}}\mathrm{sin}\chi \right)e^{i\varphi },`$ $`\widehat{f}_z=\sqrt{3}\mathrm{cos}\delta \mathrm{cos}\chi ,`$ (78)
$`\zeta _{\pm 2}`$ $`={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\sqrt{3}}{2}}\pm \mathrm{cos}\delta \mathrm{cos}\chi +{\displaystyle \frac{i}{2}}\mathrm{sin}\chi \right)e^{i\varphi },\zeta _{\pm 1}={\displaystyle \frac{1}{2}}e^{i\varphi }\mathrm{sin}\delta \mathrm{cos}\chi ,\zeta _0={\displaystyle \frac{1}{\sqrt{2}}}\left({\displaystyle \frac{1}{2}}+i{\displaystyle \frac{\sqrt{3}}{2}}\mathrm{sin}\chi \right)e^{i\varphi },`$ $`\widehat{f}_z=\sqrt{3}\mathrm{cos}\delta \mathrm{cos}\chi ,`$ (79)
where a global phase $`\varphi `$, which is chosen to be a particular value in Eq.(A9), is recovered. While these solutions include three parameters, the condition (71) leaves only two parameters free. It should be noted that these three solutions give the same ground-state energy, and hence are equally possible unless magnetization exceeds $`\sqrt{3}`$. When it exceeds $`\sqrt{3}`$, only solution (77) is possible.

A remark is here in order. As the representation (A16) is obtained by assuming the full isotropy of space, it does not cover the whole order parameter space in the presence of magnetic field. We should therefore keep in mind that the solutions (77)-(79) do not necessarily exhaust the whole order parameter space of the cyclic phase.

In the following discussions we shall focus on the solution (77). Making the absolute square of each amplitude yields

$`|\zeta _{\pm 2}|^2={\displaystyle \frac{1}{4}}(1\pm \widehat{f}_z/2)^2,|\zeta _0|^2={\displaystyle \frac{1}{2}}(1\widehat{f}_z^2/4).`$ (80)
Hence we obtain

$`\zeta _{\pm 2}={\displaystyle \frac{1}{2}}\left(1\pm {\displaystyle \frac{\widehat{f}_z}{2}}\right)e^{i\varphi _{\pm 2}},\zeta _{\pm 1}=0,\zeta _0=\sqrt{{\displaystyle \frac{1\widehat{f}_z^2/4}{2}}}e^{i\varphi _0},`$ (81)
where

$`\mathrm{tan}\varphi _{\pm 2}={\displaystyle \frac{\mathrm{sin}\chi \mathrm{sin}\theta }{\mathrm{cos}\theta \pm \mathrm{cos}\chi }},\mathrm{tan}\varphi _0=\mathrm{sin}\chi \mathrm{cot}\theta .`$ (82)
It follows from this or by direct calculation that

$`\varphi _2+\varphi _22\varphi _0=\pm \pi .`$ (83)
Because of this restriction, the ground state of the cyclic phase is degenerate with respect to at least two continuous phase variables. We therefore expect to have at least two Goldstone modes, as will be shown in Sec.V D.

### D The phase boundaries

In the absence of external magnetic field, the ground-state energies for the three phases are given from Eqs. (58), (62), and (68) by

$`ϵ_0^\mathrm{F}=4c_1,ϵ_0^{\mathrm{AF}}={\displaystyle \frac{c_2}{5}},ϵ_0^\mathrm{C}=0.`$ (84)
It follows that each phase is specified by .

$`\mathrm{ferromagnetic}`$ $`c_1<0\mathrm{and}c_1c_2/20<0,`$ (85)
$`\mathrm{antiferromagnetic}`$ $`c_2<0\mathrm{and}c_1c_2/20>0,`$ (86)
$`\mathrm{cyclic}`$ $`c_1>0\mathrm{and}c_2>0.`$ (87)
In the presence of external magnetic field we define each phase as follows:

$`\mathrm{ferromagnetic}`$ $`\widehat{f}_z=2,`$
$`\mathrm{antiferromagnetic}`$ $`\widehat{s}_{}0,`$
$`\mathrm{cyclic}`$ $`\widehat{s}_{}=0\mathrm{and}\widehat{f}_z<2.`$
By directly comparing the energies in Eqs. (62), (58), and (68), we find that each phase is specified by the following conditions:

$`\mathrm{ferromagnetic}`$ $`c_1\stackrel{~}{p}/4\mathrm{and}c_1c_2/20\stackrel{~}{p}/4,`$ (88)
$`\mathrm{antiferromagnetic}`$ $`c_2<0\mathrm{and}c_1c_2/20>\stackrel{~}{p}/4,`$ (89)
$`\mathrm{cyclic}`$ $`c_1>\stackrel{~}{p}/4\mathrm{and}c_2>0.`$ (90)
## IV Many-body spin correlations and magnetic response

In this section we study the case in which the system is so tightly confined that the coordinate part of the ground-state wave function $`\varphi _0(𝒓)`$ is independent of the spin state and solely determined by $`\widehat{H}_{\mathrm{KE}}`$, $`\widehat{H}_{\mathrm{PE}}`$, and the spin-independent part of $`\widehat{V}`$; that is, $`\varphi _0(𝒓)`$ is the solution $`\varphi `$ to the equation

$`\left[{\displaystyle \frac{\mathrm{}^2\mathbf{}^2}{2M}}+U_{\mathrm{trap}}+c_0(N1)|\varphi |^2\right]\varphi =ϵ\varphi `$ (91)
with the lowest eigenvalue $`ϵ=ϵ_0`$. This assumption is justified if the second lowest eigenvalue $`ϵ_1`$ satisfies

$`ϵ_1ϵ_0|p|,|c_1|N/V^{\mathrm{eff}},|c_2|N/V^{\mathrm{eff}},`$ (92)
where $`V^{\mathrm{eff}}(𝑑𝒓|\varphi _0|^4)^1`$ is an effective volume which coincides with $`V`$ in Eq. (39) for the spatially uniform case (i.e., $`U_{\mathrm{trap}}=0`$). When the condition (92) is met, the field operator $`\widehat{\mathrm{\Psi }}_m`$ may be approximated as $`\widehat{\mathrm{\Psi }}_m\widehat{a}_m\varphi _0`$, where $`\widehat{a}_m`$ is the annihilation operator of the bosons that are specified by the spin component $`m`$ and by the coordinate wave function $`\varphi _0`$. The spin-dependent part of the Hamiltonian can then be written as

$$\widehat{H}=\frac{c_1}{2V^{\mathrm{eff}}}:\widehat{𝓕}^2:+\frac{2c_2}{5V^{\mathrm{eff}}}\widehat{𝒮}_+\widehat{𝒮}_{}p\widehat{}_z,$$
(93)
where

$`\widehat{𝓕}{\displaystyle \underset{mn}{}}𝒇_{mn}\widehat{a}_m^{}\widehat{a}_n,\widehat{𝒮}_{}=(\widehat{𝒮}_+)^{}{\displaystyle \frac{1}{2}}{\displaystyle \underset{m}{}}(1)^m\widehat{a}_m\widehat{a}_m,\widehat{}_z{\displaystyle \underset{m}{}}m\widehat{a}_m^{}\widehat{a}_m.`$ (94)
### A Spectrum and degeneracy

We first make some remarks on the properties of the operators $`\widehat{𝒮}_{}=\widehat{𝒮}_+^{}=(\widehat{a}_0)^2/2\widehat{a}_1\widehat{a}_1+\widehat{a}_2\widehat{a}_2`$. The operator $`\widehat{𝒮}_+`$, when applied to the vacuum, creates a pair of bosons in the spin-singlet state. This pair, however, should not be regarded as a single composite boson because $`\widehat{𝒮}_+`$ does not satisfy the commutation relations of bosons. The operator $`\widehat{𝒮}_+`$ instead satisfies the $`SU(1,1)`$ commutation relations together with $`\widehat{𝒮}_z(2\widehat{N}+5)/4`$, namely,

$`[\widehat{𝒮}_z,\widehat{𝒮}_\pm ]=\pm \widehat{𝒮}_\pm ,[\widehat{𝒮}_+,\widehat{𝒮}_{}]=2\widehat{𝒮}_z,`$ (95)
where the minus sign in the last equation is the only distinction from the usual spin commutation relations. As a consequence, the Casimir operator $`\widehat{𝒮}^2`$ that commutes with $`\widehat{𝒮}_\pm `$ and $`\widehat{𝒮}_z`$ is given by

$`\widehat{𝒮}^2\widehat{𝒮}_+\widehat{𝒮}_{}+\widehat{𝒮}_z^2\widehat{𝒮}_z.`$ (96)
Consider an eigenspace $`_\nu `$ of $`\widehat{𝒮}^2`$ with an eigenvalue $`\nu `$. The requirement that $`\widehat{𝒮}_+\widehat{𝒮}_{}=\widehat{𝒮}_z^2\widehat{𝒮}_z\widehat{𝒮}^2`$ must be positive semidefinite means that, in $`_\nu `$, the eigenvalues of the operator $`\widehat{𝒮}_z`$ has a minimum value $`(2N_0+5)/4`$, where $`N_0`$ is a nonnegative integer. For a state $`|\varphi `$ that belongs to the minimum eigenvalue, the norm of $`\widehat{𝒮}_{}|\varphi `$ must vanish; hence $`\nu =S(S1)`$ with $`S=(2N_0+5)/4`$. We thus obtain the allowed combinations of eigenvalues $`\{S(S1),S_z\}`$ for $`\widehat{𝒮}^2`$ and $`\widehat{𝒮}_z`$ such that

$$S=(2N_0+5)/4(N_0=0,1,2,\mathrm{})$$
(97)
and

$$S_z=S+N_\mathrm{S}(N_\mathrm{S}=0,1,2,\mathrm{}).$$
(98)
Here we have introduced quantum numbers $`N_\mathrm{S}`$ and $`N_0`$, where the operator $`\widehat{𝒮}_+`$ raises $`N_\mathrm{S}`$ by one and the relation

$$N=2N_\mathrm{S}+N_0$$
(99)
holds. We may thus interpret $`N_\mathrm{S}`$ as the number of spin-singlet ‘pairs’, and $`N_0`$ as that of the remaining bosons.

Exact energy eigenvalues of Hamiltonian (93) can be obtained as follows. The operators $`\widehat{𝒮}_\pm `$ are invariant under any rotation of the system, namely, they commute with the total spin operator $`\widehat{𝓕}`$. The energy eigenstates can thus be classified according to quantum numbers $`N_0`$ and $`N_\mathrm{S}`$, total spin $`F`$, and magnetic quantum number $`F_z`$. We thus denote the eigenstates as $`|N_0,N_\mathrm{S},F,F_z;\lambda `$, where $`\lambda =1,2,\mathrm{},g_{N_0,F}`$ labels orthonormal degenerate states, that is,

$`N_0,N_\mathrm{S},F,F_z;\lambda ^{}|N_0,N_\mathrm{S},F,F_z;\lambda =\delta _{\lambda ^{},\lambda }.`$ (100)
The number of degenerate states $`g_{N_0,F}`$ for a given set of $`\{N_0,N_\mathrm{S},F,F_z\}`$ will be referred to as the size of the eigenspaces for $`\{N_0,N_\mathrm{S},F,F_z\}`$. It will be shown to be independent of $`N_\mathrm{S}`$ and $`F_z`$ below. The energy eigenvalue for the state $`|N_0,N_\mathrm{S},F,F_z;\lambda `$ is given by

$$E=\frac{c_1}{2V^{\mathrm{eff}}}[F(F+1)6N]+\frac{c_2}{5V^{\mathrm{eff}}}N_\mathrm{S}(N+N_0+3)pF_z,$$
(101)
where the relation $`2N_\mathrm{S}+N_0=N`$ is used.

The degeneracy $`g_{N_0,F}`$ can be calculated as follows. First we show that $`g_{N_0,F}`$ is independent of $`N_\mathrm{S}`$ and $`F_z`$. This is seen by the following relations,

$$N_0,N_\mathrm{S},F,F_z;\lambda ^{}|\widehat{}_\pm \widehat{}_{}|N_0,N_\mathrm{S},F,F_z;\lambda =(F\pm F_z)(FF_z+1)\delta _{\lambda ^{},\lambda },$$
(102)
where $`\widehat{}_\pm \widehat{}_x\pm i\widehat{}_y`$,

$$N_0,N_\mathrm{S},F,F_z;\lambda ^{}|\widehat{𝒮}_+\widehat{𝒮}_{}|N_0,N_\mathrm{S},F,F_z;\lambda =N_\mathrm{S}(N_0+N_\mathrm{S}+3/2)\delta _{\lambda ^{},\lambda },$$
(103)
and

$$N_0,N_\mathrm{S},F,F_z;\lambda ^{}|\widehat{𝒮}_{}\widehat{𝒮}_+|N_0,N_\mathrm{S},F,F_z;\lambda =(N_\mathrm{S}+1)(N_0+N_\mathrm{S}+5/2)\delta _{\lambda ^{},\lambda }.$$
(104)
These relations implies that the sizes of the eigenspaces for $`\{N_0,N_\mathrm{S}\pm 1,F,F_z\pm 1\}`$ are not smaller than the size of the eigenspace for $`\{N_0,N_\mathrm{S},F,F_z\}`$. The degeneracy thus depends only on $`N_0`$ and $`F`$. Next, we introduce a generating function of $`g_{N_0,F}`$ defined by

$$G(x,y)\underset{N_0=0}{\overset{\mathrm{}}{}}\underset{F=0}{\overset{\mathrm{}}{}}g_{N_0,F}x^{N_0}y^F.$$
(105)
Let $`h_{N,F_z}`$ be the total number of states with a fixed number of bosons $`N`$ and a fixed magnetic quantum number $`F_z`$. This is given by the total number of combinations of nonnegative integers $`\{n_2,n_1,n_0,n_1,n_2\}`$ that satisfy $`n_2+n_1+n_0+n_1+n_2=N`$ and $`2n_2n_1+n_1+2n_2=F_z`$. It follows that

$$\underset{N=0}{\overset{\mathrm{}}{}}\underset{F_z=\mathrm{}}{\overset{\mathrm{}}{}}h_{N,F_z}z^Ny^{2N+F_z+1}=\underset{\{n_j\}}{}\underset{j=2}{\overset{2}{}}z^{n_j}y^{(j+2)n_j+1}=y\underset{j=2}{\overset{2}{}}(1zy^{j+2})^1,$$
(106)
where we assume $`|y|<1`$ and $`|z|<1`$ to ensure the convergence of the series. Let $`\stackrel{~}{h}_{N,F,F_z}`$ be the total number of states for given $`N`$, $`F`$, and $`F_z`$. Because $`\stackrel{~}{h}_{N,F,F_z}`$ is independent of the value of $`F_z`$, we shall denote it simply as $`\stackrel{~}{h}_{N,F}`$.

The quantity $`h_{N,F_z}`$ is written in terms of the sum of $`\stackrel{~}{h}_{N,F}`$ as

$$h_{N,F_z}=\underset{F(|F_z|)}{}\stackrel{~}{h}_{N,F},$$
(107)
and hence $`\stackrel{~}{h}_{N,F}=h_{N,F}h_{N,F+1}`$. Let us extend the definition of $`\stackrel{~}{h}_{N,F}`$ to the negative values of $`F`$ through this relation. It follows then from Eq. (106) that

$$\underset{N=0}{\overset{\mathrm{}}{}}\underset{F=\mathrm{}}{\overset{\mathrm{}}{}}\stackrel{~}{h}_{N,F}z^Ny^{2N+F+1}=(y1)\underset{j=2}{\overset{2}{}}(1zy^{j+2})^1.$$
(108)
The right hand side of this equation can be written as the sum of two fractions $`G_1(z,y)+G_2(z,y)`$, where

$$G_1(z,y)=\frac{y(1zy^3+z^2y^6)}{(1zy^3)(1zy^4)(1z^2y^4)(1z^3y^6)}$$
(109)
and

$$G_2(z,y)=\frac{1zy+z^2y^2}{(1z)(1zy)(1z^2y^4)(1z^3y^6)}.$$
(110)
Making Maclaurin expansions of $`G_1`$ and $`G_2`$ around $`z=y=0`$ and regrouping them in terms of the form $`z^Ny^{2N+F+1}`$, we find that $`G_1`$ consists only of the terms with $`F0`$, and $`G_2`$ of those with $`F<0`$. We thus obtain

$$\underset{N=0}{\overset{\mathrm{}}{}}\underset{F=0}{\overset{\mathrm{}}{}}\stackrel{~}{h}_{N,F}z^Ny^{2N+F+1}=G_1(z,y).$$
(111)
The quantity $`\stackrel{~}{h}_{N,F}`$ is written by the sum of the degeneracy $`g_{N_0,F}`$ as

$$\stackrel{~}{h}_{N,F}=\underset{N_\mathrm{S}(N/2)}{}g_{N2N_\mathrm{S},F},$$
(112)
and hence we can write $`g_{N_0,F}=\stackrel{~}{h}_{N_0,F}\stackrel{~}{h}_{N_02,F}`$, where we assume that $`\stackrel{~}{h}_{1,F}=\stackrel{~}{h}_{2,F}=0`$. It follows then from Eq. (111) that

$$\underset{N_0=0}{\overset{\mathrm{}}{}}\underset{F=\mathrm{}}{\overset{\mathrm{}}{}}g_{N_0,F}z^{N_0}y^{2N_0+F+1}=(1z^2y^4)G_1(z,y),$$
(113)
and we finally obtain an explicit form of the generating function $`G(x,y)`$ defined by Eq. (105) as

$$\underset{N_0=0}{\overset{\mathrm{}}{}}\underset{F=0}{\overset{\mathrm{}}{}}g_{N_0,F}x^{N_0}y^F=y^1(1x^2)G_1(xy^2,y)=\frac{1xy+x^2y^2}{(1xy)(1xy^2)(1x^3)}.$$
(114)
The total spin $`F`$ can, in general, take integer values in the range $`0F2N_0`$. However, from Eq. (114) we find that there are some forbidden values. That is, $`F=1,2,5,2N_01`$ are not allowed when $`N_0=3k(k𝒁)`$, and $`F=0,1,3,2N_01`$ are forbidden when $`N_0=3k\pm 1`$. An easier way to find the forbidden values is discussed at the end of Sec. IV B.

### B Energy eigenstates

The energy eigenstates $`|N_0,N_\mathrm{S},F,F_z;\lambda `$ can be constructed using one-, two- and three-boson creation operators. Let us define the operator $`\widehat{A}_f^{(n)}^{}`$ such that it creates $`n`$ bosons in the state with total spin $`F=f`$ and magnetic quantum number $`F_z=f`$ when applied to the vacuum. Such states are unique when $`n3`$. Among possible operators $`\widehat{A}_f^{(n)}`$, we choose the following five operators for constructing the eigenstates:

$`\widehat{A}_2^{(1)}`$ $`=`$ $`\widehat{a}_2^{}`$ (115)
$`\widehat{A}_0^{(2)}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{10}}}[(\widehat{a}_0^{})^22\widehat{a}_1^{}\widehat{a}_1^{}+2\widehat{a}_2^{}\widehat{a}_2^{}]=\sqrt{{\displaystyle \frac{2}{5}}}\widehat{𝒮}_+`$ (116)
$`\widehat{A}_2^{(2)}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{14}}}[2\sqrt{2}\widehat{a}_2^{}\widehat{a}_0^{}\sqrt{3}(\widehat{a}_1^{})^2]`$ (117)
$`\widehat{A}_0^{(3)}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{210}}}[\sqrt{2}(\widehat{a}_0^{})^33\sqrt{2}\widehat{a}_1\widehat{a}_0^{}\widehat{a}_1^{}+3\sqrt{3}(\widehat{a}_1^{})^2\widehat{a}_2^{}+3\sqrt{3}\widehat{a}_2^{}(\widehat{a}_1^{})^26\sqrt{2}\widehat{a}_2^{}\widehat{a}_0^{}\widehat{a}_2^{}]`$ (118)
$`\widehat{A}_3^{(3)}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{20}}}[(\widehat{a}_1^{})^3\sqrt{6}\widehat{a}_2^{}\widehat{a}_1^{}\widehat{a}_0^{}+2(\widehat{a}_2^{})^2\widehat{a}_1^{}].`$ (119)
Note that $`\widehat{A}_1^{(2)}`$ and $`\widehat{A}_1^{(3)}`$ do not exist because of the Bose symmetry. Note also that the operators $`\widehat{A}_f^{(n)}`$ commute with $`\widehat{F}_+`$.

Consider a set $``$ of unnormalized states,

$$|n_{12},n_{20},n_{22},n_{30},n_{33}(\widehat{a}_2^{})^{n_{12}}(\widehat{A}_0^{(2)})^{n_{20}}(\widehat{A}_2^{(2)})^{n_{22}}(\widehat{A}_0^{(3)})^{n_{30}}(\widehat{A}_3^{(3)})^{n_{33}}|\mathrm{vac}$$
(120)
with $`n_{12},n_{20},n_{22},n_{30}=0,1,2,\mathrm{},\mathrm{}`$ and $`n_{33}=0,1`$. The state $`|n_{12},n_{20},n_{22},n_{30},n_{33}`$ has the total number of bosons $`N=n_{12}+2(n_{20}+n_{22})+3(n_{30}+n_{33})`$ and the total spin $`F=F_z=2(n_{12}+n_{22})+3n_{33}`$. When $`N`$ and $`F`$ are given, $`n_{33}`$ is uniquely determined through the parity of $`F`$, namely,

$`n_{33}=F\mathrm{mod}\mathrm{\hspace{0.33em}2}.`$ (121)
If we introduce the following two parameters

$`\mu `$ $``$ $`2n_{20}+3n_{30},`$ (122)
$`\nu `$ $``$ $`n_{20}+n_{30},`$ (123)
$`n_{20}`$ and $`n_{30}`$ are uniquely specified by them, and the remaining $`n_{12}`$ and $`n_{22}`$ are also determined as

$`n_{12}`$ $`=`$ $`FN+\mu ,`$ (124)
$`n_{22}`$ $`=`$ $`N{\displaystyle \frac{F}{2}}\mu {\displaystyle \frac{3}{2}}n_{33}.`$ (125)
The parameter set $`\{n_{12},n_{20},n_{22},n_{30},n_{33}\}`$ is thus uniquely specified by the set $`\{N,F,\mu ,\nu \}`$. Let us consider the two commuting observables

$`\widehat{n}_\mu `$ $``$ $`2\widehat{a}_2^{}\widehat{a}_2+{\displaystyle \frac{3}{2}}(\widehat{a}_1^{}\widehat{a}_1\widehat{n}_{33})`$ (126)
$`\widehat{n}_\nu `$ $``$ $`\widehat{a}_2^{}\widehat{a}_2+{\displaystyle \frac{1}{2}}(\widehat{a}_1^{}\widehat{a}_1\widehat{n}_{33})`$ (127)
with $`\widehat{n}_{33}\widehat{F}\mathrm{mod}\mathrm{\hspace{0.33em}2}`$. Let $`\widehat{P}(n_\mu ,n_\nu )`$ be the projection operator onto the simultaneous eigenspace of $`\widehat{n}_\mu `$ and $`\widehat{n}_\nu `$ corresponding to eigenvalues $`n_\mu `$ and $`n_\nu `$, respectively. Noting that $`\widehat{A}_0^{(2)}`$ includes the term $`\widehat{a}_2\widehat{a}_2`$ and $`\widehat{A}_0^{(3)}`$ includes $`\widehat{a}_2(\widehat{a}_1)^2`$, it can be seen from Eqs. (115)-(120) that

$$\widehat{P}(n_\mu ,n_\nu )|n_{12},n_{20},n_{22},n_{30},n_{33}\mathrm{𝟎}\mathrm{if}n_\mu =\mu \mathrm{and}n_\nu =\nu ,$$
(128)
and also

$$\widehat{P}(n_\mu ,n_\nu )|n_{12},n_{20},n_{22},n_{30},n_{33}=\mathrm{𝟎}\mathrm{if}n_\mu >\mu \mathrm{or}n_\nu >\nu .$$
(129)
If we order the pair $`(\mu ,\nu )`$ in the lexicographic way, namely, $`(\mu ,\nu )<(\mu ^{},\nu ^{})`$ if $`\mu <\mu ^{}`$ or $`(\mu =\mu ^{}`$ and $`\nu <\nu ^{})`$, the above equations imply that a state in $``$ specified by $`\{N,F,\mu ,\nu \}`$ is linearly independent of the set of states specified by $`\{N,F,\mu ^{},\nu ^{}\}`$ with $`(\mu ^{},\nu ^{})<(\mu ,\nu )`$. The set $``$ is thus a linearly independent set of states. Let $`\stackrel{~}{h}_{N,F}^{}`$ be the total number of states in $``$ with $`N`$ bosons, total spin $`F`$, and magnetic quantum number $`F_z=F`$. A generating function of $`\stackrel{~}{h}_{N,F}^{}`$ is calculated as

$`{\displaystyle \underset{N=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{F=0}{\overset{\mathrm{}}{}}}\stackrel{~}{h}_{N,F}^{}x^Ny^F`$ $`=`$ $`{\displaystyle \underset{n_{12},n_{20},n_{22},n_{30}=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{n_{33}=0}{\overset{1}{}}}x^{n_{12}+2(n_{20}+n_{22})+3(n_{30}+n_{33})}y^{2(n_{12}+n_{22})+3n_{33}}`$ (130)
$`=`$ $`{\displaystyle \frac{1+x^3y^3}{(1xy^2)(1x^2)(1x^2y^2)(1x^3)}}=y^1G_1(xy^2,y),`$ (131)
where $`G_1`$ is defined by Eq. (109). Compared to Eq. (111), we have $`\stackrel{~}{h}_{N,F}^{}=\stackrel{~}{h}_{N,F}`$. This implies that $``$ is complete, namely, the set $``$ forms a nonorthogonal basis of the subspace $`_{(F_z=F)}`$ in which magnetic quantum number $`F_z`$ is equal to total spin $`F`$.

The energy eigenstates can be obtained by partially applying the method of Schmidt’s orthogonalization to the nonorthogonal basis $``$. Here by “partially” we mean that eigenstates corresponding to different energies are orthogonal, but that those corresponding to the same energy are not always so. Let us consider a series of subspaces $`_{(F_z=F)}=^{(0)}^{(1)}\mathrm{}`$, where $`^{(j)}`$ is spanned by the states with quantum number $`N_\mathrm{S}`$ \[see Eq. (98)\] satisfying $`N_\mathrm{S}j`$. Let $`\widehat{P}^{(j)}`$ be the projection operator onto $`^{(j)}`$. From the relation $`\widehat{𝒮}_+\widehat{P}^{(j)}=\widehat{P}^{(j+1)}\widehat{𝒮}_+`$, we have $`\widehat{P}^{(j)}(\widehat{A}_0^{(2)})^j=(\widehat{A}_0^{(2)})^j\widehat{P}^{(0)}`$. Since $`\widehat{P}^{(0)}|\varphi =|\varphi `$ for any state $`|\varphi `$, we have

$$(\widehat{P}^{(0)}\widehat{P}^{(n_{20})})|n_{12},n_{20},n_{22},n_{30},n_{33}=\mathrm{𝟎},$$
(132)
implying that $`^{(j)}`$ is spanned by all the states $`\{|n_{12},n_{20},n_{22},n_{30},n_{33}\}`$ satisfying $`n_{20}j`$. We can then construct a new basis $`^{}`$ made up of the states of the form

$$(\widehat{P}^{(0)}\widehat{P}^{(n_{20}+1)})|n_{12},n_{20},n_{22},n_{30},n_{33}=(\widehat{A}_0^{(2)})^{n_{20}}\widehat{P}_{(N_\mathrm{S}=0)}(\widehat{a}_2^{})^{n_{12}}(\widehat{A}_2^{(2)})^{n_{22}}(\widehat{A}_0^{(3)})^{n_{30}}(\widehat{A}_3^{(3)})^{n_{33}}|\mathrm{vac},$$
(133)
where $`\widehat{P}_{(N_\mathrm{S}=0)}\widehat{P}^{(0)}\widehat{P}^{(1)}`$ is the projection onto the subspace with $`N_\mathrm{S}=0`$ (the kernel of $`\widehat{𝒮}_{}`$). It is easy to see that the states belonging to $`^{}`$ are simultaneous eigenstates of $`\{\widehat{N},\widehat{𝒮}_+\widehat{𝒮}_{},\widehat{𝓕}^2,\widehat{}_z\}`$, and hence energy eigenstates. The energy eigenstates with $`F_z<F`$ can be constructed by applying $`(\widehat{}_{})^{FF_z}`$ to the states of $`^{}`$.

To summarize, the energy eigenstates can be represented as

$$(\widehat{}_{})^{\mathrm{\Delta }F}(\widehat{A}_0^{(2)})^{n_{20}}\widehat{P}_{(N_\mathrm{S}=0)}(\widehat{a}_2^{})^{n_{12}}(\widehat{A}_2^{(2)})^{n_{22}}(\widehat{A}_0^{(3)})^{n_{30}}(\widehat{A}_3^{(3)})^{n_{33}}|\mathrm{vac},$$
(134)
with $`n_{12},n_{20},n_{22},n_{30}=0,1,2,\mathrm{},\mathrm{}`$, $`n_{33}=0,1`$, and $`\mathrm{\Delta }F=0,1,\mathrm{},2F`$. These parameters are related to $`\{N_0,N_\mathrm{S},F,F_z\}`$ as

$`N_0`$ $`=`$ $`n_{12}+2n_{22}+3n_{30}+3n_{33}`$ (135)
$`N_\mathrm{S}`$ $`=`$ $`n_{20}`$ (136)
$`F`$ $`=`$ $`2n_{12}+2n_{22}+3n_{33}`$ (137)
$`F_z`$ $`=`$ $`F\mathrm{\Delta }F,`$ (138)
and the corresponding eigenenergy is given by Eq. (101). Note that the states defined in (134) are unnormalized, and the states having the same energy (i.e., those belonging to the same set of parameter values $`\{N_0,N_\mathrm{S},F,F_z\}`$) are nonorthogonal.

The representation (134) of the energy eigenstates utilizes the operator $`\widehat{A}_f^{(n)}`$ that creates correlated $`n`$ bosons having total spin $`f`$. It might be tempting to envisage a physical picture that the system is, like in 4He, made up of $`n_{nf}`$ composite bosons whose creation operator is given by $`\widehat{A}_f^{(n)}`$. However, this picture is oversimplified. First of all, the operator $`\widehat{A}_f^{(n)}`$ does not obey the boson commutation relation. In addition, the projection operator $`\widehat{P}_{(N_\mathrm{S}=0)}`$ in (134) imposes many-body spin correlations such that the spin correlation between any two bosons must have vanishing spin-singlet component. Note that two bosons with independently fluctuating spins have a nonzero overlap with the spin-singlet state in general. The many-body spin correlations of the energy eigenstates are thus far more complicated than what an intuitive picture of composite bosons suggests. On the other hand, as long as quantities such as the number of bosons, magnetization, and energy are concerned, the above simplified picture is quite helpful. By way of illustration, we provide an alternative explanation for the existence of forbidden values for the total spin $`F`$, which were found earlier using the generating function (114). For example, to construct a state with $`F=0`$ or $`F=3`$, composite particles with total spin 2 must be avoided, namely, $`n_{12}=n_{22}=0`$. Then we have $`N_0=3(n_{30}+n_{33})`$, implying that $`F=0`$ or $`F=3`$ is only possible when $`N_0=3k(k𝒁)`$. For a state with $`F=2`$ or $`F=5`$, we have $`n_{12}+n_{22}=1`$ and $`N_0=1+n_{22}+3(n_{30}+n_{33})`$ implying that $`N_03k(k𝒁)`$. The above simplified picture is also helpful when we consider the magnetic response as discussed below.

### C Magnetic response

We consider here how the ground state and the magnetization $`F_z`$ respond to the applied magnetic field $`p`$. From Eq. (101), we see that the minimum energy states always satisfy $`F_z=F`$ when $`p>0`$. The problem thus reduces to minimizing the function

$$E(F_z,N_s)=\frac{c_1}{2V^{\mathrm{eff}}}[F_z(F_z+1)6N]+\frac{c_2}{5V^{\mathrm{eff}}}N_\mathrm{S}(2N2N_\mathrm{S}+3)pF_z.$$
(139)
For this purpose, it is convenient to consider the cases $`c_2>0`$ and $`c_2<0`$ separately.

#### 1 Case of $`c_2>0`$

Let us first consider the case $`c_2>0`$.

When $`c_1<0`$, the energy (139) is minimized when $`N_\mathrm{S}=0`$, $`N_0=N`$, and $`F=F_z=2N`$, and the ground state is given by $`(\widehat{a}_2^{})^N|\mathrm{vac}`$, that is, the system is ferromagnetic. This result is consistent with that obtained from MFT.

When $`c_1>0`$, let us rewrite the energy as

$$E(F_z,N_s)=\frac{c_1}{2V^{\mathrm{eff}}}\left(F_z\frac{pV^{\mathrm{eff}}}{c_1}+\frac{1}{2}\right)^2\frac{2c_2}{5V^{\mathrm{eff}}}\left(N_\mathrm{S}\frac{2N+3}{4}\right)^2+\mathrm{const}.$$
(140)
The energy is thus lower when $`F_z`$ is closer to $`pV^{\mathrm{eff}}/c_11/2`$ and when $`N_\mathrm{S}`$ is smaller. In most of the parameter space, the ground state is $`|N_\mathrm{S}=0,N_0=N,F,F_z=F;\lambda `$ with $`F`$ taking the allowed integer closest to $`pV^{\mathrm{eff}}/c_11/2`$. Since these states satisfy $`\widehat{S}_+\widehat{S}_{}=0`$, they belong to the ferromagnetic phase or the cyclic phase. The separatrix between the two phases is given by $`pV^{\mathrm{eff}}/c_11/2=2N1`$. We thus find that the ground state is ferromagnetic if $`c_1<pV^{\mathrm{eff}}/(2N1/2)pV^{\mathrm{eff}}/(2N)`$ and cyclic otherwise. This classification is consistent with the mean-field analysis given in Sec. III D. As seen in Sec. IV A, the above ground state is highly degenerated in general; this may originate from the continuous symmetry that leaves free at least two parameters characterizing the order parameter of the cyclic phase as shown in Sec. III C. According to the discussions in Sec. IV A, the degeneracy of the states $`|N_\mathrm{S}=0,N_0=N,F,F_z=F;\lambda `$ is equal to the number of the combinations of $`\{n_{12},n_{22},n_{30}\}`$ satisfying $`n_{12}+2n_{22}+3n_{30}=N3n_{33}`$ and $`2n_{12}+2n_{22}=F3n_{33}`$. The number of trios, $`n_{30}`$, can take values in the range $`(NF)/3n_{30}(NF/2)/3`$. When the magnetic field is nearly zero and $`F0`$, there is little degeneracy and $`n_{30}N/3`$, namely, almost all bosons form trios.

When $`pV^{\mathrm{eff}}/c_11/2`$ is close to the forbidden values of $`F`$ ($`F=2,5`$ if $`N_0=0`$ (mod 3), and $`F=0,3`$ otherwise) for the above state with $`N_\mathrm{S}=0`$, $`F_z`$ may take those values at the cost of increasing $`N_\mathrm{S}`$ to 1 or 2, since any of the three values (0,1,2) of $`N_0`$ mod 3 is realized by setting $`N_\mathrm{S}`$ as 0, 1, or 2 \[recall the relation $`N=2N_\mathrm{S}+N_0`$\]. Whether or not the states $`|N_\mathrm{S}=1,N_0=N2,F,F_z=F;\lambda `$ and $`|N_\mathrm{S}=2,N_0=N4,F,F_z=F;\lambda `$ can be the lowest-energy state depends on the ratio $`c_2/c_1`$. In Fig 1, we give diagrams of the ground states for small $`pV^{\mathrm{eff}}/c_1`$.

#### 2 Case of $`c_2<0`$

In this case, it is convenient to introduce a new parameter

$$c_1^{}c_1\frac{c_2}{20},$$
(141)
and consider the energy as a function of $`F_z`$ and $`l2N_0F_z`$:

$$E(F_z,l)=\frac{c_1^{}}{2V^{\mathrm{eff}}}\left[F_z\frac{V^{\mathrm{eff}}}{c_1^{}}\left(p+\frac{c_2}{8V^{\mathrm{eff}}}\right)+\frac{1}{2}\right]^2\frac{c_2}{40V^{\mathrm{eff}}}l(l+2F+6)+\mathrm{const}.$$
(142)
Since $`c_2<0`$, we see that $`E(F_z,l)`$ is an increasing function of $`l`$, namely,

$$E(F_z,l)>E(F_z,l^{})\mathrm{for}l>l^{}.$$
(143)
Let us consider the cases $`c_1^{}<0`$ and $`c_1^{}>0`$ separately.

(a) $`c_1^{}<0`$ — In this parameter region, MFT predicts that the system is ferromagnetic, namely, the lowest-energy state always shows $`F_z=2N`$ regardless of the magnitude of magnetic field $`p`$. In the exact solution considered here, the magnetic response is different because of the offset term $`c_2/(8V^{\mathrm{eff}})`$ in Eq. (142). Since $`c_2<0`$, this term counteracts the applied magnetic field. It is thus expected that magnetization is suppressed when magnetic field is weak. The exact ground state is derived as follows. When $`N`$ is even, the state $`(\widehat{A}_0^{(2)})^{N/2}|\mathrm{vac}`$ has energy $`E(F_z=0,l=0)`$, and the state $`(\widehat{a}_2^{})^N|\mathrm{vac}`$ has energy $`E(F_z=2N,l=0)`$. Any other set $`\{F_z,l\}`$ gives an energy higher than one of these states. $`(\widehat{A}_0^{(2)})^{N/2}|\mathrm{vac}`$ is thus the ground state when $`E(F_z=0,l=0)<E(F_z=2N,l=0)`$, or equivalently,

$$pV^{\mathrm{eff}}<c_1^{}\left(N+\frac{1}{2}\right)\frac{c_2}{8}$$
(144)
and otherwise the ground state is $`(\widehat{a}_2^{})^N|\mathrm{vac}`$ \[see Fig 2\]. When $`N`$ is odd, $`F_z=0`$ is attained only when $`N_03`$, $`F_z=1`$ is forbidden, and the state $`\widehat{a}_2^{}(\widehat{A}_0^{(2)})^{(N1)/2}|\mathrm{vac}`$ has energy $`E(F_z=2,l=0)`$. It is easy to confirm that $`E(F_z=2,l=0)<E(F_z=0,l=6)`$ always holds. Therefore, $`\widehat{a}_2^{}(\widehat{A}_0^{(2)})^{(N1)/2}|\mathrm{vac}`$ is the ground state when $`E(F_z=2,l=0)<E(F_z=2N,l=0)`$, or equivalently,

$$pV^{\mathrm{eff}}<c_1^{}\left(N+\frac{3}{2}\right)\frac{c_2}{8},$$
(145)
and otherwise the ground state is $`(\widehat{a}_2^{})^N|\mathrm{vac}`$ \[see Fig 2\]. These results indicate that in the parameter region of $`c_28Nc_1^{}<0`$, magnetization of the ground state jumps from 0 or 2 to $`2N`$. Such a huge discontinuity does not appear in MFT with a linear Zeeman effect. (However, in the presence of a quadratic Zeeman effect, such a jump occurs also in MFT .)

(b) $`c_1^{}>0`$ — Given $`F_z6`$, the minimum allowed value of $`l`$ is determined as follows. Note that $`l=2NF_z4N_\mathrm{S}`$ is minimized when the number of singlets $`N_\mathrm{S}`$ is maximized. For $`F_z=2N4k`$ ($`k`$ is an integer), the state $`(\widehat{A}_0^{(2)})^k(\widehat{a}_2^{})^{N2k}|\mathrm{vac}`$ gives $`l=0`$. To increase $`F_z`$ by one ($`F_z=2N4k+1`$), one singlet pair must be broken and the minimum of $`l`$ is $`l=3`$ given by the state $`(\widehat{A}_0^{(2)})^{k1}(\widehat{a}_2^{})^{N2k1}\widehat{A}_3^{(3)}|\mathrm{vac}`$. Keeping the singlet part, $`F_z`$ can be further increased to $`F_z=2N4k+2`$ by the state $`(\widehat{A}_0^{(2)})^{k1}(\widehat{a}_2^{})^{N2k}\widehat{A}_2^{(2)}|\mathrm{vac}`$ with $`l=2`$. Since $`F_z=F=2N_01`$ is forbidden, $`F_z=2N4k+3`$ requires one more singlet pair to break up, resulting in $`l=5`$ with the state $`(\widehat{A}_0^{(2)})^{k2}(\widehat{a}_2^{})^{N2k1}\widehat{A}_2^{(2)}\widehat{A}_3^{(3)}|\mathrm{vac}`$. When $`(pV^{\mathrm{eff}}+c_2/8)/c_1^{}1/2`$ falls between $`2N4k`$ and $`2N4(k1)`$, the lowest energy is the minimum of $`E(F_z=2N4k,l=0)`$, $`E(F_z=2N4k+1,l=3)`$, $`E(F_z=2N4k+2,l=3)`$, $`E(F_z=2N4k+3,l=5)`$, and $`E(F_z=2N4k+4,l=0)`$. From Eq. (142), we expect that when $`|c_2|/c_1^{}`$ is large, nonzero $`l`$ pushes up the energy significantly and cannot be the ground state, so that $`F_z`$ increases stepwise with $`\mathrm{\Delta }F_z=4`$. When $`|c_2|/c_1^{}`$ is small, $`F_z`$ will increase with the step size of $`\mathrm{\Delta }F_z=1`$. This is confirmed by explicitly calculating $`E(F_z,l)`$ using Eq. (142), and we obtain the diagrams in Fig. 3.

In the region $`pV^{\mathrm{eff}}/c_1^{}<f_1(|c_2|/c_1^{})`$ with $`f_1(x)=(x20)(9x4)(32x)^1`$ (see broken curves in Fig. 3), $`F_z`$ increases by taking every integer. When $`pV^{\mathrm{eff}}/c_1^{}>f_1(|c_2|/c_1^{})`$, the values of $`F_z=2N4k+3`$ are suppressed. In the region $`pV^{\mathrm{eff}}/c_1^{}>f_2(|c_2|/c_1^{})`$ with $`f_2(x)=(x20)(13x20)(80x)^1`$, the values $`F_z=2N4k+1`$ are suppressed, and when $`pV^{\mathrm{eff}}/c_1^{}>f_3(|c_2|/c_1^{})`$ with $`f_3(x)=(x20)(x8)(8x)^1`$, the values $`F_z=2N4k+2`$ are further suppressed and the step size becomes 4.

While the averaged slope $`\mathrm{\Delta }F_z/\mathrm{\Delta }pV^{\mathrm{eff}}/c_1^{}`$ coincides with that in MFT, the offset term $`c_2/(8V^{\mathrm{eff}})`$ in Eq. (142) (see also the broken lines in Fig. 3) makes a qualitative distinction from MFT, namely, the onset of the magnetization displaces from $`p=0`$ to $`p=|c_2|/(8V^{\mathrm{eff}})`$. Note that the slope $`V^{\mathrm{eff}}/c_1^{}`$ and the offset $`|c_2|/(8V^{\mathrm{eff}})`$ are determined by independent parameters. A typical behavior of the magnetic response when $`|c_2|c_1^{}`$ is shown in Fig. 4.

### D Property of ground states for $`c_2<0`$

We now calculate the Zeeman-level populations of the ground states for $`c_2<0`$. In MFT, the lowest-energy states for $`c_2<0`$ have vanishing population in the $`m=0,\pm 1`$ levels. In contrast, the exact ground states derived in the preceding subsection, $`(\widehat{A}_0^{(2)})^{N_\mathrm{S}}(\widehat{a}_2^{})^{n_{12}}(\widehat{A}_2^{(2)})^{n_{22}}(\widehat{A}_3^{(3)})^{n_{33}}|\mathrm{vac}`$ with $`n_{22}=0,1`$ and $`n_{33}=0,1`$, have nonzero populations in the $`m=0,\pm 1`$ levels. The exact forms for the averaged population $`\widehat{a}_m^{}\widehat{a}_m`$ are calculated as follows. The above ground states have the form of $`(\widehat{A}_0^{(2)})^{N_\mathrm{S}}|\varphi (\widehat{𝒮}_+)^{N_\mathrm{S}}|\varphi `$ with $`|\varphi `$ being a state with a fixed number ($`sn_{12}+2n_{22}+3n_{33}`$) of bosons satisfying $`\widehat{𝒮}_{}|\varphi =0`$. The average Zeeman population for the ground states,

$$\widehat{a}_m^{}\widehat{a}_m\frac{\varphi |(\widehat{𝒮}_{})^{N_\mathrm{S}}\widehat{a}_m^{}\widehat{a}_m(\widehat{𝒮}_+)^{N_\mathrm{S}}|\varphi }{\varphi |(\widehat{𝒮}_{})^{N_\mathrm{S}}(\widehat{𝒮}_+)^{N_\mathrm{S}}|\varphi },$$
(146)
is then simply related to the average Zeeman populations for the state $`|\varphi `$ as

$$\widehat{a}_m^{}\widehat{a}_m=\widehat{a}_m^{}\widehat{a}_m_0+\frac{N_\mathrm{S}}{s+5/2}(\widehat{a}_m^{}\widehat{a}_m_0+\widehat{a}_m^{}\widehat{a}_m_0+1),$$
(147)
where $`\widehat{a}_m^{}\widehat{a}_m_0\varphi |\widehat{a}_m^{}\widehat{a}_m|\varphi /\varphi |\varphi `$. The derivation of the formula (147) is given in Appendix B. The formula implies that when $`N_\mathrm{S}s`$, the Zeeman populations of the ground states is sensitive to the form of $`|\varphi `$.

With this formula, it is a straightforward task to calculate average Zeeman-level populations for four types of ground states, $`(\widehat{A}_0^{(2)})^{N_\mathrm{S}}(\widehat{a}_2^{})^{n_{12}}(\widehat{A}_2^{(2)})^{n_{22}}(\widehat{A}_3^{(3)})^{n_{33}}|\mathrm{vac}`$ with $`n_{22}=0,1`$ and $`n_{33}=0,1`$. The exact result will be given in Appendix B. A striking feature appears in the leading terms under the condition $`1n_{12}N_\mathrm{S}`$. The results are summarized as

$$\widehat{a}_1^{}\widehat{a}_1\widehat{a}_1^{}\widehat{a}_1N_\mathrm{S}(1+n_{33})/n_{12}$$
(148)
and

$$\widehat{a}_0^{}\widehat{a}_0N_\mathrm{S}(1+2n_{22})/n_{12}.$$
(149)
As seen in the preceding subsection, with the increase of magnetic field, the ground state alternates among the four types of states. While this change causes a very small difference in magnetization, it leads to large changes in the average Zeeman-level populations, by a factor of 2 or 3. The origin of this drastic change may be explained as follows. Let us first consider the state $`(\widehat{a}_2^{})^{n_{12}}|\mathrm{vac}`$ with $`n_{12}0`$. This state has no population in Zeeman levels $`m=0,\pm 1`$. When $`\widehat{A}_2^{(2)}`$ is applied to this state, the operator $`\widehat{a}_2^{}`$ that appears in $`\widehat{A}_2^{(2)}`$ has effectively a large amplitude of the order of $`\sqrt{n_{21}}`$. Hence the term $`\widehat{a}_2^{}\widehat{a}_0^{}`$ is dominant, and it approximately adds one boson to the $`m=2`$ level and one boson to the $`m=0`$ level. Hence the $`m=0`$ population of the state $`\widehat{A}_2^{(2)}(\widehat{a}_2^{})^{n_{12}}|\mathrm{vac}`$ is close to unity. This change is then amplified by a factor of $`N_\mathrm{S}/n_{12}`$ according to the formula (147), leading to Eq. (149). Similarly, applying $`\widehat{A}_3^{(3)}`$ effectively results in adding of two bosons to the $`m=2`$ level and one boson to the $`m=1`$ level through the dominant term $`(\widehat{a}_2^{})^2\widehat{a}_1^{}`$.

## V Low-lying excitation spectra

In this section we study the low-lying excitation spectrum of spin-2 BECs in the thermodynamic limit using the Bogoliubov approximation. We shall see that the symmetry of each ground state discussed in Sec. III is reflected in the excitation spectrum.

### A Effective Hamiltonian

In the center-of-mass frame of the system BEC occurs in the $`𝒌=\mathrm{𝟎}`$ state. We therefore decompose the operators that appear in Eq. (35) into the $`𝒌=\mathrm{𝟎}`$ components and the $`𝐤\mathrm{𝟎}`$ ones. The first term on the right-hand side (rhs) of Eq. (41) is rewritten as

$`{\displaystyle \underset{𝒌}{}}:\widehat{\rho }_𝒌^{}\widehat{\rho }_𝒌:=:\widehat{\rho }_\mathrm{𝟎}^{}\widehat{\rho }_\mathrm{𝟎}:+{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}:\widehat{\rho }_𝒌^{}\widehat{\rho }_𝒌:=:\widehat{N}^2:+{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}:\widehat{\rho }_𝒌^{}\widehat{\rho }_𝒌:.`$ (150)
If we ignore the terms that do not include the $`𝒌=\mathrm{𝟎}`$ components, we may approximate $`\widehat{\rho }_{𝒌\mathrm{𝟎}}`$ as

$`\widehat{\rho }_{𝒌\mathrm{𝟎}}{\displaystyle \underset{m}{}}(\widehat{a}_{\mathrm{𝟎},m}^{}\widehat{a}_{𝒌,m}+\widehat{a}_{𝒌,m}^{}\widehat{a}_{\mathrm{𝟎},m})=\widehat{D}_𝒌+\widehat{D}_𝒌^{},`$ (151)
where

$`\widehat{D}_𝒌{\displaystyle \underset{m}{}}\widehat{a}_{\mathrm{𝟎},m}^{}\widehat{a}_{𝒌,m}.`$ (152)
We thus obtain

$`{\displaystyle \underset{𝒌}{}}:\widehat{\rho }_𝒌^{}\widehat{\rho }_𝒌::\widehat{N}^2:+{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}:(2\widehat{D}_𝒌^{}\widehat{D}_𝒌+\widehat{D}_𝒌\widehat{D}_𝐤+\widehat{D}_𝒌^{}\widehat{D}_𝐤^{}):.`$ (153)
Similarly, we may approximate the second term on the rhs of Eq. (41) as

$`{\displaystyle \underset{𝒌}{}}:\widehat{𝑭}_𝒌^{}\widehat{𝑭}_𝒌:`$ $``$ $`:\widehat{𝑭}^2:+2{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}{\displaystyle \underset{mn}{}}\widehat{𝑭}𝒇_{mn}\widehat{a}_{𝒌,m}^{}\widehat{a}_{𝒌,n}+{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}{\displaystyle \underset{ijmn}{}}𝒇_{ij}𝒇_{mn}`$ (155)
$`\times (2\widehat{a}_{\mathrm{𝟎},i}^{}\widehat{a}_{\mathrm{𝟎},n}\widehat{a}_{𝒌,m}^{}\widehat{a}_{𝒌,j}+\widehat{a}_{\mathrm{𝟎},i}^{}\widehat{a}_{\mathrm{𝟎},m}^{}\widehat{a}_{𝒌,j}\widehat{a}_{𝒌,n}+\widehat{a}_{\mathrm{𝟎},j}\widehat{a}_{\mathrm{𝟎},n}\widehat{a}_{𝒌,i}^{}\widehat{a}_{𝒌,m}^{}),`$
where

$`\widehat{𝑭}={\displaystyle \underset{mn}{}}𝒇_{mn}\widehat{a}_{\mathrm{𝟎},m}^{}\widehat{a}_{\mathrm{𝟎},n}.`$ (156)
The third term on the rhs of Eq. (41) is decomposed into

$`{\displaystyle \underset{𝒌}{}}\widehat{A}_𝒌^{}\widehat{A}_𝒌{\displaystyle \frac{4}{5}}\widehat{S}_+\widehat{S}_{}+{\displaystyle \frac{2}{5}}(\widehat{S}_+{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}{\displaystyle \underset{m}{}}(1)^m\widehat{a}_{𝒌,m}\widehat{a}_{𝐤,m}+\mathrm{H}.\mathrm{c}.)+{\displaystyle \frac{4}{5}}{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}{\displaystyle \underset{mn}{}}(1)^{m+n}\widehat{a}_{\mathrm{𝟎},m}^{}\widehat{a}_{\mathrm{𝟎},n}\widehat{a}_{𝒌,m}^{}\widehat{a}_{𝒌,n},`$ (157)
where

$`\widehat{S}_{}=\widehat{S}_+^{}={\displaystyle \frac{1}{2}}{\displaystyle \underset{m}{}}(1)^m\widehat{a}_{\mathrm{𝟎},m}\widehat{a}_{𝒌,m},`$ (158)
and H.c. denotes Hermitian conjugates of the preceding terms. Substituting Eqs. (151), (155) and (157) into Eq. (41), we obtain

$`\widehat{V}`$ $``$ $`{\displaystyle \frac{c_0}{2V}}:\widehat{N}^2:+{\displaystyle \frac{c_1}{2V}}:\widehat{𝑭}^2:+{\displaystyle \frac{2c_2}{5V}}\widehat{S}_+\widehat{S}_{}+{\displaystyle \frac{c_0}{2V}}{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}:(2\widehat{D}_𝒌^{}\widehat{D}_𝒌+\widehat{D}_𝒌\widehat{D}_𝒌+\widehat{D}_𝒌^{}\widehat{D}_𝐤^{}):`$ (161)
$`+{\displaystyle \frac{c_1}{V}}{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}{\displaystyle \underset{mn}{}}\widehat{𝑭}𝒇_{mn}\widehat{a}_{𝒌,m}^{}\widehat{a}_{𝒌,n}+{\displaystyle \frac{c_1}{2V}}{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}{\displaystyle \underset{ijmn}{}}𝒇_{ij}𝒇_{mn}(2\widehat{a}_{\mathrm{𝟎},i}^{}\widehat{a}_{\mathrm{𝟎},n}\widehat{a}_{𝒌,m}^{}\widehat{a}_{𝒌,j}+\widehat{a}_{\mathrm{𝟎},i}^{}\widehat{a}_{\mathrm{𝟎},m}^{}\widehat{a}_{𝒌,j}\widehat{a}_{𝒌,n}+\widehat{a}_{\mathrm{𝟎},j}\widehat{a}_{\mathrm{𝟎},n}\widehat{a}_{𝒌,i}^{}\widehat{a}_{𝒌,m}^{})`$
$`+{\displaystyle \frac{c_2}{5V}}(\widehat{S}_+{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}{\displaystyle \underset{m}{}}(1)^m\widehat{a}_{𝒌,m}\widehat{a}_{𝒌,m}+\mathrm{H}.\mathrm{c}.)+{\displaystyle \frac{2c_2}{5V}}{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}{\displaystyle \underset{mn}{}}(1)^{m+n}\widehat{a}_{\mathrm{𝟎},m}^{}\widehat{a}_{\mathrm{𝟎},n}\widehat{a}_{𝒌,m}^{}\widehat{a}_{𝒌,n}.`$
In the Bogoliubov approximation we replace operators $`\widehat{a}_{\mathrm{𝟎},m}`$ by c-numbers $`\sqrt{N^{\mathrm{BEC}}}\zeta _m`$, where $`\zeta _m`$’s denote the complex mean-field amplitudes introduced in Sec. III and $`N^{\mathrm{BEC}}`$ is the number of condensate bosons. Since $`N^{\mathrm{BEC}}`$ is smaller than the total number of bosons $`N`$ due to the interparticle interactions, we take into account the conservation of the total number of bosons through the relation

$`N^{\mathrm{BEC}}+{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}{\displaystyle \underset{m}{}}\widehat{n}_{𝒌,m}=N.`$ (162)
Then Eq. (40) becomes

$`\widehat{H}_0=pN\widehat{f}_z+{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}{\displaystyle \underset{m}{}}\left[ϵ_𝒌p(m\widehat{f}_z)\right]\widehat{n}_{𝒌,m},`$ (163)
and Eq. (161) becomes

$`\widehat{V}`$ $``$ $`N\left({\displaystyle \frac{c_0n}{2}}+{\displaystyle \frac{c_1n}{2}}\widehat{𝒇}^2+{\displaystyle \frac{2c_2n}{5}}|\widehat{s}_{}|^2\right)\left(c_1n\widehat{𝒇}^2+{\displaystyle \frac{4c_2n}{5}}|\widehat{s}_{}|^2\right){\displaystyle \underset{𝒌\mathrm{𝟎}}{}}{\displaystyle \underset{m}{}}\widehat{n}_{𝒌,m}`$ (167)
$`+{\displaystyle \frac{c_0n}{2}}{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}{\displaystyle \underset{mn}{}}(2\zeta _m^{}\zeta _n\widehat{a}_{𝒌,n}^{}\widehat{a}_{𝒌,m}+\zeta _m^{}\zeta _n^{}\widehat{a}_{𝒌,m}\widehat{a}_{𝒌,n}+\zeta _m\zeta _n\widehat{a}_{𝒌,m}^{}\widehat{a}_{𝒌,n})`$
$`+c_1n{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}{\displaystyle \underset{mn}{}}\widehat{𝒇}𝒇_{mn}\widehat{a}_{𝒌,m}^{}\widehat{a}_{𝒌,n}+{\displaystyle \frac{c_1n}{2}}{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}{\displaystyle \underset{ijmn}{}}𝒇_{ij}𝒇_{mn}(2\zeta _i^{}\zeta _n\widehat{a}_{𝒌,m}^{}\widehat{a}_{𝒌,j}+\zeta _i^{}\zeta _m^{}\widehat{a}_{𝒌,j}\widehat{a}_{𝒌,n}+\zeta _j\zeta _n\widehat{a}_{𝒌,i}^{}\widehat{a}_{𝒌,m}^{})`$
$`+{\displaystyle \frac{c_2n}{5}}(\widehat{s}_{}^{}{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}{\displaystyle \underset{m}{}}(1)^m\widehat{a}_{𝒌,m}\widehat{a}_{𝒌,m}+\mathrm{H}.\mathrm{c}.)+{\displaystyle \frac{2c_2n}{5}}{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}{\displaystyle \underset{mn}{}}(1)^{m+n}\zeta _m^{}\zeta _n\widehat{a}_{𝒌,m}^{}\widehat{a}_{𝒌,n},`$
where $`nN/V`$, and the definitions of $`\widehat{f}_z`$, $`\widehat{𝒇}`$, and $`\widehat{s}_{}`$ are given in Eqs. (48)–(50). Equations (163) and (167) constitute our basic Hamiltonian in the following discussions. We use this Hamiltonian to examine low-lying excitation spectra for each phase.

### B Excitation spectrum of a ferromagnetic BEC

Let us first examine the excitation spectrum of a ferromagnetic phase in which BEC occurs in the $`m=2`$ state. Then $`\zeta _m=2\delta _{m,2}`$ and Eqs. (163) and (167) become

$`\widehat{H}_0=2p\widehat{N}+{\displaystyle \underset{𝒌\mathrm{𝟎},m}{}}\left[ϵ_𝒌+(2m)p\right]\widehat{n}_{𝒌,m},`$ (168)
and

$`\widehat{V}`$ $``$ $`{\displaystyle \frac{1}{2}}g_4n\widehat{N}+n{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}\left[c_0\widehat{n}_{𝒌,2}+2c_1(2\widehat{n}_{𝒌,2}2\widehat{n}_{𝒌,0}3\widehat{n}_{𝒌,1}4\widehat{n}_{𝒌,2})+{\displaystyle \frac{2}{5}}c_2\widehat{n}_{𝒌,2}+{\displaystyle \frac{g_4}{2}}(\widehat{a}_{𝒌,2}\widehat{a}_{𝒌,2}+\widehat{a}_{𝒌,2}^{}\widehat{a}_{𝒌,2}^{})\right],`$ (169)
respectively, where $`g_4=c_0+4c_1`$. The total Hamiltonian is therefore given by

$`\widehat{H}^\mathrm{F}`$ $`=`$ $`{\displaystyle \frac{1}{2}}g_4n\widehat{N}2p\widehat{N}`$ (173)
$`+{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}\left[\left(ϵ_𝒌+g_4n\right)\widehat{n}_{𝒌,2}+{\displaystyle \frac{1}{2}}g_4n(\widehat{a}_{𝒌,2}\widehat{a}_{𝒌,2}+\widehat{a}_{𝒌,2}^{}\widehat{a}_{𝒌,2}^{})\right]`$
$`+{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}[(ϵ_𝒌+p)\widehat{n}_{𝒌,1}+(ϵ_𝒌+2p4c_1n)\widehat{n}_{𝒌,0}`$
$`+(ϵ_𝒌+3p6c_1n)\widehat{n}_{𝒌,1}+(ϵ_𝒌+4p8c_1n+2c_2n/5)\widehat{n}_{𝒌,2}].`$
The second line gives the Bogoliubov spectrum

$`E_{m=2,𝒌}^\mathrm{F}=\sqrt{ϵ_𝒌^2+2g_4nϵ_𝐤},`$ (174)
while other terms give single-particle spectra:

$`E_{m=1,𝒌}^\mathrm{F}`$ $`=`$ $`ϵ_𝒌+p`$ (175)
$`E_{m=0,𝒌}^\mathrm{F}`$ $`=`$ $`ϵ_𝒌+2p4c_1n`$ (176)
$`E_{m=1,𝒌}^\mathrm{F}`$ $`=`$ $`ϵ_𝒌+3p6c_1n`$ (177)
$`E_{m=2,𝒌}^\mathrm{F}`$ $`=`$ $`ϵ_𝒌+4p8c_1n+2c_2n/5`$ (178)
For the Bogoliubov excitation energy to be positive, we must have

$`g_4={\displaystyle \frac{4\pi \mathrm{}^2}{M}}a_4>0.`$ (179)
That is, the s-wave scattering length for the total spin-4 channel must be positive. This condition is the same as that required for the ferromagnetic mean field to be stable, that is, the first term on the rhs of Eq. (173) being positive. For the single-particle excitation energies to be positive, we must have $`p>2c_1n`$ and $`p>(2c_1c_2/10)n`$. These conditions are the same as those in (88) for which the ferromagnetic phase is the lowest-energy mean field (note that $`\stackrel{~}{p}2p/n`$).

We note that the Bogoliubov spectrum (174) is independent of applied magnetic field and remains massless in its presence. This Goldstone mode is a consequence of the global U(1) gauge invariance due to the conservation of the total number of bosons, as discussed in Sec. III A.

### C Excitation spectrum of an antiferromagnetic BEC

Let us next examine the excitation spectrum of an antiferromagnetic phase in which the order parameter is given by (64). Making the replacements

$`\zeta _2=e^{i\varphi _2}\sqrt{{\displaystyle \frac{1}{2}}+{\displaystyle \frac{\widehat{f}_z}{4}}},\zeta _2=e^{i\varphi _2}\sqrt{{\displaystyle \frac{1}{2}}{\displaystyle \frac{\widehat{f}_z}{4}}},\zeta _0=\zeta _{\pm 1}=0,`$ (180)
and substituting Eq. (180) into Eq. (167), together with Eq. (163), we obtain the total Hamiltonian of an antiferromagnetic BEC:

$`\widehat{H}^{\mathrm{AF}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}(c_0+c_2/5)nN+{\displaystyle \frac{1}{2}}(c_1c_2/20)nN\widehat{f}_z^2pN\widehat{f}_z`$ (188)
$`+{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}\{[ϵ_𝒌+p(\widehat{f}_z2)+c_0n|\zeta _2|^2+c_1n(4|\zeta _2|^2+2\widehat{f}_z\widehat{f}_z^2)+{\displaystyle \frac{2c_2n}{5}}(|\zeta _2|^22|\zeta _2\zeta _2|^2)]\widehat{n}_{𝒌,2}`$
$`+\left[ϵ_𝒌+p(\widehat{f}_z+2)+c_0n|\zeta _2|^2+c_1n(4|\zeta _2|^22\widehat{f}_z\widehat{f}_z^2)+{\displaystyle \frac{2c_2n}{5}}(|\zeta _2|^22|\zeta _2\zeta _2|^2)\right]\widehat{n}_{𝒌,2}`$
$`+[{\displaystyle \frac{1}{2}}g_4n(\zeta _2^2\widehat{a}_{𝒌,2}\widehat{a}_{𝒌,2}+\zeta _2^2\widehat{a}_{𝒌,2}\widehat{a}_{𝒌,2})+(c_04c_1+2c_2/5)n(\zeta _2^{}\zeta _2^{}\widehat{a}_{𝒌,2}\widehat{a}_{𝒌,2}+\zeta _2\zeta _2^{}\widehat{a}_{𝒌,2}^{}\widehat{a}_{𝒌,2})+\mathrm{H}.\mathrm{c}.]\}`$
$`+{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}\{[ϵ_𝒌+p(\widehat{f}_z1)+c_1n(2|\zeta _2|^2+\widehat{f}_z\widehat{f}_z^2){\displaystyle \frac{4c_2n}{5}}|\zeta _2\zeta _2|^2]\widehat{n}_{𝒌,1}`$
$`+\left[ϵ_𝒌+p(\widehat{f}_z+1)+c_1n(2|\zeta _2|^2\widehat{f}_z\widehat{f}_z^2){\displaystyle \frac{4c_2n}{5}}|\zeta _2\zeta _2|^2\right]\widehat{n}_{𝒌,1}`$
$`+2(c_1c_2/5)n(\zeta _2^{}\zeta _2^{}\widehat{a}_{𝒌,1}\widehat{a}_{𝒌,1}+\mathrm{H}.\mathrm{c}.)\}`$
$`+{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}\{[ϵ_𝐤+p\widehat{f}_zc_1n\widehat{f}_z^2{\displaystyle \frac{4c_2n}{5}}|\zeta _2\zeta _2|^2]\widehat{n}_{𝒌,0}+{\displaystyle \frac{c_2n}{5}}(\zeta _2^{}\zeta _2^{}\widehat{a}_{𝒌,0}\widehat{a}_{𝒌,0}+\mathrm{H}.\mathrm{c}.)\}.`$
This Hamiltonian may be simplified using the relation $`p\widehat{f}_z(c_1c_2/20)n`$, giving

$`\widehat{H}^{\mathrm{AF}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}(c_0+c_2/5)nN{\displaystyle \frac{1}{2}}pN\widehat{f}_z`$ (195)
$`+{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}\{(ϵ_𝒌+g_4n|\zeta _2|^2)\widehat{n}_{𝒌,2}+(ϵ_𝒌+g_4n|\zeta _2|^2)\widehat{n}_{𝒌,2}`$
$`+[{\displaystyle \frac{g_4n}{2}}(\zeta _2^2\widehat{a}_{𝒌,2}\widehat{a}_{𝒌,2}+\zeta _2^2\widehat{a}_{𝒌,2}\widehat{a}_{𝒌,2})+(c_04c_1+2c_2/5)n(\zeta _2^{}\zeta _2^{}\widehat{a}_{𝒌,2}\widehat{a}_{𝒌,2}+\zeta _2\zeta _2^{}\widehat{a}_{𝒌,2}^{}\widehat{a}_{𝒌,2})+\mathrm{H}.\mathrm{c}.]\}`$
$`+{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}\{[ϵ_𝒌+(c_1c_2/5)n+{\displaystyle \frac{1}{2}}(c_1+c_2/10)n\widehat{f}_z]\widehat{n}_{𝒌,1}+[ϵ_𝒌+(c_1c_2/5)n{\displaystyle \frac{1}{2}}(c_1+c_2/10)n\widehat{f}_z]\widehat{n}_{𝒌,1}`$
$`+2(c_1c_2/5)n(\zeta _2^{}\zeta _2^{}\widehat{a}_{𝒌,1}\widehat{a}_{𝒌,1}+\mathrm{H}.\mathrm{c}.)\}`$
$`+{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}\{(ϵ_𝒌c_2n/5)\widehat{n}_{𝒌,0}+{\displaystyle \frac{c_2n}{5}}(\zeta _2^{}\zeta _2^{}\widehat{a}_{𝒌,0}\widehat{a}_{𝒌,0}+\mathrm{H}.\mathrm{c}.)\}.`$
This result shows that the eigenmodes are classified into three categories: the $`m=0`$ mode, the coupled $`m=\pm 1`$ modes, and the coupled $`m=\pm 2`$ modes. Below we analyze each of them.

#### 1 The spin-0 quasiparticles

The Hamiltonian (195) shows that the $`m=0`$ mode is decoupled from other modes even in the presence of magnetic field. This part of the Hamiltonian can readily be diagonalized to give

$`E_{𝒌,0}=\sqrt{ϵ_𝒌^2+{\displaystyle \frac{2|c_2|n}{5}}ϵ_𝒌+\left({\displaystyle \frac{c_2n}{10}}\right)^2\widehat{f}_z^2}.`$ (196)
We note that the spectrum (196) becomes massive in the presence of magnetic field.

#### 2 The spin-1 quasiparticles

The $`m=1`$ and $`m=1`$ modes are coupled in the Hamiltonian (195), and the eigenenergies can be obtained by diagonalizing the following Hamiltonian:

$`\widehat{h}_1^{\mathrm{AF}}={\displaystyle \underset{𝒌\mathrm{𝟎}}{}}\{ϵ_{𝒌,1}\widehat{n}_{𝒌,1}+ϵ_{𝒌,1}\widehat{n}_{𝒌,1}+(\delta \widehat{a}_{𝒌,1}\widehat{a}_{𝒌,1}+\mathrm{H}.\mathrm{c}.)\},`$ (197)
where $`\delta 2(c_1c_2/5)n\zeta _2^{}\zeta _2^{}`$ and

$`ϵ_{𝒌,\pm 1}ϵ_𝒌+(c_1c_2/5)n\pm {\displaystyle \frac{1}{2}}(c_1+c_2/10)n\widehat{f}_z`$ (198)
The dispersion relation can be found by writing down the equations of motion for $`\widehat{a}_{𝒌,1}`$ and $`\widehat{a}_{𝒌,1}^{}`$ and seeking for the solution of the form $`\mathrm{exp}(iE_{𝒌,\pm 1}t/\mathrm{})`$, with the result

$`E_{𝒌,\pm 1}`$ $`=`$ $`\pm {\displaystyle \frac{1}{2}}\left(c_1+{\displaystyle \frac{c_2}{10}}\right)n\widehat{f}_z+\sqrt{ϵ_𝒌^2+2\left(c_1{\displaystyle \frac{c_2}{5}}\right)nϵ_𝐤+{\displaystyle \frac{1}{4}}\left(c_1{\displaystyle \frac{c_2}{5}}\right)^2n^2\widehat{f}_z^2}.`$ (199)
The excitation spectra (199) become massive in the presence of magnetic field. The positivity of this energy is guaranteed by the conditions $`c_1c_2/20>p/2n>0`$ and $`c_2<0`$ which are required for the antiferromagnetic phase to be the lowest-lying mean field (see condition (89)).

#### 3 The spin-2 quasiparticles

The $`m=2`$ and $`m=2`$ modes are coupled in the Hamiltonian (195), and the relevant part of the Hamiltonian reads

$`\widehat{h}_2`$ $`=`$ $`{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}\{ϵ_{𝒌,2}\widehat{n}_{𝒌,2}+ϵ_{𝒌,2}\widehat{n}_{𝒌,2}+{\displaystyle \frac{g_4n}{2}}(\zeta _2^2\widehat{a}_{𝒌,2}\widehat{a}_{𝒌,2}+\zeta _2^2\widehat{a}_{𝒌,2}\widehat{a}_{𝒌,2}+\mathrm{H}.\mathrm{c}.)`$ (201)
$`+(c_04c_1+2c_2/5)n(\zeta _2^{}\zeta _2^{}\widehat{a}_{𝒌,2}\widehat{a}_{𝒌,2}+\zeta _2\zeta _2^{}\widehat{a}_{𝒌,2}^{}\widehat{a}_{𝒌,2}+\mathrm{H}.\mathrm{c}.)\},`$
where

$`ϵ_{𝒌,2}=ϵ_𝒌+g_4n|\zeta _2|^2,ϵ_{𝒌,2}=ϵ_𝒌+g_4n|\zeta _2|^2.`$ (202)
By unitary transformations $`\widehat{a}_{\pm 𝒌,2}\widehat{a}_{\pm 𝒌,2}e^{i\varphi _2}`$ and $`\widehat{a}_{\pm 𝒌,2}\widehat{a}_{\pm 𝒌,2}e^{i\varphi _2}`$, Eq. (201) reduces to

$`\widehat{h}_2^{\mathrm{AF}}`$ $`=`$ $`{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}\{ϵ_{𝒌,2}\widehat{n}_{𝒌,2}+ϵ_{𝒌,2}\widehat{n}_{𝒌,2}+{\displaystyle \frac{1}{2}}(\alpha \widehat{a}_{𝒌,2}\widehat{a}_{𝒌,2}+\beta \widehat{a}_{𝒌,2}\widehat{a}_{𝒌,2}+\mathrm{H}.\mathrm{c}.)+\gamma (\widehat{a}_{𝒌,2}\widehat{a}_{𝒌,2}+\widehat{a}_{𝒌,2}^{}\widehat{a}_{𝒌,2}+\mathrm{H}.\mathrm{c}.)\},`$ (203)
where $`\alpha g_4n|\zeta _2|^2`$, $`\beta g_4n|\zeta _2|^2`$, and $`\gamma (c_04c_1+2c_2/5)n|\zeta _2\zeta _2|`$. This Hamiltonian can be diagonalized by writing down the Heisenberg equations of motion for $`\widehat{a}_{𝒌,2}+\widehat{a}_{𝒌,2}^{}`$ and $`\widehat{a}_{𝒌,2}+\widehat{a}_{𝒌,2}^{}`$ as

$`(i\mathrm{})^2{\displaystyle \frac{d^2}{dt^2}}(\widehat{a}_{𝒌,2}+\widehat{a}_{𝒌,2}^{})=ϵ_𝒌(ϵ_𝒌+2\alpha )(\widehat{a}_{𝒌,2}+\widehat{a}_{𝒌,2}^{})+2\gamma ϵ_𝒌(\widehat{a}_{𝒌,2}+\widehat{a}_{𝒌,2}^{}),`$ (204)
$`(i\mathrm{})^2{\displaystyle \frac{d^2}{dt^2}}(\widehat{a}_{𝒌,2}+\widehat{a}_{𝐤,2}^{})=ϵ_𝒌(ϵ_𝒌+\beta )(\widehat{a}_{𝒌,2}+\widehat{a}_{𝒌,2}^{})+2\gamma ϵ_𝒌(\widehat{a}_{𝒌,2}+\widehat{a}_{𝒌,2}^{}).`$ (205)
By assuming that $`\widehat{a}_{𝒌,\pm 2}+\widehat{a}_{𝒌,\pm 2}^{}\mathrm{exp}(iE_{𝒌,\pm 2}^{\mathrm{AF}}t/\mathrm{})`$, we obtain the following dispersion relations:

$`(E_{𝒌,\pm 2}^{\mathrm{AF}})^2=ϵ_𝒌\left[ϵ_𝒌+g_4n\pm g_4n\sqrt{{\displaystyle \frac{\widehat{f}_z^2}{4}}+\left[1{\displaystyle \frac{8}{g_4}}\left(c_1{\displaystyle \frac{c_2}{20}}\right)\right]^2\left(1{\displaystyle \frac{\widehat{f}_z^2}{4}}\right)}\right].`$ (206)
The positivity of this energy is met if the conditions $`c_1c_2/20>0`$ and $`c_0+c_2/5>0`$ are satisfied. The former condition is met whenever the antiferromagnetic phase is the lowest-energy state (see (89)), while the latter condition is required for the antiferromagnetic phase to be mechanically stable, that is, the first term on the rhs of Eq. (195) is positive. We note that the dispersion relations (206) are massless even in the presence of the magnetic field. They are the Goldstone modes associated with the U(1) gauge symmetry and the relative gauge symmetry (the rotational symmetry about the direction of the applied magnetic field) that are manifest in the mean-field solution discussed in Sec. III B.

At zero magnetic field, Eqs. (196), (199) and (206) reduces to

$`(E_{𝒌,0}^{\mathrm{AF}})^2`$ $`=`$ $`ϵ_𝒌\left[ϵ_𝒌(2c_2/5)n\right],`$ (207)
$`(E_{𝒌,\pm 1}^{\mathrm{AF}})^2`$ $`=`$ $`ϵ_𝒌\left[ϵ_𝒌+2(c_1c_2/5)n\right],`$ (208)
$`(E_{𝒌,\pm 2}^{\mathrm{AF}})^2`$ $`=`$ $`\{\begin{array}{cc}ϵ_𝒌\left[ϵ_𝒌+2(c_0+c_2/5)n\right],\hfill & \\ ϵ_𝒌\left[ϵ_𝒌+8(c_1c_2/20)n\right],\hfill & \end{array}`$ (209)
implying that all the five excitations are Goldstone modes. This reflects the fact that in the absence of the magnetic field the ground state is degenerate with respect to five continuous variables (see Eq. (67)).

### D Excitation spectrum of a cyclic BEC

We consider the case of Eq. (81), namely,

$`\zeta _{\pm 2}={\displaystyle \frac{1}{2}}\left(1\pm a\right)e^{i\varphi _{\pm 2}},\zeta _{\pm 1}=0,\zeta _0=\sqrt{{\displaystyle \frac{1a^2}{2}}}e^{i\varphi _0},`$ (210)
where $`a\widehat{f}_z/2`$ and $`\varphi _2+\varphi _22\varphi _0=\pm \pi `$. Because $`\widehat{s}_{}=0`$ and $`\widehat{f}_\pm =0`$, the interaction Hamiltonian (167) reduces to

$`\widehat{V}`$ $``$ $`{\displaystyle \frac{c_0n}{2}}+{\displaystyle \frac{1}{2}}pN\widehat{f}_zc_1n\widehat{f}_z^2{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}{\displaystyle \underset{m}{}}\widehat{n}_{𝒌,m}+{\displaystyle \frac{c_0}{2V}}{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}:(\widehat{D}_𝒌\widehat{D}_𝒌+\widehat{D}_𝒌^{}\widehat{D}_𝒌+\mathrm{H}.\mathrm{c}.):+c_1n\widehat{f}_z{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}{\displaystyle \underset{m}{}}m\widehat{n}_{𝒌,m}`$ (213)
$`+{\displaystyle \frac{c_1}{2V}}{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}{\displaystyle \underset{ijmn}{}}𝒇_{ij}𝒇_{mn}(\widehat{a}_{\mathrm{𝟎},i}^{}\widehat{a}_{\mathrm{𝟎},n}\widehat{a}_{𝒌,m}^{}\widehat{a}_{𝒌,j}+\widehat{a}_{\mathrm{𝟎},i}^{}\widehat{a}_{\mathrm{𝟎},m}^{}\widehat{a}_{𝒌,j}\widehat{a}_{𝒌,n}+\mathrm{H}.\mathrm{c}.)`$
$`+{\displaystyle \frac{2c_2}{5V}}{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}{\displaystyle \underset{mn}{}}(1)^{m+n}\widehat{a}_{\mathrm{𝟎},m}^{}\widehat{a}_{\mathrm{𝟎},n}\widehat{a}_{𝒌,m}^{}\widehat{a}_{𝒌,n}.`$
Substituting Eq. (210) into this, performing unitary transformations $`\widehat{a}_{\pm 𝒌,m}\widehat{a}_{\pm 𝒌,m}e^{i\varphi _m}`$, where $`\zeta _m=|\zeta _m|e^{i\varphi _m}`$ and $`\varphi _2+\varphi _22\varphi _0=\pi `$, and combining the result with Eq. (163), we obtain

$`\widehat{H}^\mathrm{C}`$ $`=`$ $`{\displaystyle \frac{1}{2}}c_0nN{\displaystyle \frac{1}{2}}pN\widehat{f}_z+{\displaystyle \underset{𝒌\mathrm{𝟎},m}{}}A_{𝒌,m}\widehat{n}_{𝒌,m}`$ (218)
$`+{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}\{{\displaystyle \frac{1}{2}}(\alpha +\beta )(|\zeta _2|^2\widehat{a}_{𝒌,2}\widehat{a}_{𝒌,2}+|\zeta _2|^2\widehat{a}_{𝒌,2}\widehat{a}_{𝒌,2})+{\displaystyle \frac{\alpha }{2}}|\zeta _0|^2\widehat{a}_{𝒌,0}\widehat{a}_{𝒌,0}+(\alpha \beta )|\zeta _2\zeta _2|\widehat{a}_{𝒌,2}\widehat{a}_{𝒌,2}`$
$`+(\alpha \beta +\gamma )|\zeta _2\zeta _2|\widehat{a}_{𝒌,2}^{}\widehat{a}_{𝒌,2}+\alpha |\zeta _0|(|\zeta _2|\widehat{a}_{𝒌,2}\widehat{a}_{𝒌,0}+|\zeta _2|\widehat{a}_{𝒌,0}\widehat{a}_{𝒌,2})`$
$`+|\zeta _0|(\alpha |\zeta _2|\gamma |\zeta _2|)\widehat{a}_{𝒌,2}^{}\widehat{a}_{𝒌,0}+|\zeta _0|(\alpha |\zeta _2|\gamma |\zeta _2|)\widehat{a}_{𝒌,0}^{}\widehat{a}_{𝒌,2}+\mathrm{H}.\mathrm{c}.\}`$
$`+{\displaystyle \frac{\beta }{4}}{\displaystyle \underset{𝒌\mathrm{𝟎}}{}}\{2\zeta _0^2\widehat{a}_{𝒌,1}\widehat{a}_{𝒌,1}+\sqrt{6}\zeta _0^{}(\zeta _2^{}\widehat{a}_{𝒌,1}\widehat{a}_{𝒌,1}+\zeta _2^{}\widehat{a}_{𝒌,1}\widehat{a}_{𝒌,1})+\sqrt{6}(\zeta _2^{}\zeta _0+\zeta _0^{}\zeta _2)\widehat{a}_{𝒌,1}^{}\widehat{a}_{𝒌,1}+\mathrm{H}.\mathrm{c}.\},`$
where $`\alpha c_0n`$, $`\beta 4c_1n`$, $`\gamma 2c_2n/5`$, and

$`A_{𝒌,\pm 2}`$ $`=`$ $`ϵ_𝒌+(\alpha +\beta )|\zeta _{\pm 2}|^2+\gamma |\zeta _2|^2,`$ (219)
$`A_{𝒌,\pm 1}`$ $`=`$ $`ϵ_𝒌+{\displaystyle \frac{\beta }{4}}(2|\zeta _{\pm 2}|^2+3|\zeta _0|^2),`$ (220)
$`A_{𝒌,0}`$ $`=`$ $`ϵ_𝒌+(\alpha +\gamma )|\zeta _0|^2.`$ (221)
It can be seen from the Hamiltonian (218) that there are two separate sets of coupled modes, that is, the $`m=\pm 1`$ modes and the $`m=0,\pm 2`$ modes.

#### 1 The $`m=\pm 1`$ coupled modes

The equations of motion governing the $`m=\pm 1`$ coupled modes are given by

$`i\mathrm{}{\displaystyle \frac{d}{dt}}\widehat{a}_{𝒌,1}`$ $`=`$ $`A_{𝒌,1}\widehat{a}_{𝒌,1}+D_2\widehat{a}_{𝒌,1}^{}+B\widehat{a}_{𝒌,1}^{}+C\widehat{a}_{𝒌,1}`$ (222)
$`i\mathrm{}{\displaystyle \frac{d}{dt}}\widehat{a}_{𝒌,1}^{}`$ $`=`$ $`A_{𝒌,1}\widehat{a}_{𝒌,1}D_2^{}\widehat{a}_{𝒌,1}B^{}\widehat{a}_{𝒌,1}C^{}\widehat{a}_{𝒌,1}^{}`$ (223)
$`i\mathrm{}{\displaystyle \frac{d}{dt}}\widehat{a}_{𝒌,1}`$ $`=`$ $`A_{𝒌,1}\widehat{a}_{𝒌,1}+D_2\widehat{a}_{𝒌,1}^{}+B\widehat{a}_{𝒌,1}^{}+C^{}\widehat{a}_{𝒌,1}`$ (224)
$`i\mathrm{}{\displaystyle \frac{d}{dt}}\widehat{a}_{𝒌,1}^{}`$ $`=`$ $`A_{𝒌,1}\widehat{a}_{𝒌,1}^{}D_2^{}\widehat{a}_{𝒌,1}B^{}\widehat{a}_{𝒌,1}C\widehat{a}_{𝒌,1}^{},`$ (225)
where $`B\beta \zeta _0^2/2`$, $`C\sqrt{6}\beta (\zeta _2\zeta _0^{}+\zeta _0\zeta _2^{})/4`$, $`D_{\pm 2}\sqrt{6}\beta \zeta _0\zeta _{\pm 2}/2`$. The eigenenergies of these modes are given by

$`(E_{\pm 1}^\mathrm{C})^2`$ $`=`$ $`ϵ_𝐤^2+\beta \left(1{\displaystyle \frac{\widehat{f}_z^2}{8}}\right)ϵ_𝒌+{\displaystyle \frac{\beta ^2}{32}}\widehat{f}_z^2\pm {\displaystyle \frac{\beta \widehat{f}_z}{2}}\left\{\left(ϵ_𝒌+{\displaystyle \frac{\beta }{4}}\right)\left[\left(1{\displaystyle \frac{3}{16}}\widehat{f}_z^2\right)ϵ_𝒌+{\displaystyle \frac{\beta }{64}}\widehat{f}_z^2\right]\right\}^{\frac{1}{2}}.`$ (226)
These excitation energies are always positive semidefinite and massive in the presence of magnetic field.

We note that Eq. (226) has one gapless (but not massless) mode in the presence of external magnetic field. Taking the limit $`|𝒌|\mathrm{𝟎}`$ of Eq. (226), we obtain

$`E_{\pm 1}^\mathrm{C}=\left({\displaystyle \frac{4}{\widehat{f}_z}}{\displaystyle \frac{\widehat{f}_z}{2}}\right)ϵ_𝒌+{\displaystyle \frac{\beta }{4}}\widehat{f}_z\mathrm{and}\left[\left(\widehat{f}_z{\displaystyle \frac{4}{\widehat{f}_z}}\right)^2+1\right]^{\frac{1}{2}}ϵ_𝒌.`$ (227)
However, both of these modes become massless in the absence of magnetic field. In fact, Eq. (226) then reduces to

$`E_{\pm 1}^\mathrm{C}`$ $`=`$ $`\sqrt{ϵ_𝒌(ϵ_𝒌+4c_1n)}.`$ (228)
#### 2 The $`m=\pm 2,0`$ coupled modes

The equations of motion governing the $`m=\pm 2,0`$ coupled modes are given by

$`i\mathrm{}{\displaystyle \frac{d}{dt}}\widehat{X}_+`$ $`=`$ $`(A_2(\alpha +\beta )\zeta _2^2)\widehat{X}_{}\gamma |\zeta _2\zeta _0|\widehat{Y}_{}+\gamma \zeta _0^2\widehat{Z}_{}/2,`$ (229)
$`i\mathrm{}{\displaystyle \frac{d}{dt}}\widehat{X}_{}`$ $`=`$ $`(A_2+(\alpha +\beta )\zeta _2^2)\widehat{X}_++(2\alpha \zeta _2\gamma \zeta _2)\zeta _0\widehat{Y}_++(\alpha \beta +\gamma /2)\zeta _0^2\widehat{Z}_+,`$ (230)
$`i\mathrm{}{\displaystyle \frac{d}{dt}}\widehat{Y}_+`$ $`=`$ $`\gamma |\zeta _2\zeta _0|\widehat{X}_{}+(A_0\alpha \zeta _0^2)\widehat{Y}_{}\gamma |\zeta _2\zeta _0|\widehat{Z}_{},`$ (231)
$`i\mathrm{}{\displaystyle \frac{d}{dt}}\widehat{Y}_{}`$ $`=`$ $`(2\alpha \zeta _2\gamma \zeta _2)\zeta _0\widehat{X}_++(A_0+\alpha \zeta _0^2)\widehat{Y}_++(2\alpha \zeta _2\gamma \zeta _2)\zeta _0\widehat{Z}_+,`$ (232)
$`i\mathrm{}{\displaystyle \frac{d}{dt}}\widehat{Z}_+`$ $`=`$ $`\gamma \zeta _0^2\widehat{X}_{}/2\gamma |\zeta _2\zeta _0|\widehat{Y}_{}+(A_2(\alpha +\beta )\zeta _2^2)\widehat{Z}_{},`$ (233)
$`i\mathrm{}{\displaystyle \frac{d}{dt}}\widehat{Z}_{}`$ $`=`$ $`(\alpha \beta +\gamma /2)\zeta _0^2\widehat{X}_++(2\alpha \zeta _2\gamma \zeta _2)\zeta _0\widehat{Y}_++(A_2+(\alpha +\beta ))\widehat{Z}_+,`$ (234)
where $`\widehat{X}_\pm \widehat{a}_{𝒌,2}\pm \widehat{a}_{𝒌,2}^{}`$, $`\widehat{Y}_\pm \widehat{a}_{𝒌,0}\pm \widehat{a}_{𝒌,0}^{}`$, and $`\widehat{Z}_\pm \widehat{a}_{𝒌,2}\pm \widehat{a}_{𝒌,2}^{}`$. Substituting the expressions for $`\widehat{X}_{},\widehat{Y}_{}`$ and $`\widehat{Z}_{}`$ into those for $`\widehat{X}_+,\widehat{Y}_+`$ and $`\widehat{Z}_+`$, we obtain the equations of motion for the latter set, which reduces to the cubic equation and therefore can be solved analytically. The result is given by

$`E_{\pm 2,0}^\mathrm{C}=ϵ_𝒌+\gamma ,\left\{ϵ_𝒌\left[ϵ_𝐤+\alpha +{\displaystyle \frac{4+\widehat{f}_z^2}{8}}\beta \pm \sqrt{\alpha ^2\left(1{\displaystyle \frac{3}{4}}\widehat{f}_z^2\right)\alpha \beta +{\displaystyle \frac{(4+\widehat{f}_z^2)^2}{64}}\beta ^2}\right]\right\}^{\frac{1}{2}},`$ (235)
where we recall that $`\alpha c_0n`$, $`\beta 4c_1n`$, $`\gamma 2c_2n/5`$. The second solutions in Eq. (235) are always positive semidefinite. The positivity of the first solution is guaranteed by the condition $`c_2>0`$. We note that the first solution is massive and independent of the applied magnetic field and that the second solutions remain massless in the presence of external magnetic field. The latter is a consequence of the fact that the mean-field solution is degenerate with respect to at least two continuous variables, as discussed in Sec. III C.

In the absence of external magnetic field, the results (235) reduce to

$`E_{\pm 2,0}^\mathrm{C}=ϵ_𝒌+\gamma ,\sqrt{ϵ_𝒌(ϵ_𝒌+2\alpha )},\sqrt{ϵ_𝒌(ϵ_𝒌+\beta )}.`$ (236)
It can be shown from the general analytic solutions that these results are valid up to the first order in magnetization $`\widehat{f}_z`$.

While we have been unable to complete the analysis of the cyclic phase except for the case of Eq. (210), we would like to point out that the excitation spectrum of the cyclic phase always includes the first solution in Eq. (235), even when the mean-field solution is not given by Eq. (210). This can be seen directly by writing down the equation of motion for $`\widehat{a}_{𝒌m}`$ and $`\widehat{a}_{𝒌m}^{}`$ and the corresponding eigenvalue equation. It can then be seen that $`ϵ_𝒌+\gamma `$ is a solution to this equation.

## VI Conclusions

In this paper we have studied quantum spin correlations and magnetic response of spin-2 Bose-Einstein condensates (BECs) in a mesoscopic regime, and low-lying excitation spectra of each phase of spin-2 BECs in the thermodynamic regime.

The ground states of spin-2 BECs have three distinct phases: ferromagnetic (FM), antiferromagnetic (AF) and cyclic (C) phases. The former two phases appear also in spin-1 BECs, while the last phase is unique to spin-2 BECs. The building block of the AF phase is spin-singlet pairs and that of the C phase is spin-singlet trios. These many-body features usually elude mean-field treatments that are based on the order parameter derived from the single-particle density matrix.

Bose symmetry restricts possible building blocks of spin-2 BECs. This can be summarized in terms of the annihilation operator $`\widehat{A}_f^{(n)}`$ of $`n`$-bosons having total spin $`f`$. The fundamental building block is not unique, but one minimal set is $`\widehat{A}_2^{(1)}`$, $`\widehat{A}_0^{(2)}`$, $`\widehat{A}_2^{(2)}`$, $`\widehat{A}_0^{(3)}`$, and $`\widehat{A}_3^{(3)}`$. Bose statistics does not allow units such as $`\widehat{A}_1^{(2)}`$ and $`\widehat{A}_1^{(3)}`$. The unit $`\widehat{A}_3^{(3)}`$ is required to represent a state with odd values of the total spin.

We have investigated quantum spin correlations and magnetic response in the mesoscopic regime. Under the assumption that the system is so tightly confined that the spatial degrees of freedom are frozen, we derived the exact many-body ground states which are expressed in terms of the minimal set of creation operators $`\widehat{A}_2^{(1)}`$, $`\widehat{A}_0^{(2)}`$, $`\widehat{A}_2^{(2)}`$, $`\widehat{A}_0^{(3)}`$, and $`\widehat{A}_3^{(3)}`$. These pairwise and trio-wise units help us understand the complicated response of the magnetization to the applied magnetic field, which stems from the fact that several values of the magnetization cannot be constructed from such units and are hence forbidden. In addition to the quantization of the magnetization to discrete values, several new features which elude mean-field treatments are found, such as a sudden jump from the minimum to the maximum magnetization, and robustness of the minimum-magnetization state against a small increase in the applied magnetic field until it starts to show a linear response. The average Zeeman level populations for the AF-phase ground states were calculated, showing that $`m=0,\pm 1`$ populations, which stay zero in MFT, vary sensitively to the applied magnetic field.

We have also examined low-lying excitation spectra using the Bogoliubov approximation. The excitation spectra of FM and AF phases are similar to those of the spin-1 case . In the FM phase, the spectrum consists of one massless mode (174) reflecting the global gauge invariance and four single-particle modes (175)-(178) whose energy gaps are generated by the Zeeman shifts as well as mean-field interactions. In the AF phase, the spectrum includes two massless modes (206) due to the global gauge invariance and the rotational symmetry about the spin quantization axis. The remaining three are also Bogoliubov modes, but they all become massive in the presence of magnetic field due to the Zeeman shifts. In the C phase, the spectrum have at least two massless modes (the second term in Eq. (235)) by the same reasons as in the AF phase. The spectrum includes one peculiar single-particle mode (the first term in Eq. (235)) whose energy gap depends solely on the spin-dependent interactions and is insensitive to the applied magnetic field. In addition, the spectrum has one gapless mode (the second term in Eq. (227)) whose mass depends only on magnetization $`\widehat{f}_z`$ and vanishes at zero magnetic field. The remaining mode (the first term in Eq. (227)) is a Bogoliubov mode which becomes massive in the presence of the magnetic field.

In the present paper we have studied only static properties of spin-2 BEC. With the very rich phenomena that we have found here, we may very well expect that much more remains to be revealed in their dynamics.

###### Acknowledgements.

This work was supported by a Grant-in-Aid for Scientific Research (Grant No. 11216204) by the Ministry of Education, Science, Sports, and Culture of Japan, and by the Toray Science Foundation. M.U. acknowledges the hospitality of the Aspen Center for Physics, where part of this work was carried out.

## A The order parameter of spin-2 BEC

The order parameter $`\zeta `$ of a spin-2 BEC has the same structure as that of the d-wave superconductor which was examined by Mermin . We here recapitulate as much of it as is relevant to the our theory. The spin part of the order parameter, $`\mathrm{\Psi }`$, of a spin-2 BEC is described as a function of the azimuthal angle $`\theta `$ and the radial angle $`\varphi `$, and may be expanded in terms of the spherical harmonics of rank 2 $`Y_2^m(\theta ,\varphi )`$ as

$`\mathrm{\Psi }={\displaystyle \underset{m=2}{\overset{2}{}}}\zeta _mY_2^m,`$ (A1)
where $`\zeta _m`$ obeys the normalization condition (46). The angle dependence of $`Y_2^m`$ may be expressed in terms of components of a three-dimensional unit vector: $`\widehat{𝐧}^T=(k_x,k_y,k_z)(\mathrm{sin}\theta \mathrm{cos}\varphi ,\mathrm{sin}\theta \mathrm{sin}\varphi ,\mathrm{cos}\theta )`$ as follows

$`Y_2^{\pm 2}=\sqrt{{\displaystyle \frac{15}{32\pi }}}(k_x\pm ik_y)^2,Y_2^{\pm 1}=\sqrt{{\displaystyle \frac{15}{8\pi }}}(k_x\pm ik_y)k_z,Y_2^0=\sqrt{{\displaystyle \frac{5}{16\pi }}}(2k_z^2k_x^2k_y^2).`$ (A2)
Substituting these into Eq. (A1), we obtain

$`\mathrm{\Psi }=\sqrt{{\displaystyle \frac{15}{8\pi }}}\widehat{𝐧}^T𝐌\widehat{𝐧},`$ (A3)
where

$`𝐌={\displaystyle \frac{1}{2}}\left(\begin{array}{ccc}\zeta _2+\zeta _2\sqrt{\frac{2}{3}}\zeta _0& i(\zeta _2\zeta _2)& \zeta _1+\zeta _1\\ i(\zeta _2\zeta _2)& \zeta _2\zeta _2\sqrt{\frac{2}{3}}\zeta _0& i(\zeta _1+\zeta _1)\\ \zeta _1+\zeta _1& i(\zeta _1+\zeta _1)& 2\sqrt{\frac{2}{3}}\zeta _0\end{array}\right).`$ (A7)
The order parameter is thus characterized by a $`3\times 3`$ traceless matrix Tr$`𝐌=0`$ with unit normalization

$`\mathrm{Tr}(𝐌^{}𝐌)=1.`$ (A8)
We may exploit the freedom of the gauge invariance to choose the overall phase so that the real part of Tr$`𝐌^2`$ vanishes.

$`\mathrm{ReTr}𝐌^2=0.`$ (A9)
Let the real and imaginary parts of $`𝐌`$ be $`𝐗`$ and $`𝐘`$, respectively. It follows from Eqs. (A8) and (A9) that

$`\mathrm{Tr}𝐗^2=\mathrm{Tr}𝐘^2={\displaystyle \frac{1}{2}}.`$ (A10)
Because $`𝐁`$ is traceless, so can be $`𝐗`$ and $`\mathrm{Y}`$ traceless. $`𝐗`$ and $`\mathrm{Y}`$ do not commute, so they cannot be diagonalized simultaneously. We follow Mermin to take a representation in which $`𝐗`$ is diagonal. Then the diagonal elements of $`𝐗`$ become

$`x_n=\sqrt{{\displaystyle \frac{1}{3}}}\mathrm{sin}\left(\theta +{\displaystyle \frac{2\pi }{3}}n\right).`$ (A11)
The matrix elements of $`𝐘`$ are given by

$`Y_{nn}=\sqrt{{\displaystyle \frac{1}{3}}}\mathrm{sin}\left(\varphi +{\displaystyle \frac{2\pi }{3}}n\right)\mathrm{sin}\chi ,Y_{23}=Y_{32}={\displaystyle \frac{1}{2}}\mathrm{sin}\delta \mathrm{cos}\psi \mathrm{cos}\chi ,`$ (A12)
$`Y_{31}=Y_{13}={\displaystyle \frac{1}{2}}\mathrm{sin}\delta \mathrm{sin}\psi \mathrm{cos}\chi ,Y_{12}=Y_{21}={\displaystyle \frac{1}{2}}\mathrm{cos}\delta \mathrm{cos}\chi .`$ (A13)
Substituting Eqs. (A10) and (A13) into $`𝐌=𝐗+i𝐘`$ and comparing it with Eq. (A7), we obtain the following characterization of the order parameter.

$`\zeta _{\pm 2}={\displaystyle \frac{1}{2}}(\mathrm{cos}\theta \pm \mathrm{cos}\delta \mathrm{cos}\chi +i\mathrm{cos}\varphi \mathrm{sin}\chi )`$ (A14)
$`\zeta _{\pm 1}={\displaystyle \frac{1}{2}}\mathrm{cos}\chi \mathrm{sin}\delta e^{\pm i\psi }`$ (A15)
$`\zeta _0={\displaystyle \frac{1}{\sqrt{2}}}(\mathrm{sin}\theta +i\mathrm{sin}\chi \mathrm{sin}\varphi ).`$ (A16)
## B Calculation of Zeeman-level populations

### 1 Derivation of Eq. (147)

Let us first write $`\widehat{a}_m^{}\widehat{a}_m`$ as

$$\widehat{a}_m^{}\widehat{a}_m=\frac{\varphi _m|(\widehat{𝒮}_{})^{N_\mathrm{S}}(\widehat{𝒮}_+)^{N_\mathrm{S}}|\varphi _m}{\varphi |(\widehat{𝒮}_{})^{N_\mathrm{S}}(\widehat{𝒮}_+)^{N_\mathrm{S}}|\varphi }1$$
(B1)
with $`|\varphi _m\widehat{a}_m^{}|\varphi `$. What we need is thus a general formula for calculating $`\psi |(\widehat{𝒮}_{})^l(\widehat{𝒮}_+)^l|\psi `$. Let us first consider the decomposition of $`|\psi `$ into a sum of eigenstates for $`\widehat{𝒮}^2`$, such that

$$|\psi =\underset{k=0}{\overset{n/2}{}}(\widehat{𝒮}_+)^k|\psi _k,$$
(B2)
where $`|\psi _k`$ is an unnormalized simultaneous eigenstate of $`\{\widehat{𝒮}^2,\widehat{𝒮}_z\}`$ with eigenvalues $`S(S1)`$ and $`S=[2(n2k)+5]/4`$, respectively. It follows from Eqs. (97) and (98) that $`N_s=0`$ and $`N_0=n2k`$. Here $`n`$ is the number of bosons in $`|\psi `$, and $`x`$ denotes the largest integer that is not larger than $`x`$. Let us define $`\omega _j`$ and $`\mu _j`$ such that

$$\omega _j\frac{\psi |(\widehat{𝒮}_+)^j(\widehat{𝒮}_{})^j|\psi }{\psi |\psi }$$
(B3)
and

$$\mu _j\frac{\psi _j|(\widehat{𝒮}_{})^j(\widehat{𝒮}_+)^j|\psi _j}{\psi |\psi }=(n2j;j)\frac{\psi _j|\psi _j}{\psi |\psi },$$
(B4)
where we have used Eq. (104) and defined the coefficient $`(a;b)`$ as

$$(a;b)\frac{b!(b+a+3/2)!}{(a+3/2)!}.$$
(B5)
By definition, $`\omega _0=1`$ and $`_j\mu _j=1`$. Substituting Eq. (B3) into Eq. (B2) yields

$$\omega _j=(n2j;j)\mu _j+\underset{k=j+1}{\overset{n/2}{}}\mu _k\frac{(n2k;k)}{(n2k;kj)},$$
(B6)
or equivalently,

$$\mu _j=\frac{\omega _j}{(n2j;j)}\underset{k=j+1}{\overset{n/2}{}}\mu _k\frac{(n2k;k)}{(n2j;j)(n2k;kj)}.$$
(B7)
Using this relation recursively, we can calculate $`\mu _j`$ as a function of $`\{\omega _k\}(k=j,j+1,\mathrm{},n/2)`$. On the other hand, multiplying $`(\widehat{𝒮}_+)^l`$ on both sides of Eq. (B2) and taking norms, we obtain

$$\frac{\psi |(\widehat{𝒮}_{})^l(\widehat{𝒮}_+)^l|\psi }{\psi |\psi }=\underset{k=0}{\overset{n/2}{}}\frac{(n2k;l+k)}{(n2k;k)}\mu _k.$$
(B8)
Using this formula, we can evaluate Eq. (B1). If we apply $`|\psi =|\varphi `$, $`n=s`$, and $`l=N_\mathrm{S}`$ to the relations (B3), (B7), and (B8), and noting that $`\widehat{𝒮}_{}|\varphi =0`$, we have $`\omega _0=\mu _0=1`$, $`\omega _j=\mu _j=0(j>0)`$, and

$$\frac{\varphi |(\widehat{𝒮}_{})^{N_\mathrm{S}}(\widehat{𝒮}_+)^{N_\mathrm{S}}|\varphi }{\varphi |\varphi }=(s;N_\mathrm{S}).$$
(B9)
For $`|\varphi _m`$, we have $`\widehat{𝒮}_{}|\varphi _m=[\widehat{𝒮}_{},\widehat{a}_m^{}]|\varphi =(1)^m\widehat{a}_m|\varphi `$ and $`(\widehat{𝒮}_{})^2|\varphi _m=0`$. Using these when we apply $`|\psi =|\varphi _m`$, $`n=s+1`$, and $`l=N_\mathrm{S}`$ to the relations (B3), (B7), and (B8), we have $`\omega _1=\varphi |\widehat{a}_m^{}\widehat{a}_m|\varphi /\varphi _m|\varphi _m`$, $`\omega _j=\mu _j=0(j>1)`$, $`\mu _1=\omega _1/(s1;1)`$, $`\mu _0=1\mu _1`$, and

$$\frac{\varphi _m|(\widehat{𝒮}_{})^{N_\mathrm{S}}(\widehat{𝒮}_+)^{N_\mathrm{S}}|\varphi _m}{\varphi _m|\varphi _m}=(s+1;N_\mathrm{S})\mu _0+\frac{(s1;N_\mathrm{S}+1)}{(s1;1)}\mu _1.$$
(B10)
Combining Eqs. (B1), (B9), and (B10), we obtain

$$\widehat{a}_m\widehat{a}_m^{}=\frac{(s+1;N_\mathrm{S})}{(s;N_\mathrm{S})}\left(\widehat{a}_m\widehat{a}_m^{}_0\frac{\widehat{a}_m^{}\widehat{a}_m_0}{(s1;1)}\right)+\frac{(s1;N_\mathrm{S}+1)}{(s;N_\mathrm{S})(s1;1)^2}\widehat{a}_m^{}\widehat{a}_m_0,$$
(B11)
where $`\mathrm{}_0\varphi |\mathrm{}|\varphi /\varphi |\varphi `$. Substituting Eq.(B5), we obtain Eq. (147).

### 2 Exact forms of Zeeman populations

Here we give the results of Zeeman populations $`\widehat{a}_m^{}\widehat{a}_m`$ for the states $`(\widehat{A}_0^{(2)})^{N_\mathrm{S}}(\widehat{a}_2^{})^{n_{12}}(\widehat{A}_2^{(2)})^{n_{22}}(\widehat{A}_3^{(3)})^{n_{33}}|\mathrm{vac}`$ with $`n_{22}=0,1`$ and $`n_{33}=0,1`$.

i) $`n_{22}=0`$ and $`n_{33}=0`$. We take $`|\varphi =(\widehat{a}_2^{})^{n_{12}}|\mathrm{vac}`$ and $`s=n_{12}`$. Then, $`\widehat{a}_2^{}\widehat{a}_2_0=n_{12}`$ and $`\widehat{a}_m^{}\widehat{a}_m_0=0`$ for $`m<2`$. Putting these into the formula (147), we obtain

$$\widehat{a}_2^{}\widehat{a}_2=n_{12}+N_\mathrm{S}\frac{3N_\mathrm{S}}{2n_{12}+5},$$
(B12)
$$\widehat{a}_1^{}\widehat{a}_1=\widehat{a}_0^{}\widehat{a}_0=\widehat{a}_1^{}\widehat{a}_1=\frac{2N_\mathrm{S}}{2n_{12}+5},$$
(B13)
$$\widehat{a}_2^{}\widehat{a}_2=N_\mathrm{S}\frac{3N_\mathrm{S}}{2n_{12}+5}.$$
(B14)
ii) $`n_{22}=1`$ and $`n_{33}=0`$. We take $`|\varphi =(\widehat{a}_2^{})^{n_{12}}[\sqrt{3}(\widehat{a}_1)^22\sqrt{2}\widehat{a}_2\widehat{a}_0]^{}|\mathrm{vac}`$ and $`s=n_{12}+2`$. Then,

$$\widehat{a}_2^{}\widehat{a}_2_0=n_{12}+1\frac{3}{4n_{12}+7},$$
(B15)
$$\widehat{a}_1^{}\widehat{a}_1_0=\frac{6}{4n_{12}+7},$$
(B16)
$$\widehat{a}_0^{}\widehat{a}_0_0=1\frac{3}{4n_{12}+7},$$
(B17)
and $`\widehat{a}_m^{}\widehat{a}_m_0=0`$ for $`m<0`$. Putting these into the formula (147), we obtain

$$\widehat{a}_2^{}\widehat{a}_2=n_{12}+N_\mathrm{S}\frac{5N_\mathrm{S}}{2n_{12}+9}\left(\frac{2N_\mathrm{S}}{2n_{12}+9}+1\right)\frac{3}{4n_{12}+7}+1,$$
(B18)
$$\widehat{a}_1^{}\widehat{a}_1=\frac{2N_\mathrm{S}}{2n_{12}+9}+\left(\frac{2N_\mathrm{S}}{2n_{12}+9}+1\right)\frac{6}{4n_{12}+7},$$
(B19)
$$\widehat{a}_0^{}\widehat{a}_0=\frac{6N_\mathrm{S}}{2n_{12}+9}\left(\frac{4N_\mathrm{S}}{2n_{12}+9}+1\right)\frac{3}{4n_{12}+7}+1,$$
(B20)
$$\widehat{a}_1^{}\widehat{a}_1=\frac{2N_\mathrm{S}}{2n_{12}+9}\left(1+\frac{6}{4n_{12}+7}\right),$$
(B21)
$$\widehat{a}_2^{}\widehat{a}_2=N_\mathrm{S}\frac{N_\mathrm{S}}{2n_{12}+9}\left(5+\frac{6}{4n_{12}+7}\right).$$
(B22)
iii) $`n_{22}=0`$ and $`n_{33}=1`$. We take $`|\varphi =(\widehat{a}_2^{})^{n_{12}}[(\widehat{a}_1)^3\sqrt{6}\widehat{a}_2\widehat{a}_1\widehat{a}_0+2(\widehat{a}_2)^2\widehat{a}_1]^{}|\mathrm{vac}`$ and $`s=n_{12}+3`$. Then,

$$\widehat{a}_2^{}\widehat{a}_2_0=n_{12}+2\frac{3(n_{12}+3)}{(n_{12}+2)(2n_{12}+5)},$$
(B23)
$$\widehat{a}_1^{}\widehat{a}_1_0=\frac{3(n_{12}+4)}{(n_{12}+2)(2n_{12}+5)},$$
(B24)
$$\widehat{a}_0^{}\widehat{a}_0_0=\frac{3(n_{12}+1)}{(n_{12}+2)(2n_{12}+5)},$$
(B25)
$$\widehat{a}_1^{}\widehat{a}_1_0=1\frac{3}{2n_{12}+5},$$
(B26)
and $`\widehat{a}_2^{}\widehat{a}_2_0=0`$. Putting these into the formula (147), we obtain

$$\widehat{a}_2^{}\widehat{a}_2=n_{12}+N_\mathrm{S}\frac{5N_\mathrm{S}}{2n_{12}+11}\left(\frac{2N_\mathrm{S}}{2n_{12}+11}+1\right)\frac{3(n_{12}+3)}{(n_{12}+2)(2n_{12}+5)}+2,$$
(B27)
$$\widehat{a}_1^{}\widehat{a}_1=\frac{4N_\mathrm{S}}{2n_{12}+11}+\left(\frac{12N_\mathrm{S}}{2n_{12}+11}+3(n_{12}+4)\right)\frac{1}{(n_{12}+2)(2n_{12}+5)},$$
(B28)
$$\widehat{a}_0^{}\widehat{a}_0=\frac{2N_\mathrm{S}}{2n_{12}+11}+\left(\frac{4N_\mathrm{S}}{2n_{12}+11}+1\right)\frac{3(n_{12}+1)}{(n_{12}+2)(2n_{12}+5)},$$
(B29)
$$\widehat{a}_1^{}\widehat{a}_1=\frac{4N_\mathrm{S}}{2n_{12}+11}+\left(\frac{12N_\mathrm{S}}{2n_{12}+11}3(n_{12}+2)\right)\frac{1}{(n_{12}+2)(2n_{12}+5)}+1,$$
(B30)
$$\widehat{a}_2^{}\widehat{a}_2=N_\mathrm{S}\frac{N_\mathrm{S}}{2n_{12}+11}\left(5+\frac{6(n_{12}+3)}{(n_{12}+2)(2n_{12}+5)}\right).$$
(B31)
iv) $`n_{22}=1`$ and $`n_{33}=1`$. We take $`|\varphi =(\widehat{a}_2^{})^{n_{12}}[\sqrt{3}(\widehat{a}_1)^55\sqrt{2}\widehat{a}_2(\widehat{a}_1)^3\widehat{a}_0+4\sqrt{3}(\widehat{a}_2)^2\widehat{a}_1(\widehat{a}_0)^2+2\sqrt{3}(\widehat{a}_2)^2(\widehat{a}_1)^2\widehat{a}_14\sqrt{2}(\widehat{a}_2)^3\widehat{a}_0\widehat{a}_1]^{}|\mathrm{vac}`$ and $`s=n_{12}+5`$. Then,

$$\widehat{a}_2^{}\widehat{a}_2_0=n_{12}+3\frac{30(n_{12}+4)^2}{(n_{12}+3)(2n_{12}+7)(4n+13)},$$
(B32)
$$\widehat{a}_1^{}\widehat{a}_1_0=\frac{9(4n_{12}^2+37n_{12}+83)}{(n_{12}+3)(2n_{12}+7)(4n+13)},$$
(B33)
$$\widehat{a}_0^{}\widehat{a}_0_0=1+\frac{18(n_{12}^2+3n_{12}3)}{(n_{12}+3)(2n_{12}+7)(4n+13)},$$
(B34)
$$\widehat{a}_1^{}\widehat{a}_1_0=1\frac{3(8n_{12}^2+49n_{12}+71)}{(n_{12}+3)(2n_{12}+7)(4n+13)},$$
(B35)
and $`\widehat{a}_2^{}\widehat{a}_2_0=0`$. Putting these into the formula (147), we obtain

$$\widehat{a}_2^{}\widehat{a}_2=n_{12}+N_\mathrm{S}\frac{7N_\mathrm{S}}{2n_{12}+15}\left(\frac{2N_\mathrm{S}}{2n_{12}+15}+1\right)\frac{30(n_{12}+4)^2}{(n_{12}+3)(2n_{12}+7)(4n+13)}+3,$$
(B36)
$$\widehat{a}_1^{}\widehat{a}_1=\frac{4N_\mathrm{S}}{2n_{12}+15}+\left(4N_\mathrm{S}(n_{12}+8)\frac{124N_\mathrm{S}}{2n_{12}+15}+3(4n_{12}^2+37n_{12}+83)\right)\frac{3}{(n_{12}+3)(2n_{12}+7)(4n+13)},$$
(B37)
$$\widehat{a}_0^{}\widehat{a}_0=\frac{6N_\mathrm{S}}{2n_{12}+15}+\left(\frac{4N_\mathrm{S}}{2n_{12}+15}+1\right)\frac{18(n_{12}^2+3n_{12}3)}{(n_{12}+3)(2n_{12}+7)(4n+13)}+1,$$
(B38)
$$\widehat{a}_1^{}\widehat{a}_1=\frac{4N_\mathrm{S}}{2n_{12}+15}+\left(4N_\mathrm{S}(n_{12}+8)\frac{124N_\mathrm{S}}{2n_{12}+15}(8n_{12}^2+49n_{12}+71)\right)\frac{3}{(n_{12}+3)(2n_{12}+7)(4n+13)}+1,$$
(B39)
$$\widehat{a}_2^{}\widehat{a}_2=N_\mathrm{S}\frac{N_\mathrm{S}}{2n_{12}+15}\left(7+\frac{60(n_{12}+4)^2}{(n_{12}+3)(2n_{12}+7)(4n+13)}\right).$$
(B40)