# Phase transitions and spin excitations of spin-1 bosons in optical lattice

Min-Jie Zhu  Shanghai Branch, National Laboratory for Physical Sciences at
Microscale and Department of Modern Physics, University of Science and
Technology of China, Hefei, Anhui 230026, China  CAS Center for Excellence and
Synergetic Innovation Center in Quantum Information and Quantum Physics,
University of Science and Technology of China, Shanghai 201315, China  CAS-
Alibaba Quantum Computing Laboratory, Shanghai 201315, China  Bo Zhao
Shanghai Branch, National Laboratory for Physical Sciences at Microscale and
Department of Modern Physics, University of Science and Technology of China,
Hefei, Anhui 230026, China  CAS Center for Excellence and Synergetic
Innovation Center in Quantum Information and Quantum Physics, University of
Science and Technology of China, Shanghai 201315, China  CAS-Alibaba Quantum
Computing Laboratory, Shanghai 201315, China

###### Abstract

We investigate ground state properties of spin-1 bosonic system trapped in
optical lattice with extended standard basis operator (SBO) method. For both
ferromagnetic ($U_{2}<0$) and antiferromagnetic ($U_{2}>0$) systems, we
analytically figure out the symmetry properties in Mott-insulator and
superfluid phases, which would provide a deeper insight into the MI-SF phase
transition process. Then by applying self-consistent approach to the method,
we include the effect of quantum and thermal fluctuations and derive the MI-SF
transition phase diagram, which is in quantitative agreement with recent
Monte-Carlo simulation at zero temperature, and at finite temperature, we find
the underestimation of finite-temperature-effect in the mean-field
approximation method. If we further consider the spin excitations in the
insulating states of spin-1 system in external field, distinct spin phases are
expected. Therefore, in the Mott lobes with $n=1$ and $n=2$ atoms per site, we
give analytical and numerical boundaries of the singlet, nematic, partially
magnetic and ferromagnetic phases in the magnetic phase diagrams.

##  I Introduction

Ultracold atoms trapped in optical lattices provide clean and highly
controllable platform for the study of strong-correlated systems Bloch et al.
(2008). Besides the well-known Mott insulator and superfluid phases Jaksch et
al. (1998); Greiner et al. (2002), recent experiments are approaching the
regime where magnetic properties could be revealed Jördens et al. (2010);
Medley et al. (2011); Jiang et al. (2016), such as the expected Neél
antiferromagnet in spin-1/2 fermions. The interplay of magnetism and
superfluidity in such systems can be realized with spinor gases trapped in
optical lattice. The presence of extra spin degrees of freedom induces
exceedingly rich physics in spinor system and allows the studies of multi-band
condensed matter Hamiltonians and symmetry and topology in quantum-ordered
materials Stamper-Kurn and Ueda (2013).

Most work in spinor systems, both theoretically and experimentally Hall et al.
(1998); Matthews et al. (1998); Widera et al. (2005, 2006); Ho (1998); Ohmi
and Machida (1998); Katsura and Tasaki (2013); Snoek and Zhou (2004); Demler
and Zhou (2002); Imambekov et al. (2003); Tsuchiya et al. (2004); Kimura et
al. (2005); Pai et al. (2008); Zhou and Semenoff (2006); Bernier et al.
(2007); Song and Zhou (2008); Jiang et al. (2016); Natu et al. (2015); So and
Ueda (2016), was focused on atoms with effective spin $F=1$, which is also
exclusively considered in this paper. In theory, the spin-1 atoms trapped in
optical lattice can be described by extended Bose-Hubbard model with spin-
dependent terms Imambekov et al. (2003); Demler and Zhou (2002); Stamper-Kurn
and Ueda (2013), in which both the MI phases and MI-SF phase transitions have
been intensively studied. For ferromagnetic interacting system ($U_{2}<0$),
the whole phase diagram, including Mott insulator and superfluid phases, is
predicted to be ferromagnetic Katsura and Tasaki (2013); de Forges de Parny et
al. (2013), and all the Mott lobes shrink, eventually disappear, as the
increase of spin-dependent interaction strength $|U_{2}|$ Batrouni et al.
(2009). While for antiferromagnetic interacting system ($U_{2}>0$), the
superfluid phase in the MI-SF transition is expected to be polarized with zero
global magnetization in the whole phase diagram Demler and Zhou (2002);
Katsura and Tasaki (2013); Snoek and Zhou (2004); Tsuchiya et al. (2004); de
Forges de Parny et al. (2013), and the Mott lobes with even filling numbers
grow at the expense of odd ones Tsuchiya et al. (2004); de Forges de Parny et
al. (2013). Then for small positive $U_{2}$, mean-field studies have pointed
out that the MI-SF transition from even filling Mott lobes is first-order
phase transition process (second-order for all other cases) Snoek and Zhou
(2004); Kimura et al. (2005), which has also been confirmed in experiment
Jiang et al. (2016) and one-dimensional quantum Monte-Carlo simulation
Batrouni et al. (2009). Moreover, spin phases including singlet, nematic and
dimerized phases are predicted in the Mott lobes for antiferromagnetic case
Demler and Zhou (2002); Imambekov et al. (2003, 2004). If the filling number
is odd, the Mott lobes are expected to form nematic state, which can be
changed to dimerized state in one-dimension system Imambekov et al. (2003);
Rizzi et al. (2005); Apaja and Syljuåsen (2006). If the filling number is
even, the Mott lobes are expected to have spin singlets in each site with a
first-order transition to nematic phase when $zU_{2}/U_{0}<0.1$
($z=2\text{dim}$) Imambekov et al. (2003).

If external magnetic field is applied to the system, numerous new physics can
be found. For example, the transitions for bosons in Mott lobes with different
filling numbers may occur into different superfluid phases with spins aligned
along or opposite to the field direction Svidzinsky and Chui (2003). In the
insulating states, 3D mean-field analysis Chung and Yip (2009) predicted
nematic to ferromagnetic (or partially magnetic) transitions, then in Ref.
Rodríguez et al. (2011), the whole magnetic phase diagram is shown in the
system at unit filling with quadratic Zeeman effect. In spite of these
findings, a general analytical theory of the field-induced quantum phases in
spin-1 bosonic system is still missing.

In this paper, we would like to propose extended standard basis operator (SBO)
method to solve the spin-1 bosonic system trapped in optical lattice. This
method, which is applied to the motion equation for Green’s functions Haley
and Erdös (1972); Sajna et al. (2015), has the capacity to analytically solve
most of the phase transition problems in spin-1 system on one hand, on the
other hand, it is possible to handle quantum fluctuations in the system by
taking self-consistent approach, showing that the method goes significantly
beyond mean-field level. Therefore, we apply it to study the symmetry
properties (polar or ferromagnetic) of superfluid and Mott-insulator phases in
both ferromagnetic ($U_{2}<0$) and antiferromagnetic ($U_{2}>0$) systems. The
MI-SF transition phase diagrams are given afterwards, and we find the self-
consistent SBO phase diagrams at zero temperature are in quantitative
agreement with quantum Monte-Carlo simulation results. Then in the Mott lobes
of the system in external magnetic field, of which the effective Hamiltonian
is derived with perturbation theory, we consider the spin excitations and
calculate the excitation spectrums in the system with SBO method, this enables
us to determine the boundaries of different phases in the magnetic phase
diagrams.

The rest of the paper is organized as follows. In Sec. II, we give a brief
introduction of the spin-1 bosonic system trapped in an optical lattice. In
Sec. III, we consider the MI-SF phase transition process with SBO method and
analyze the symmetry properties in the insulating and superfluid phases, and
we also give the self-consistent SBO phase diagrams in this section. In Sec.
IV, we analyze the magnetic phase diagram in the first ($n=1$) Mott lobe with
linear and quadratic Zeeman effect. In Sec. V, we analyze the magnetic phase
diagram in the second ($n=2$) Mott lobe with linear Zeeman effect. Finally, in
Sec. VI, we summarize our results and give some conclusions. Some of the
technical calculation details are presented in the appendix section.

##  II spin-1 model

Spin-1 bosons trapped in a spin-independent optical lattice can be described
with Bose-Hubbard Hamiltonian supplemented by spinor interaction in each
lattice site Imambekov et al. (2003); Demler and Zhou (2002); Stamper-Kurn and
Ueda (2013). The Hamiltonian is given by

$\displaystyle\hat{H}$ $\displaystyle=$ $\displaystyle-t\sum_{\langle
ij\rangle,\sigma}(\hat{a}^{\dagger}_{i\sigma}\hat{a}_{j\sigma}+\hat{a}^{\dagger}_{j\sigma}\hat{a}_{i\sigma})+\frac{U_{0}}{2}\sum_{i}\hat{n}_{i}(\hat{n}_{i}-1)$
(1)
$\displaystyle+\frac{U_{2}}{2}\sum_{i}(\vec{S}_{i}^{2}-2\hat{n}_{i})-\mu\sum_{i}\hat{n}_{i},$

where $\hat{a}_{i\sigma}$ is the annihilation operator of one boson in site
$i$ with spin-component $\sigma\in\\{-1,0,1\\}$,
$\hat{n}_{i}=\sum_{\sigma}\hat{a}^{\dagger}_{i\sigma}\hat{a}_{i\sigma}$, and
$\mu$ is the chemical potential that controls the total number of bosons in
the system. The total spin in site $i$ is defined as

$\displaystyle\vec{S}_{i}$ $\displaystyle=$
$\displaystyle\sum_{\sigma\sigma^{\prime}}\hat{a}^{\dagger}_{i\sigma}\vec{F}_{\sigma\sigma^{\prime}}\hat{a}_{i\sigma^{\prime}},$
(2)

where $\vec{F}_{\sigma\sigma^{\prime}}$ are the spin operators for spin-1
bosons.

The first term in Eq. (1) is the kinetic energy term that allows particles to
hop between neighbour sites, the second term describes spin-independent
interaction while the third term describes spin-exchange interaction which is
due to the difference in scattering lengths for $S=0$ and $S=2$ channels Ho
(1998). The spin-independent interaction is always repulsive ($U_{0}>0$), and
$U_{2}<0$ case (e.g. 87Rb) is referred to as ”ferromagnetic” case whereas the
$U_{2}>0$ case (e.g. 23Na) is referred to as ”antiferromagnetic”.

The on-site Hamiltonian in Eq. (1) (neglecting the hopping term) can be
described with eigenstate $|S_{i},m_{i};n_{i}\rangle$, which represents total
spin ($S_{i}$), spin-projection on z-axis ($m_{i}$) and total number of bosons
($n_{i}$) in site $i$ Law et al. (1998); Tsuchiya et al. (2004); Svidzinsky
and Chui (2003). The corresponding eigenenergy is given by

$\displaystyle E^{(0)}(S_{i},n_{i})$ $\displaystyle=$ $\displaystyle-\mu
n_{i}+\frac{U_{0}}{2}n_{i}(n_{i}-1)$ (3)
$\displaystyle+\frac{U_{2}}{2}\left[S_{i}(S_{i}+1)-2n_{i}\right].$

The states are degenerate with respect to different $m_{i}$. For each lattice
site, the eigenstates are constrained by $0\leq S_{i}\leq n_{i}$, $-S_{i}\leq
m_{i}\leq S_{i}$ and $S_{i}+n_{i}$ is an even number (due to the symmetry of
spin-wave function Wu (1996)).

In the MI phase, the ground state is determined by minimizing the on-site
energy $E^{(0)}(S_{i},n_{i})$ (Eq. (3)). For antiferromagnetic case
($U_{2}>0$), the ground state in MI phase can be $|0,0;n\rangle$ (for even
$n$) or $|1,m_{S};n\rangle$ (for odd $n$). While for ferromagnetic case
($U_{2}<0$), the ground state is $|n,m_{S};n\rangle$ for both even and odd
filling number $n$.

##  III MI-SF phase transition

For finite hopping energy $t$, the Mott insulator-superfluid phase transition
is induced by change in fluctuation of particle number in individual lattice
site. The transition process of spin-1 boson system can be calculated with
standard basis operator (SBO) method Haley and Erdös (1972); Sheshadri et al.
(1993); Sajna et al. (2015).

We start by expanding Hamiltonian in Eq. (1) on the basis of eigenstates of
on-site Hamiltonian (denoted as $|i\mu\rangle$ in expression)

$\displaystyle\hat{H}$ $\displaystyle=$
$\displaystyle\sum_{i\mu}E_{\mu}^{i}\hat{L}_{\mu\mu}^{i}-t\sum_{\langle
ij\rangle}\sum_{\mu\mu^{\prime}\nu\nu^{\prime}}T^{ij}_{\mu\mu^{\prime}\nu\nu^{\prime}}\hat{L}^{i}_{\mu\mu^{\prime}}\hat{L}^{j}_{\nu\nu^{\prime}},$
(4)

where $\hat{L}^{i}_{\mu\mu^{\prime}}=|i\mu\rangle\langle i\mu^{\prime}|$ and
$T^{ij}_{\mu\mu^{\prime}\nu\nu^{\prime}}=\sum_{\sigma}(c^{i\sigma}_{\mu\mu^{\prime}}d^{j\sigma}_{\nu\nu^{\prime}}+d^{i\sigma}_{\mu\mu^{\prime}}c^{j\sigma}_{\nu\nu^{\prime}})$
with the definitions $c^{i\sigma}_{\mu\mu^{\prime}}=\langle
i\mu|\hat{a}_{i\sigma}|i\mu^{\prime}\rangle$ and
$d^{i\sigma}_{\mu\mu^{\prime}}=\langle
i\mu|\hat{a}^{\dagger}_{i\sigma}|i\mu^{\prime}\rangle$. $E^{i}_{\mu}$ in Eq.
(4) is the engenenergy of the eigenstate $|i\mu\rangle$ in site $i$.

Then the retarded Green’s function can also be expanded on the basis of
$|i\mu\rangle$

$\displaystyle G^{i\alpha j\beta}(t-t^{\prime})$ $\displaystyle=$
$\displaystyle-i\Theta(t-t^{\prime})\langle[\hat{a}_{i\alpha}(t),\hat{a}^{\dagger}_{j\beta}(t^{\prime})]\rangle$
(5) $\displaystyle=$
$\displaystyle\sum_{\mu\mu^{\prime}\nu\nu^{\prime}}c^{i\alpha}_{\mu\mu^{\prime}}d^{j\beta}_{\nu\nu^{\prime}}G^{ij}_{\mu\mu^{\prime},\nu\nu^{\prime}}(t-t^{\prime}),$

where we have defined
$G^{ij}_{\mu\mu^{\prime},\nu\nu^{\prime}}(t-t^{\prime})=-i\Theta(t-t^{\prime})\langle[\hat{L}^{i}_{\mu\mu^{\prime}}(t),\hat{L}^{j}_{\nu\nu^{\prime}}(t^{\prime})]\rangle$.

For the homogeneous lattice with translational symmetry, we can eliminate
site-number and give the motion equation for Green’s function (see appendix.A)

$\displaystyle(\omega+E_{\mu}-E_{\mu^{\prime}})G_{\mu\mu^{\prime},\nu\nu^{\prime}}(\mathbf{k},\omega)-D_{\mu\mu^{\prime}}\delta_{\mu\nu^{\prime}}\delta_{\mu^{\prime}\nu}$
(6) $\displaystyle=$ $\displaystyle
D_{\mu\mu^{\prime}}\epsilon(\mathbf{k})\sum_{\xi\xi^{\prime}}T_{\mu^{\prime}\mu\xi\xi^{\prime}}G_{\xi\xi^{\prime},\nu\nu^{\prime}}(\mathbf{k},\omega),$

where $G_{\mu\mu^{\prime},\nu\nu^{\prime}}(\mathbf{k},\omega)$ is derived from
$G^{ij}_{\mu\mu^{\prime},\nu\nu^{\prime}}(t-t^{\prime})$ via Fourier
transformation, $D_{\mu\mu^{\prime}}=D_{\mu}-D_{\mu^{\prime}}$ with definition
$D_{\mu}=\langle\hat{L}_{\mu\mu}\rangle$ and
$\epsilon(\mathbf{k})=-2t\sum_{s=1}^{\text{dim}}\cos{k_{s}l}$ ($l$ is the
length between neighbour sites in the lattice).

Then the Green’s function which describes annihilation of spin-$\alpha$ and
construction of spin-$\beta$ bosons can be derived from the Fourier
transformation form of Eq. (5)

$\displaystyle G^{\alpha\beta}(\mathbf{k},\omega)$ $\displaystyle=$
$\displaystyle\sum_{\mu\mu^{\prime}\nu\nu^{\prime}}c_{\mu\mu^{\prime}}^{\alpha}d_{\nu\nu^{\prime}}^{\beta}G_{\mu\mu^{\prime},\nu\nu^{\prime}}(\mathbf{k},\omega).$
(7)

In the calculation of $G^{\alpha\beta}(\mathbf{k},\omega)$, we can treat it as
matrix element of a $3\times 3$ matrix $\mathbf{G}(\mathbf{k},\omega)$ with
the definition $G^{\alpha\beta}=[\mathbf{G}]_{\alpha\beta}$. Then matrix
$\mathbf{G}(\mathbf{k},\omega)$ can be solved (in appendix.A) as

$\displaystyle\mathbf{G}(\mathbf{k},\omega)$ $\displaystyle=$
$\displaystyle[\mathbf{I}_{3\times
3}-\epsilon(\mathbf{k})\mathbf{\Pi}(\mathbf{k},\omega)]^{-1}\mathbf{\Pi}(\mathbf{k},\omega),$
(8) $\displaystyle\mathbf{\Pi}(\mathbf{k},\omega)$ $\displaystyle=$
$\displaystyle\mathbf{N}_{11}(\omega)+\epsilon(\mathbf{k})\mathbf{\Gamma}(\mathbf{k},\omega),$
(9) $\displaystyle\mathbf{\Gamma}(\mathbf{k},\omega)$ $\displaystyle=$
$\displaystyle\mathbf{N}_{12}(\omega)[\mathbf{I}_{3\times
3}-\epsilon(\mathbf{k})\mathbf{N}_{22}(\omega)]^{-1}\mathbf{N}_{21}(\omega),$
(10)

where $\mathbf{I}_{3\times 3}$ represents $3\times 3$ identity matrix and
matrices $\mathbf{N}_{mn}$ are defined as

$\displaystyle\begin{pmatrix}\mathbf{N}_{11}^{\alpha,\beta}(\omega)&\mathbf{N}_{12}^{\alpha,\beta}(\omega)\\\
\mathbf{N}_{21}^{\alpha,\beta}(\omega)&\mathbf{N}_{22}^{\alpha,\beta}(\omega)\end{pmatrix}$
(11) $\displaystyle=$
$\displaystyle\sum_{\mu\mu^{\prime}}\frac{D_{\mu\mu^{\prime}}}{\omega+E_{\mu}-E_{\mu^{\prime}}}\begin{pmatrix}c_{\mu\mu^{\prime}}^{\alpha}d_{\mu^{\prime}\mu}^{\beta}&c_{\mu\mu^{\prime}}^{\alpha}c_{\mu^{\prime}\mu}^{\beta}\\\
d_{\mu\mu^{\prime}}^{\alpha}d_{\mu^{\prime}\mu}^{\beta}&d_{\mu\mu^{\prime}}^{\alpha}c_{\mu^{\prime}\mu}^{\beta}\end{pmatrix}.$

Usually we have $\mathbf{N}_{12}=\mathbf{N}_{21}=0$, so the Green’s function
matrix is simplified to be

$\displaystyle\mathbf{G}(\mathbf{k},\omega)=\frac{\mathbf{N}_{11}(\omega)}{\mathbf{I}_{3\times
3}-\epsilon(\mathbf{k})\mathbf{N}_{11}(\omega)}.$ (12)

In the MI-SF phase transition process, we assume the one-particle excitation
is described by operator $\hat{A}_{i}^{\dagger}=\sum_{\sigma=\pm
1,0}\chi_{\sigma}\hat{a}_{i\sigma}^{\dagger}$ with
$\sum_{\sigma}|\chi_{\sigma}|^{2}=1$, which indicates that the superfluid
order parameter can be expressed as $\psi_{\sigma}=\sqrt{n_{s}}\chi_{\sigma}$
($\sigma=0,\pm 1$ and $n_{s}$ is the number of condensate atoms in each site).
Then the most probable excitation is determined by maximizing
$|\langle\langle[A,A^{\dagger}]\rangle\rangle|_{\mathbf{k}=\mathbf{0},\omega=0}=\sum_{\sigma_{1}\sigma_{2}}\chi_{\sigma_{1}}^{*}\chi_{\sigma_{2}}|G^{\sigma_{1}\sigma_{2}}(\mathbf{k}=\mathbf{0},\omega=0)|$,
this is equivalent to diagonalizing matrix
$\mathbf{G}(\mathbf{k}=\mathbf{0},\omega=0)$ (or matrix
$\mathbf{N}_{11}(\omega=0)$) and keeping the eigenvector with maximum absolute
eigenvalue. The phase boundary is determined by the poles of the diagonalized
Green’s function matrix.

###  III.1 Ferromagnetic case ($U_{2}<0$)

For ferromagnetic case ($U_{2}<0$), theorem in Ref. Katsura and Tasaki (2013)
shows that the ground state universally exhibits saturated ferromagnetic in
superfluid and Mott insulator phases, and the same result is predicted in
numerical calculations with mean-field approximation Pai et al. (2008) and
1D-QMC simulations Batrouni et al. (2009). In this subsection, we will give
detailed analytical analysis of this behavior with SBO method.

The ground state in the MI phase is $|n,m_{S};n\rangle$, which is degenerate
for $m_{s}=-n,-n+1,...,n$. Without loss of generality, we assume that system
in MI phase is in superposition state
$|\psi_{0}\rangle=\sum_{s=-n}^{n}c_{s}|n,s;n\rangle$
($\sum_{s=-n}^{n}|c_{s}|^{2}=1$) with $D_{\psi_{0}}=1$. For the convenience of
following discussions, we also define

$\displaystyle A^{(\alpha)s,s+\alpha}_{n,n\pm 1}$ $\displaystyle=$
$\displaystyle\langle n\pm
1,s+\alpha;n+1|\hat{a}_{\alpha}^{\dagger}|n,s;n\rangle,$ (13) $\displaystyle
B^{(\alpha)s,s-\alpha}_{n,n-1}$ $\displaystyle=$ $\displaystyle\langle
n-1,s-\alpha;n-1|\hat{a}_{\alpha}|n,s;n\rangle,$ (14)

and

$\displaystyle\Omega_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{E^{(0)}(n-1,n+1)-E^{(0)}(n,n)},$ (15)
$\displaystyle\Omega_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{E^{(0)}(n+1,n+1)-E^{(0)}(n,n)},$ (16)
$\displaystyle\Omega_{3}$ $\displaystyle=$
$\displaystyle\frac{1}{E^{(0)}(n-1,n-1)-E^{(0)}(n,n)},$ (17)

then matrix $\mathbf{N}_{11}(0)$ can be derived

$\displaystyle\mathbf{N}_{11}^{0,0}$ $\displaystyle=$
$\displaystyle\sum_{s}|c_{s}|^{2}\Big{[}\Omega_{1}(A^{(0)s,s}_{n,n-1})^{2}+\Omega_{2}(A^{(0)s,s}_{n,n+1})^{2}$
(18) $\displaystyle+\Omega_{3}(B^{(0)s,s}_{n,n-1})^{2}\Big{]},$
$\displaystyle\mathbf{N}_{11}^{1,1}$ $\displaystyle=$
$\displaystyle\sum_{s}|c_{s}|^{2}\Big{[}\Omega_{1}(A^{(1)s,s+1}_{n,n-1})^{2}+\Omega_{2}(A^{(1)s,s+1}_{n,n+1})^{2}$
(19) $\displaystyle+\Omega_{3}(B^{(1)s,s-1}_{n,n-1})^{2}\Big{]},$
$\displaystyle\mathbf{N}_{11}^{-1,-1}$ $\displaystyle=$
$\displaystyle\sum_{s}|c_{s}|^{2}\Big{[}\Omega_{1}(A^{(-1)s,s-1}_{n,n-1})^{2}+\Omega_{2}(A^{(-1)s,s-1}_{n,n+1})^{2}$
(20) $\displaystyle+\Omega_{3}(B^{(-1)s,s+1}_{n,n-1})^{2}\Big{]},$
$\displaystyle\mathbf{N}_{11}^{0,1}=\sum_{s}c_{s}c_{s+1}^{*}\Big{[}\Omega_{1}A^{(1)s,s+1}_{n,n-1}A^{(0)s+1,s+1}_{n,n-1}$
$\displaystyle+\Omega_{2}A^{(1)s,s+1}_{n,n+1}A^{(0)s+1,s+1}_{n,n+1}+\Omega_{3}B^{(1)s+1,s}_{n,n-1}B^{(0)s,s}_{n,n-1}\Big{]},$
$\displaystyle\mathbf{N}_{11}^{0,-1}=\sum_{s}c_{s}c_{s-1}^{*}\Big{[}\Omega_{1}A^{(-1)s,s-1}_{n,n-1}A^{(0)s-1,s-1}_{n,n-1}$
$\displaystyle+\Omega_{2}A^{(-1)s,s-1}_{n,n+1}A^{(0)s-1,s-1}_{n,n+1}+\Omega_{3}B^{(-1)s-1,s}_{n,n-1}B^{(0)s,s}_{n,n-1}\Big{]},$
$\displaystyle\mathbf{N}_{11}^{1,-1}=\sum_{s}c_{s+1}c_{s-1}^{*}\Big{[}\Omega_{1}A^{(-1)s+1,s}_{n,n-1}A^{(1)s-1,s}_{n,n-1}$
$\displaystyle+\Omega_{2}A^{(-1)s+1,s}_{n,n+1}A^{(1)s-1,s}_{n,n+1}+\Omega_{3}B^{(-1)s-1,s}_{n,n-1}B^{(1)s+1,s}_{n,n-1}\Big{]},$

and

$\displaystyle\mathbf{N}_{11}^{1,0}$ $\displaystyle=$
$\displaystyle(\mathbf{N}_{11}^{0,1})^{*},$ (24)
$\displaystyle\mathbf{N}_{11}^{-1,0}$ $\displaystyle=$
$\displaystyle(\mathbf{N}_{11}^{0,-1})^{*},$ (25)
$\displaystyle\mathbf{N}_{11}^{-1,1}$ $\displaystyle=$
$\displaystyle(\mathbf{N}_{11}^{1,-1})^{*}.$ (26)

Diagonalizing the matrix $\mathbf{N}_{11}$ directly is a tough work,
especially when the filling number $n$ is very large. So here we come up with
a more effective method to solve the problem.

When we set $c_{n}=1$ (or $c_{-n}=1$) in the ground state $|\psi_{0}\rangle$,
matrix $\mathbf{N}_{11}$ is diagonal with maximum eigenvalue
$\lambda_{m}=\Omega_{2}(A^{(1)n,n+1}_{n,n+1})^{2}+\Omega_{3}(B^{(1)n,n-1}_{n,n-1})^{2}=(n+1)\Omega_{2}+n\Omega_{3}$.
Then we prove that the matrix $\lambda_{m}\mathbf{I}_{3\times
3}-\mathbf{N}_{11}$ can be expanded as summation of a series of positive-
semidefinite matrices, the details of the provement are presented in
appendix.B. Thus the eigenvalues $\lambda$ of matrix $\mathbf{N}_{11}$ must
satisfy the condition $\lambda\leq\lambda_{m}$, and the conditions for
$\lambda=\lambda_{m}$ are given as

$\displaystyle\frac{P^{(\alpha_{1})}_{B}(s)}{P^{(\alpha_{2})}_{B}(s)}$
$\displaystyle=$ $\displaystyle R_{\alpha_{1},\alpha_{2}},$ (27)
$\displaystyle R_{1,0}R_{-1,0}=\frac{1}{2},$ (28)

where $R_{\alpha_{1},\alpha_{2}}$ is s-independent and we have defined

$\displaystyle P^{(\alpha)}_{A}(s)$ $\displaystyle=$ $\displaystyle
c_{s-\alpha}A^{(\alpha)s-\alpha,s}_{n,n-1},$ (29) $\displaystyle
P^{(\alpha)}_{B}(s)$ $\displaystyle=$ $\displaystyle
c_{s+\alpha}B^{(\alpha)s+\alpha,s}_{n,n-1}.$ (30)

In the ground state $|\psi_{0}\rangle$ (Mott-insulator phase), we have
$\langle\vec{S}\rangle^{2}=n^{2}$ (see appendix.B), so the MI phase is in
ferromagnetic state. Meanwhile the excitation described with the eigenvector
$(1,R_{1,0},R_{-1,0})$ indicates that the superfluid phase involved in the
phase transition process is also in ferromagnetic state (because
$1-2R_{1,0}R_{-1,0}=0$). Therefore the system is ferromagnetic in both
superfluid and Mott-insulator phases.

The phase boundary for ferromagnetic case is given by

$\displaystyle zt_{c}[(n+1)\Omega_{2}+n\Omega_{3}]=1.$ (31)

###  III.2 Antiferromagnetic case ($U_{2}>0$)

For antiferromagnetic case ($U_{2}>0$), perturbative mean-field approximation
(PMFA) Tsuchiya et al. (2004), density-matrix renormalization group (DMRG)
Rizzi et al. (2005) and 1D-QMC Batrouni et al. (2009); Apaja and Syljuåsen
(2006) studies have shown that the MI phase for even filling lobe is
considerably more stable against the superfluid phase than that for odd
filling lobe. The MI phase for both even and odd filling lobes is transited to
polar superfluid phase with zero global magnetization. In this subsection, we
use SBO method to analytically analyze the MI-SF transition process for
antiferromagnetic case, of which the ground state energy in the MI phase is
minimized by taking minimum total spin $S_{i}$ in Eq. (3).

For odd filling number case, the ground state is $|1,m_{S};n\rangle$, which is
degenerate for $m_{S}=\pm 1,0$, the system in MI phase can be in any
superposition states of $|1,m_{S};n\rangle$. Similarly, we assume the system
is in state $|\psi_{0}\rangle=\sum_{s=\pm 1,0}c_{s}|1,s;n\rangle$
($\sum_{s=\pm 1,0}|c_{s}|^{2}=1$) with $D_{\psi_{0}}=1$, then we calculate the
corresponding matrix $\mathbf{N}_{11}(0)$ (see appendix.B). Diagonalize matrix
$\mathbf{N}_{11}(0)$ and the resulting eigenvalues are

$\displaystyle\lambda_{0}=3(\Upsilon_{2}+\Upsilon_{4}),$ (32)
$\displaystyle\lambda_{\pm}=\frac{1}{2}(\Upsilon_{1}+\Upsilon_{3}+7\Upsilon_{2}+7\Upsilon_{4})$
$\displaystyle\pm\frac{1}{2}\sqrt{(K_{1}-K_{2})^{2}+4K_{1}K_{2}|c_{0}^{2}-2c_{1}c_{-1}|^{2}},$
(33)

where $K_{1}=3\Upsilon_{2}+\Upsilon_{3}-2\Upsilon_{4}$ and
$K_{2}=\Upsilon_{1}-2\Upsilon_{2}+3\Upsilon_{4}$ with $\Upsilon_{m}$ defined
as

$\displaystyle\Upsilon_{1}$ $\displaystyle=$
$\displaystyle\frac{n+1}{3}\frac{1}{E^{(0)}(0,n+1)-E^{(0)}(1,n)},$ (34)
$\displaystyle\Upsilon_{2}$ $\displaystyle=$
$\displaystyle\frac{n+4}{15}\frac{1}{E^{(0)}(2,n+1)-E^{(0)}(1,n)},$ (35)
$\displaystyle\Upsilon_{3}$ $\displaystyle=$
$\displaystyle\frac{n+2}{3}\frac{1}{E^{(0)}(0,n-1)-E^{(0)}(1,n)},$ (36)
$\displaystyle\Upsilon_{4}$ $\displaystyle=$
$\displaystyle\frac{n-1}{15}\frac{1}{E^{(0)}(2,n-1)-E^{(0)}(1,n)}.$ (37)

The maximum eigenvalue can be obtained by choosing $\lambda_{+}$ and take the
maximum value of $|c_{0}^{2}-2c_{1}c_{-1}|^{2}_{\text{max}}=1$, the
corresponding eigenvector $(\chi_{0},\chi_{1},\chi_{-1})$, which describes the
symmetry of superfluid order in the excitation, satisfies the condition
$|\chi_{0}^{2}-2\chi_{1}\chi_{-1}|^{2}=1$. This result indicates that the
superfluid phase is in polar state, which is in consistent with perturbation
mean-field approximation (PMFA) method Tsuchiya et al. (2004). Then obviously,
the phase boundary for the odd filling number case is given by

$\displaystyle
zt_{c}\left(\Upsilon_{1}+4\Upsilon_{2}+\Upsilon_{3}+4\Upsilon_{4}\right)=1.$
(38)

For even filling number case, the ground state is $|0,0;n\rangle$, so we can
set $D_{|\mu\rangle=|0,0;n\rangle}=1$ if fluctuations are neglected. Then it
can be directly derived that matrix $\mathbf{N}_{11}(0)$ is diagonal with
$\mathbf{N}_{11}^{\alpha,\beta}=N_{0}\delta_{\alpha\beta}$, the excitation
spectrum is degenerate for spin-$\sigma$ ($\sigma=\pm 1,0$) excitations. The
phase boundary is derived by setting $1+zt_{c}N_{0}=0$ and have

$\displaystyle zt_{c}\left(\frac{n/3+1}{\Delta E_{1}}+\frac{n/3}{\Delta
E_{2}}\right)=1,$ (39)

where

$\displaystyle\Delta E_{1}$ $\displaystyle=$ $\displaystyle
E^{(0)}(1,n+1)-E^{(0)}(0,n),$ (40) $\displaystyle\Delta E_{2}$
$\displaystyle=$ $\displaystyle E^{(0)}(1,n-1)-E^{(0)}(0,n).$ (41)

To decide the symmetry of superfluid order parameter in the MI-SF phase
transition process, the PMFA method has taken fourth-order perturbation into
account Tsuchiya et al. (2004). In the SBO method, we can consider the two-
particle excitation, which is described with operator
$(\hat{A}^{\dagger}_{i})^{2}=\sum_{\sigma_{1}\sigma_{2}}\chi_{\sigma_{1}}\chi_{\sigma_{2}}\hat{a}_{i\sigma_{1}}^{\dagger}\hat{a}_{i\sigma_{2}}^{\dagger}$.
Similar approach can be applied to the two-particle excitation process and
gives

$\displaystyle G^{ij}_{\hat{A}^{2}}(t-t^{\prime})$ $\displaystyle=$
$\displaystyle-i\Theta(t-t^{\prime})\langle[\hat{A}_{i}(t)^{2},(\hat{A}_{j}^{\dagger}(t^{\prime}))^{2}]\rangle$
(42) $\displaystyle=$
$\displaystyle\sum_{\mu\mu^{\prime}\nu\nu^{\prime}}G^{ij}_{\mu\mu^{\prime},\nu\nu^{\prime}}(t-t^{\prime})c^{\hat{A}^{2}}_{\mu\mu^{\prime}}d^{\hat{A}^{2}}_{\nu\nu^{\prime}},$

where
$c^{\hat{A}^{2}}_{\mu\mu^{\prime}}=\langle\mu|\hat{A}^{2}|\mu^{\prime}\rangle$
and
$d^{\hat{A}^{2}}_{\nu\nu^{\prime}}=\langle\nu|(\hat{A}^{\dagger})^{2}|\nu^{\prime}\rangle$.
By Taking the conditions $D_{|\mu\rangle=|0,0;n\rangle}=1$ and
$D_{|\mu\rangle\neq|0,0;n\rangle}=0$ as well as the expression in Eq. (6), the
problem is equivalent to maximising absolute value of
$G_{\hat{A}^{2}}(\mathbf{k}=0,\omega=0)$, which can be calculated as

$\displaystyle|G_{\hat{A}^{2}}(\mathbf{k}=\mathbf{0},\omega=0)|=2(\Lambda_{1}+\Lambda_{3})$
$\displaystyle+(\Lambda_{2}+\Lambda_{4}-\frac{2}{3}\Lambda_{1}-\frac{2}{3}\Lambda_{3})|\chi_{0}^{2}-2\chi_{1}\chi_{-1}|^{2},$

where

$\displaystyle\Lambda_{1}$ $\displaystyle=$
$\displaystyle\frac{n(n-2)}{15}\frac{1}{E^{(0)}(2,n-2)-E^{(0)}(0,n)},$ (44)
$\displaystyle\Lambda_{2}$ $\displaystyle=$
$\displaystyle\frac{n(n+1)}{9}\frac{1}{E^{(0)}(0,n-2)-E^{(0)}(0,n)},$ (45)
$\displaystyle\Lambda_{3}$ $\displaystyle=$
$\displaystyle\frac{(n+3)(n+5)}{15}\frac{1}{E^{(0)}(2,n+2)-E^{(0)}(0,n)},$
(46) $\displaystyle\Lambda_{4}$ $\displaystyle=$
$\displaystyle\frac{(n+2)(n+3)}{9}\frac{1}{E^{(0)}(0,n+2)-E^{(0)}(0,n)}.$ (47)

It is clear that
$\Lambda_{2}+\Lambda_{4}-\frac{2}{3}\Lambda_{1}-\frac{2}{3}\Lambda_{3}$ is
positive, then we need to take the maximum value of
$|\chi_{0}^{2}-2\chi_{1}\chi_{-1}|_{\text{max}}=1$, which indicates that the
superfluid phase involved in the excitation is in polar state Tsuchiya et al.
(2004). In physics, we can understand the favorable of polar state in the
phase transition process is due to the hopping of singlet pairs. The forming
of such two-atom dimers in the system can greatly reduce one-particle
excitation in even filling number case, even for ferromagnetic interaction
($U_{2}<0$). This phenomenon is obvious in the SBO phase diagrams shown in the
next subsection.

###  III.3 Phase diagrams of self-consistent SBO method

The phase diagrams of MI-SF transition in spin-1 bosonic system have been
studied with several approaches, including mean-field approximation Tsuchiya
et al. (2004); Kimura et al. (2005); Pai et al. (2008); Svidzinsky and Chui
(2003), QMC Batrouni et al. (2009); de Forges de Parny et al. (2013), DMRG
Rizzi et al. (2005) and strong-coupling expansion method Kimura (2013). In
this subsection, we will give the phase diagram with self-consistent SBO
approach.

In the above two subsections, we have assumed that all atoms stay in the same
ground state in MI phase. However, fluctuations exist in the phase diagram
even at zero temperature, especially near the phase boundary. The fluctuation
term can be calculated with self-consistent SBO approach Sajna et al. (2015),
which uses the spectral theorem Haley and Erdös (1972)

$\displaystyle\langle\hat{L}^{j}_{\nu\nu^{\prime}}(t)\hat{L}^{i}_{\mu\mu^{\prime}}\rangle$
$\displaystyle=$ $\displaystyle\int_{-\infty}^{\infty}\frac{d\omega}{2\pi
i}e^{i\omega
t}f(\omega)\big{[}G^{ij}_{\mu\mu^{\prime},\nu\nu^{\prime}}(\omega-i0^{+})$
(48)
$\displaystyle-G^{ij}_{\mu\mu^{\prime},\nu\nu^{\prime}}(\omega+i0^{+})\big{]},$

where $f(\omega)=1/(e^{\omega/T}-1)$ is the Bose-statistical function (with
$k_{B}=1$).

Let’s take the antiferromagnetic case with odd filling number as an example,
for ground state $|\mu_{0}\rangle=|1,0;n\rangle$ and excited state
$|\mu\rangle=|0,0;n-1\rangle$, we apply similar method introduced in spinless
model Sajna et al. (2015) and have

$\displaystyle D_{\mu}$ $\displaystyle=$
$\displaystyle\frac{1}{N_{s}}\sum_{\mathbf{k}}\int_{-\infty}^{\infty}\frac{d\omega}{2\pi
i}f(\omega)\Big{[}F(\mathbf{k},\omega-i0^{+})$ (49)
$\displaystyle-F(\mathbf{k},\omega+i0^{+})\Big{]},$ $\displaystyle
F(\mathbf{k},\omega)$ $\displaystyle=$ $\displaystyle
G_{\mu_{0}\mu,\mu\mu_{0}}(\mathbf{k},\omega)+\sqrt{\frac{n+1}{n+2}}G_{\mu_{1}\mu_{0},\mu\mu_{0}}(\mathbf{k},\omega)$
(50)
$\displaystyle+2\sqrt{\frac{n+4}{5(n+2)}}G_{\mu_{2}\mu_{0},\mu\mu_{0}}(\mathbf{k},\omega),$

where $|\mu_{1}\rangle=|0,0;n+1\rangle$ and $|\mu_{2}\rangle=|2,0;n+1\rangle$,
$N_{s}$ is the total number of lattice sites.

The Green’s functions in the integration can also be found in the motion
equation in Eq. (6). Then the combined self-consistent equations can be solved
afterwards.

Figure 1: Zero-temperature phase diagram for 2-d antiferromagnetic
($U_{2}=0.1U_{0}$) case. The mean-field (MF) and Monte-Carlo (MC) results are
obtained from Ref. Pai et al. (2008); de Forges de Parny et al. (2013).

In Fig. 1, we give the phase diagrams for antiferromagnetic case with several
methods. Compared to the spinless boson system, we find that the even lobes
grow at the expense of the odd lobes, which would entirely disappear when
$\frac{U_{2}}{U_{0}}\geq 0.5$. And for the even lobes, atoms in each site tend
to form singlet pairs, so the fluctuations in states $|1,m_{S};n\pm 1\rangle$
(for even $n$) have been neglected in the SBO phase diagram.

In this figure, we find that the self-consistent SBO zero-temperature phase
diagram is in quantitative agreement with quantum Monte-Carlo simulation
result de Forges de Parny et al. (2013), which indicates better description of
quantum fluctuation in SBO method than in mean-field approximation Tsuchiya et
al. (2004); Pai et al. (2008).

Figure 2: Zero-temperature phase diagram for 2-d ferromagnetic
($U_{2}=-0.1U_{0}$) case.

For ferromagnetic case, we have pointed out that the system universally
exhibits ferromagnetism in both MI and SF phases. The phase diagrams derived
with self-consistent SBO, mean-field approximation Pai et al. (2008) and
Monte-Carlo de Forges de Parny et al. (2013) methods are presented in Fig. 2.
Similarly, we can conclude that the SBO method gives more precise phase
diagram than mean-field approximation method. For even lobes, the difference
between self-consistent SBO and Monte-Carlo results can be explained as due to
dimer-effect in the system (consider two-particle excitation), which tends to
decrease quantum fluctuations in states $|S,m_{S};n\pm 1\rangle$ (for even
$n$) and induces lower critical hopping energy in the phase diagram.

If external magnetic field is applied to the system, the Hamiltonian is
changed by linear Zeeman effect with
$\hat{H}^{\prime}=\hat{H}-\eta\sum_{i}\hat{S}_{i}^{z}$. The external field-
induced term tends to change both Mott-insulator and superfluid phases in the
phase transition process Svidzinsky and Chui (2003). The result can be derived
with similar procedures introduced for non-external-field case.

At finite temperature, besides the states involved in the excitations from
ground state, we need to further take the fluctuations in other states into
consideration. For example, for the single-atom filling ($n=1$) case of
antiferromagnetic system with external field, the ground state is
$|1,1;1\rangle$, spin-1 excitation takes place in the MI-SF phase transition
process. Then besides the states $|0,0;0\rangle$ and $|2,2;2\rangle$, which
are involved in the spin-1 excitation from $|1,1;1\rangle$, we need further
consider the fluctuations in states $|1,\alpha;1\rangle$ ($\alpha=-1,0$) and
$|0,0;2\rangle$, whereas other states are neglected because the energy
differences are much larger than temperature ($k_{B}=1$). In Fig. 3, we give
the phase diagrams for finite temperature system with self-consistent SBO
method and mean-field approximation Ohliger (2008). By comparison, we find
that the SBO results are much more sensitive to the change of temperature than
mean-field results. Because when temperature increases but is still low
enough, the fluctuations in the spin-1 excitation-involved states can be
enhanced much greater than those in other states in consideration.

Figure 3: Finite temperature phase diagrams for single-atom filling
antiferromagnetic system in external field. Parameters are chosen as
$U_{2}=0.04U_{0}$ and $\eta=0.05U_{0}$.

##  IV insulating state with one atom per site

The particle excitation properties of the spin-1 bosonic system are shown in
Fig. 1 and Fig. 2. However, the inner spin properties in the Mott lobes have
also attracted many attentions, even before the study of cold atom system Chen
and Levy (1973).

In second-order perturbation theory in $t$, the system in the Mott lobe with
one atom per site is given by super-exchange processes, which can be described
with an effective bilinear-biquadratic spin Hamiltonian Imambekov et al.
(2003); Yip (2003); Rodríguez et al. (2011); Chung and Yip (2009)

$\displaystyle\hat{H}_{\text{eff}}$ $\displaystyle=$
$\displaystyle-\sum_{\langle
ij\rangle}\left[J_{1}\vec{S}_{i}\cdot\vec{S}_{j}+J_{2}(\vec{S}_{i}\cdot\vec{S}_{j})^{2}\right]$
(51) $\displaystyle+\sum_{i}\left(\hat{H}_{i}^{L}+\hat{H}_{i}^{Q}\right).$

The absence of higher order terms, such as
$(\vec{S}_{i}\cdot\vec{S}_{j})^{3}$, is due to the fact that the product of
any three spin operators for an $S=1$ particle can be expressed via lower
order terms. The parameters in the Hamiltonian can be calculated with
perturbation theory, which gives $J_{1}=\frac{2t^{2}}{g_{2}}$ and
$J_{2}=\frac{4t^{2}}{3g_{0}}+\frac{2t^{2}}{3g_{2}}$ with
$g_{0,2}=4\pi\hbar^{2}a_{0,2}/m_{a}$ representing the scattering strengths of
spin-0 and spin-2 channels. In the cold-atom system, we assume both of them
are positive. For antiferromagnetic case ($a_{2}>a_{0}$), $J_{2}>J_{1}$, and
for ferromagnetic case ($a_{0}>a_{2}$), $J_{1}>J_{2}$. $\hat{H}_{i}^{L}$ and
$\hat{H}_{i}^{Q}$ in Eq. (51) represent linear Zeeman effect (LZE) and
quadratic Zeeman effect (QZE) respectively when a magnetic field is applied to
the system,

$\displaystyle\hat{H}_{i}^{L}$ $\displaystyle=$
$\displaystyle-\lambda\hat{S}_{i}^{z},$ (52) $\displaystyle\hat{H}_{i}^{Q}$
$\displaystyle=$ $\displaystyle q(\hat{S}_{i}^{z})^{2}.$ (53)

In the linear term we have included the Lagrangian multiplier due to
magnetization conservation Stenger et al. (1999). In our following discussion,
we assume $q>0$, and the result is symmetrical for $\pm\lambda$.

We find the Hamiltonian in Eq. (51) is similar to that of spin-1 Heisenberg
model with biquadratic exchange, which has been studied in Ref. Chen and Levy
(1973); Micnas (1976); Papanicolaou (1988). Then to solve the system with
standard basis operator (SBO) method, we expand the Hamiltonian in Eq. (51) on
the basis $|m\rangle=|1,m;1\rangle$ and obtain the expression in appendix A
(Eq. (A.1)) with

$\displaystyle V^{i}_{mm^{\prime}}$ $\displaystyle=$
$\displaystyle\delta_{mm^{\prime}}V_{m}=\delta_{mm^{\prime}}(-\lambda
m+qm^{2}),$ (54)

and

$\displaystyle T^{ij}_{mm^{\prime}nn^{\prime}}$ $\displaystyle=$
$\displaystyle-J_{1}\langle m|_{i}\langle
n|_{j}\vec{S}_{i}\cdot\vec{S}_{j}|m^{\prime}\rangle_{i}|n^{\prime}\rangle_{j}$
(55) $\displaystyle-J_{2}\langle m|_{i}\langle
n|_{j}(\vec{S}_{i}\cdot\vec{S}_{j})^{2}|m^{\prime}\rangle_{i}|n^{\prime}\rangle_{j}.$

Substitute the expressions in Eq. (54) and Eq. (55) into the motion equation
in the appendix (Eq. (A.5)), we derive

$\displaystyle(\omega+E_{m}-E_{m^{\prime}})G_{mm^{\prime},nn^{\prime}}(\mathbf{k},\omega)-\delta_{mn^{\prime}}\delta_{nm^{\prime}}D_{mm^{\prime}}$
(56) $\displaystyle=$
$\displaystyle-\zeta(\mathbf{k})D_{mm^{\prime}}\sum_{rr^{\prime}}T_{m^{\prime}mrr^{\prime}}G_{rr^{\prime},nn^{\prime}}(\mathbf{k},\omega),$

where we have defined
$\zeta(\mathbf{k})=-2\sum_{s=1}^{\text{dim}}\cos{k_{s}l}$ and
$E_{m}=V_{m}+z\sum_{n}D_{n}T_{mmnn}$.

In our system, the lattice is initially deep enough that the hopping energy
can be neglected ($t=J_{1}=J_{2}=0$), so the ground state, which can be
$|0\rangle$ (nematic) or $|\pm 1\rangle$ (ferromagnetic), is determined by
$\lambda$ and $q$, and the boundary is given by $q=|\lambda|$. Then increase
the hopping energy by reducing the lattice depth, we have $J_{1}>0$ and
$J_{2}>0$, the system will be excited from ground state, the spin-excitation
can be described with SBO method.

When $\lambda>0$ and $\lambda<q$, the ground state is $|0\rangle$, and it will
be excited to state $|1\rangle$ first, this process is described with Green’s
function $G_{01,nn^{\prime}}(\mathbf{k},\omega)$. Based on the motion equation
in Eq. (56), we can write the combined equations related to
$G_{01,10}(\mathbf{k},\omega)$

$\displaystyle\mathbf{M}\begin{pmatrix}G_{01,10}(\mathbf{k},\omega)\\\
G_{-10,10}(\mathbf{k},\omega)\end{pmatrix}=\begin{pmatrix}D_{01}\\\
0\end{pmatrix},$ (57)

with

$\displaystyle\mathbf{M}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}-J_{1}D_{01}\zeta(\mathbf{k})&(J_{2}-J_{1})D_{01}\zeta(\mathbf{k})\\\
(J_{2}-J_{1})D_{-10}\zeta(\mathbf{k})&-J_{1}D_{-10}\zeta(\mathbf{k})\end{pmatrix}$
(58) $\displaystyle+\begin{pmatrix}\omega+\Delta E_{1}&0\\\ 0&\omega+\Delta
E_{2}\end{pmatrix},$

where $\Delta E_{1}=E_{0}-E_{1}$ and $\Delta E_{2}=E_{-1}-E_{0}$.

The excitation spectrum is determined by the solution of equation
$\det{\mathbf{M}}=0$. If we neglect the fluctuations and assume $D_{0}=1$ and
$D_{\pm 1}=0$, the equation corresponds to

$\displaystyle(\omega+\lambda)^{2}$ (59) $\displaystyle=$
$\displaystyle(q+zJ_{2}+J_{1}\zeta(\mathbf{k}))^{2}-(J_{2}-J_{1})^{2}\zeta(\mathbf{k})^{2}.$

We find that in Eq. (59), if $q$ is sufficiently large, the solution of
$\omega$ in the equation is always non-zero, leading to positive excitation
gap in the excitation, then the phase transition in the system is denied and
the the system is stable in state $|0\rangle$ (nematic state). To determine
the boundary, we analyze the minimum value of the expression on the right side
of Eq. (59) (see appendix.C).

For antiferromagnetic case ($J_{2}>J_{1}$), we have
$q>0>zJ_{2}-z\frac{J_{2}^{2}}{J_{1}}$, so the minimum (we label it as Min) is
located at the point $\mathbf{k}=\mathbf{0}$. If $\lambda^{2}<\text{Min}$, the
excitation gap turns to be positive, thus the excitation is denied. So the
boundary is given by

$\displaystyle\lambda^{2}=(q+zJ_{2}-zJ_{1})^{2}-z^{2}(J_{2}-J_{1})^{2}.$ (60)

While for ferromagnetic case ($J_{1}>J_{2}$), if $q\geq
zJ_{2}-z\frac{J_{2}^{2}}{J_{1}}$, the minimum is located at the point
$\mathbf{k}=\mathbf{0}$, and the phase boundary is also given by Eq. (60). If
$0<q<zJ_{2}-z\frac{J_{2}^{2}}{J_{1}}$, the minimum seems to be located at the
global minimum point with $\mathbf{k}\neq\mathbf{0}$. However, we note that,
for any $q=Q_{0}>0$, when we increase hopping energy $t$ from 0 to positive,
it will always pass the point $2q_{0}$ before reaching $q_{1}$ (see Fig.
4(b)), so the excitations happen before $q<zJ_{2}-zJ_{2}^{2}/J_{1}$, we can
still take $\mathbf{k}=0$ in calculating the excitation gap. So in this case,
the boundary of large-$q$ (nematic) phase is still given by Eq. (60).

Figure 4: Magnetic phase diagrams in $n=1$ Mott lobe for antiferromagnetic
case (a) and ferromagnetic case (b). The bold black line on the $q$-axis
represents XY-FM phase Rodríguez et al. (2011), and the black arrow along $-q$
direction indicates the change of the relative location for $q=Q_{0}$ when we
increase the hopping energy $t$ from zero to finite. In the figure we have
defined $q_{0}=z|J_{2}-J_{1}|$ and $q_{1}=z(J_{2}-J_{2}^{2}/J_{1})$. We note
that $q_{1}<q_{0}$ in the ferromagnetic case ($J_{1}>J_{2}$).

Similarly, when $\lambda>q>0$, we start from $|1\rangle$ and consider the
excitation to $|0\rangle$ (excite to $|0\rangle$ first), the excitation
spectrum can be solved as
$\omega_{\mathbf{k}}=\lambda-q+zJ_{1}+J_{1}\zeta(\mathbf{k})$, which is always
positive.

In Fig. 4, we divide the $\lambda$-$q$ plane (with solid blue curves) into
several regimes, corresponding to several different phases. Each phase
represents the ground state and excitation properties in the regime when we
increase hopping energy $t$ from $0$ to finite. For example, the nematic
(ferromagnetic) phase indicates that the ground state is $|0\rangle$ ($|\pm
1\rangle$) and the spin excitation is denied when $t$ increases, while the
partially magnetic phase indicates that the system initially in state
$|0\rangle$ can be excited to state $|\pm 1\rangle$ when $t$ increases.

We can take a look at the results in some special cases. For example, when
$\lambda=q=0$, states $|\pm 1\rangle$ are symmetric by rotation, so we have
$D_{1}=D_{-1}$. Then we can solve the excitation spectrum from $|0\rangle$ to
$|\pm 1\rangle$ as

$\displaystyle\omega_{\mathbf{k}}=D_{01}\sqrt{(zJ_{2}+J_{1}\zeta(\mathbf{k}))^{2}-(J_{2}-J_{1})^{2}\zeta(\mathbf{k})^{2}},$
(61)

the result in Eq. (61) is the same as that derived with systematic $1/n$
expansion method Papanicolaou (1988) besides an additional fluctuation-
correction term $D_{01}$.

When $\lambda=0$ and $q>0$, states $|\pm 1\rangle$ are still symmetric by
rotation, and the PM phase in Fig. 4 turns to be XY-FM phase Rodríguez et al.
(2011), which fulfills $\langle S_{i}^{z}\rangle=0$ but presents a non-zero
transversal magnetization.

From Fig. 4, we read out that the boundary is $q=2z(J_{1}-J_{2})$ for
ferromagnetic case, more detailed calculation including the fluctuations of
occupation probabilities gives the boundary as

$\displaystyle q_{c}=2z(D_{0}-D_{1})(J_{1}-J_{2}).$ (62)

Then $D_{0}$ and $D_{1}$ are determined by the self-consistent equations

$\displaystyle D_{1}$ $\displaystyle=$ $\displaystyle
D_{01}\frac{1}{N_{s}}\sum_{\mathbf{k}}\left[F_{+}(\mathbf{k})+F_{-}(\mathbf{k})\right],$
(63) $\displaystyle D_{0}$ $\displaystyle=$ $\displaystyle 1-2D_{1},$ (64)

with

$\displaystyle
F_{\pm}(\mathbf{k})=\frac{\pm\omega_{\mathbf{k},q}+q+zJ_{2}D_{01}+J_{2}D_{01}\zeta(\mathbf{k})}{2\omega_{\mathbf{k},q}}f(\pm\omega_{\mathbf{k},q}),$

where $\omega_{\mathbf{k},q}$ is defined as

$\displaystyle\omega_{\mathbf{k},q}$ $\displaystyle=$
$\displaystyle\sqrt{[q+D_{01}(J_{1}\zeta(\mathbf{k})+zJ_{2})]^{2}-(J_{2}-J_{1})^{2}D_{01}^{2}\zeta(\mathbf{k})^{2}}.$

In Fig. 5, we give the analytical and SC-SBO phase separation diagrams as well
as the result estimated in Ref. Rodríguez et al. (2011). We find that when
$J_{1}\approx J_{2}$ ($\theta\approx-0.75\pi$), the results of SBO method and
Ref. Rodríguez et al. (2011) are very close.

Figure 5: Phase separation in the case $\lambda=0$. The solid curves represent
analytical solutions without fluctuations, the dot-dashed curves represent SC-
SBO results and the dot-segment-dashed curves represent the results calculated
in Ref. Rodríguez et al. (2011). Here we have defined $J_{1}=-J\cos{\theta}$,
$J_{2}=-J\sin{\theta}$ with $J=\sqrt{J_{1}^{2}+J_{2}^{2}}$.

##  V insulating state with two atoms per site

For insulating state with two atoms per site, the spin value in each site can
be $S=0$ or $S=2$. For antiferromagnetic case ($U_{2}>0$), in the limit $t=0$,
the on-site energy is minimized when $S=0$, corresponding to the singlet
state. But when $t$ is of order $\sqrt{U_{0}U_{2}}$, we need to consider the
excitation to $S=2$, of which the exchange energy is of order $t^{2}/U_{0}$.

Similar to insulating state with one atom case, we derive the effective
Hamiltonian with perturbation theory ($t\ll U_{0}$), which can be expressed as

$\displaystyle\hat{H}_{\text{eff},2}$ $\displaystyle=$
$\displaystyle\sum_{\langle
ij\rangle}\sum_{\alpha\alpha^{\prime}\beta\beta^{\prime}}H^{ij}_{\alpha\alpha^{\prime},\beta\beta^{\prime}}\hat{L}^{i}_{\alpha\alpha^{\prime}}\hat{L}^{j}_{\beta\beta^{\prime}},$
(67)

where $|\alpha\rangle$-$|\beta^{\prime}\rangle$ belong to the sets $\\{S=0\\}$
and $\\{S=2,S_{z}=-2,...,2\\}$.
$H_{\alpha\alpha^{\prime},\beta\beta^{\prime}}$ can be derived with
perturbation theory Imambekov et al. (2003). For simplicity, in the following
consideration in this section, we label $|s\rangle=|0,0;2\rangle$ and
$|m\rangle=|2,m;2\rangle$. In appexdix.D, we give detailed expressions for
$H_{\alpha\alpha^{\prime},\beta\beta^{\prime}}$.

When $t=0$, the system is initially in the singlet state $|s\rangle$, then we
need to consider the $|s\rangle$ to $|m\rangle$ ($m=\pm 2,\pm 1,0$)
excitations when $t$ is nonzero. We find the excitation spectrums are the same
for different $m$ due to global SU(2) symmetry, so here we consider
$|s\rangle$ to $|0\rangle$ excitation. Neglecting quantum fluctuations, the
equations connecting to Green’s functions $G_{s0,0s}(\mathbf{k},\omega)$ can
be written as

$\displaystyle\mathbf{M}_{s\rightarrow
0}\begin{pmatrix}G_{s0,0s}(\mathbf{k},\omega)\\\
G_{0s,0s}(\mathbf{k},\omega)\end{pmatrix}=\begin{pmatrix}D_{s0}\\\
0\end{pmatrix},$ (68)

with

$\displaystyle\mathbf{M}_{s\rightarrow 0}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}\zeta(\mathbf{k})H_{0s,s0}&\zeta(\mathbf{k})H_{0s,0s}\\\
-\zeta(\mathbf{k})H_{s0,s0}&-\zeta(\mathbf{k})H_{s0,0s}\end{pmatrix}$ (69)
$\displaystyle+\begin{pmatrix}\omega+E_{s}-E_{0}&0\\\
0&\omega+E_{0}-E_{s}\end{pmatrix},$

where $E_{s}=zH_{ss,ss}$ and $E_{0}=zH_{00,ss}$. Then the excitation spectrum
is given by

$\displaystyle\omega(\mathbf{k})$ $\displaystyle=$
$\displaystyle\sqrt{9U_{2}^{2}+16\frac{U_{2}t^{2}}{U_{0}}\zeta(\mathbf{k})},$
(70)

which is the same as that derived in Imambekov et al. (2003). And by setting
excitation gap to zero, we obtain the critical hopping energy as
$t_{c}=\frac{3}{4}\sqrt{\frac{U_{0}U_{2}}{z}}$.

If we further include linear Zeeman effect (LZE) when external magnetic field
is applied to the system Imambekov et al. (2004), the effective Hamiltonian is
updated to be

$\displaystyle\hat{H}_{\text{eff},2}=\sum_{i\alpha}V^{i}_{\alpha}\hat{L}^{i}_{\alpha\alpha}+\sum_{\langle
ij\rangle}\sum_{\alpha\alpha^{\prime}\beta\beta^{\prime}}H^{ij}_{\alpha\alpha^{\prime},\beta\beta^{\prime}}\hat{L}^{i}_{\alpha\alpha^{\prime}}\hat{L}^{j}_{\beta\beta^{\prime}},$
(71)

with $V^{i}_{\alpha}=-\lambda S^{z}_{i\alpha}$. The result is symmetric for
$\pm\lambda$, thus we restrict $\lambda>0$ in the following consideration.

In the limit of zero hopping energy $t=0$, the ground states can be
$|s\rangle$ (singlet phase) when $\lambda<3/2U_{2}$ or $|2\rangle$
(ferromagnetic phase) when $\lambda>3/2U_{2}$. Start from singlet phase, the
system can be excited to nematic phase when the hopping energy is larger than
a critical value. However, when $\lambda$ is large enough, the system is
stable in FM phase without transition to other phases.

If $\lambda<3/2U_{2}$, and start from singlet phase, the system excite to
state $|2\rangle$ first. Neglecting quantum fluctuations and assume
$D_{\mu}=0$ for $\mu\neq s$, the critical hopping energy can be solved (see
appendix.D) as

$\displaystyle t_{c}=\sqrt{\frac{U_{0}}{16z}(9U_{2}^{2}-4\lambda^{2})}.$ (72)

Similarly, start from FM phase, the system will be excited to state
$|s\rangle$ and $|0\rangle$ first, and the boundary of the excitation without
fluctuations is given by

$\displaystyle\frac{zt^{2}}{U_{0}}=\frac{3\lambda
U_{2}-2\lambda^{2}}{8(\lambda-U_{2})}.$ (73)

Figure 6: Magnetic phase diagram in the Mott lobe with two atoms per site. The
solid curve represent analytical result without quantum fluctuations, and the
dashed curve represents the self-consistent result for singlet to nematic
transition including quantum fluctuation. We have defined
$t_{0}=\sqrt{U_{0}U_{2}/z}$ in the figure.

In Fig. 6, we give the magnetic phase diagram for insulating phase with two
atoms per site with LZE. The phase diagram can be interpreted as that when the
external magnetic field is small, the system can be singlet or nematic phase,
and there is also a small ferromagnetic component along the axis of the
magnetic field, and as magnetic field goes up, the magnetization also
increases and finally saturates in the FM regime. The nematic phase in the
phase diagram is referred as canted nematic as the magnetization changes
continuously in the phase Imambekov et al. (2004). We also note that for
$\lambda=0$ in Fig. (6), the self-consistent result give the critical hopping
energy as $t_{c}\approx 0.7t_{0}$, which is very close to the global
transition point given by mean-field Imambekov et al. (2003) and Monte-Carlo
simulation de Forges de Parny et al. (2013).

##  VI Conclusion

In this paper, we solve spin-1 bosonic system trapped in optical lattice with
extended standard basis operator method. Both the MI-SF transition and
magnetic phase diagrams in the insulating states are studied.

For ferromagnetic interacting case ($U_{2}<0$), we analytically calculate that
both superfluid and Mott-insulator phases in the phase diagram are
ferromagnetic. While for antiferromagnetic case ($U_{2}>0$), the superfluid
phase is in polar state and the global magnetization in the whole phase
diagram is zero.

Then we take quantum and thermal fluctuations of occupation probabilities in
the system into account with self-consistent SBO approach and derive self-
consistent SBO phase diagrams, which are in relevant agreement with Monte-
Carlo simulation results at zero temperature. At finite temperature, we also
compare the self-consistent SBO phase diagrams with mean-field phase diagrams
and conclude that the finite-temperature-effect has been underestimated in
mean-field approximation.

For $n=1$ Mott lobe with linear and quadratic Zeeman effect, we describe the
system with an effective bilinear-biquadratic spin Hamiltonian. By considering
the spin-excitations in this system with SBO method, we can decide the spin
phase separations in the magnetic phase diagrams. For antiferromagnetic case,
only nematic and ferromagnetic phases are expected, while for ferromagnetic
case, another partial-nematic phase exists between nematic and ferromagnetic
phases.

For $n=2$ Mott lobe with linear Zeeman effect, we apply similar approach and
derive the magnetic phase diagram consisting of singlet, ferromagnetic and
canted-nematic phases.

Besides the spin-1 system, the SBO method used in this paper can be extended
to solve arbitrary spin-$S$ systems with similar procedures. However, there
are still some cases, such as system in one dimensional lattice with broken
translational symmetry in a regime, where the SBO method fails to explain the
phenomenon.

##  VII appendix

###  VII.1 motion equation of Green’s function

In this section, we will show how to derive the motions equation in SBO
method. Then we will also calculate the Green’s function matrix involved in
Sec. III with the motion equations.

For an arbitrary Hamiltonian with the form

$\displaystyle\hat{H}=\sum_{i\mu\mu^{\prime}}V_{\mu\mu^{\prime}}^{i}\hat{L}^{i}_{\mu\mu^{\prime}}+\sum_{\langle
ij\rangle}\sum_{\mu\mu^{\prime}\nu\nu^{\prime}}T^{ij}_{\mu\mu^{\prime}\nu\nu^{\prime}}\hat{L}^{i}_{\mu\mu^{\prime}}\hat{L}^{j}_{\nu\nu^{\prime}},$
(A.1)

we can give the motion equation for Green’s function
$G^{ij}_{\mu\mu^{\prime},\nu\nu^{\prime}}(t-t^{\prime})=-i\Theta(t-t^{\prime})\langle[\hat{L}^{i}_{\mu\mu^{\prime}}(t),\hat{L}^{j}_{\nu\nu^{\prime}}(t^{\prime})]\rangle$.

Differentiating $G^{ij}_{\mu\mu^{\prime},\nu\nu^{\prime}}(t-t^{\prime})$ with
respect to $t$ and Fourier transforming into $\omega$, we obtain the motion
equation

$\displaystyle\omega
G^{ij}_{\mu\mu^{\prime},\nu\nu^{\prime}}(\omega)=\langle[\hat{L}^{i}_{\mu\mu^{\prime}},\hat{L}^{j}_{\nu\nu^{\prime}}]\rangle+\langle\langle[\hat{L}^{i}_{\mu\mu^{\prime}},H]|\hat{L}^{j}_{\nu\nu^{\prime}}\rangle\rangle_{\omega}.$

The last term in Eq. (VII.1) involves higher-order Green’s function. For
further calculation, we need to decouple these terms with two-operator Green’s
functions, so we take random-phase approximation (RPA) Haley and Erdös (1972),
with which the three-operator Green’s function can be expanded as

$\displaystyle\langle\langle\hat{L}^{n}_{\mu\mu^{\prime}}\hat{L}^{l}_{\xi\xi^{\prime}}|\hat{L}^{m}_{\nu\nu^{\prime}}\rangle\rangle$
(A.3) $\displaystyle=$
$\displaystyle\delta_{\mu\mu^{\prime}}\langle\hat{L}^{n}_{\mu\mu}\rangle
G^{lm}_{\xi\xi^{\prime},\nu\nu^{\prime}}(\omega)+\delta_{\xi\xi^{\prime}}\langle\hat{L}^{l}_{\xi\xi}\rangle
G^{nm}_{\mu\mu^{\prime},\nu\nu^{\prime}}(\omega)$
$\displaystyle+\delta^{nl}[\delta_{\xi\mu^{\prime}}G^{nm}_{\mu\xi^{\prime},\nu\nu^{\prime}}(\omega)-\delta_{\mu\xi^{\prime}}G^{nm}_{\xi\mu^{\prime},\nu\nu^{\prime}}(\omega)],$

then combined with the equation

$\displaystyle[\hat{L}^{m}_{\mu\mu^{\prime}},\hat{L}^{n}_{\nu\nu^{\prime}}]=\delta^{mn}(\delta_{\nu\mu^{\prime}}\hat{L}^{n}_{\mu\nu^{\prime}}-\delta_{\mu\nu^{\prime}}\hat{L}^{n}_{\nu\mu^{\prime}}),$
(A.4)

The motion equation in Eq. (VII.1) can be rewritten as

$\displaystyle\omega G^{ij}_{\mu\mu^{\prime},\nu\nu^{\prime}}(\omega)$
$\displaystyle=$
$\displaystyle\delta^{ij}\delta_{\mu\nu^{\prime}}\delta_{\nu\mu^{\prime}}D_{\mu\mu^{\prime}}^{i}+\sum_{l}\sum_{\xi\xi^{\prime}}D_{\mu\mu^{\prime}}^{i}T^{il}_{\mu^{\prime}\mu\xi\xi^{\prime}}G^{lj}_{\xi\xi^{\prime},\nu\nu^{\prime}}(\omega)+\sum_{\alpha}\Big{[}\Big{(}V^{i}_{\mu^{\prime}\alpha}+\sum_{l}\sum_{\beta}D^{l}_{\beta}T^{il}_{\mu^{\prime}\alpha\beta\beta}\Big{)}G^{ij}_{\mu\alpha,\nu\nu^{\prime}}(\omega)$
(A.5)
$\displaystyle-\Big{(}V^{i}_{\alpha\mu}+\sum_{l}\sum_{\beta}D^{l}_{\beta}T^{il}_{\alpha\mu\beta\beta}\Big{)}G^{ij}_{\alpha\mu^{\prime},\nu\nu^{\prime}}(\omega)\Big{]}.$

If we take $V^{i}_{\mu\mu^{\prime}}=\delta_{\mu\mu^{\prime}}V^{i}_{\mu}$ and
$T^{ij}_{\mu\mu^{\prime}\nu\nu}=0$ for $\mu\neq\mu^{\prime}$ in Eq. (A.5) (MI-
SF transition and insulating states with one atom per site), the equation can
be simplified by setting
$E^{i}_{\mu}=V^{i}_{\mu}+z\sum_{\nu}D_{\nu}^{i}T_{\mu\mu\nu\nu}$.

Then to calculate the Green’s function matrix $\mathbf{G}(\mathbf{k},\omega)$
introduced in Sec. III, we combine Eq. (6) and Eq. (7) and have

$\displaystyle\mathbf{G}(\mathbf{k},\omega)-\epsilon(\mathbf{k})\mathbf{N}_{11}(\omega)\mathbf{G}(\mathbf{k},\omega)$
$\displaystyle=$
$\displaystyle\mathbf{N}_{11}(\omega)+\epsilon(\mathbf{k})\mathbf{N}_{12}(\omega)\mathbf{G}^{\prime}(\mathbf{k},\omega),$
$\displaystyle\mathbf{G}^{\prime}(\mathbf{k},\omega)-\epsilon(\mathbf{k})\mathbf{N}_{22}(\omega)\mathbf{G}^{\prime}(\mathbf{k},\omega)$
$\displaystyle=$
$\displaystyle\mathbf{N}_{21}(\omega)+\epsilon(\mathbf{k})\mathbf{N}_{21}(\omega)\mathbf{G}(\mathbf{k},\omega),$
(A.7)

where the matrix $\mathbf{G}^{\prime}(\mathbf{k},\omega)$ is defined as

$\displaystyle[\mathbf{G}^{\prime}]_{\alpha\beta}(\mathbf{k},\omega)=\sum_{\mu\mu^{\prime}\nu\nu^{\prime}}d^{\alpha}_{\mu\mu^{\prime}}d^{\beta}_{\nu\nu^{\prime}}G_{\mu\mu^{\prime},\nu\nu^{\prime}}(\mathbf{k},\omega).$
(A.8)

So finally we can solve $\mathbf{G}(\mathbf{k},\omega)$ as Eq. (8).

###  VII.2 determination of the symmetry properties in the MI-SF phase
transition

In this section, we will give detailed analytical calculations to determine
the symmetry properties in the MI-SF phase transition process. Obviously, we
need to consider ferromagnetic and antiferromagnetic cases separately.

####  VII.2.1 ferromagnetic case ($U_{2}<0$)

We find the coefficients satisfy the relationships

$\displaystyle\sum_{\alpha}\sum_{s}|c_{s}|^{2}\left[(A^{(\alpha)s,s+\alpha}_{n,n-1})^{2}+(A^{(\alpha)s,s+\alpha}_{n,n+1})^{2}\right]$
(B.1) $\displaystyle=$
$\displaystyle\sum_{\alpha}\langle\psi_{0}|\hat{a}_{\alpha}\hat{a}_{\alpha}^{\dagger}|\psi_{0}\rangle$
$\displaystyle=$ $\displaystyle n+3,$
$\displaystyle\sum_{\alpha}\sum_{s}|c_{s}|^{2}\left[(B^{(\alpha)s,s-\alpha}_{n,n-1})^{2}\right]$
(B.2) $\displaystyle=$
$\displaystyle\sum_{\alpha}\langle\psi_{0}|\hat{a}_{\alpha}^{\dagger}\hat{a}_{\alpha}|\psi_{0}\rangle$
$\displaystyle=$ $\displaystyle n,$

and

$\displaystyle\sum_{n^{\prime}=n\pm
1}\sum_{s}c_{s-\alpha_{1}}c_{s-\alpha_{2}}^{*}A^{(\alpha_{1})s-\alpha_{1},s}_{n,n^{\prime}}A^{(\alpha_{2})s-\alpha_{2},s}_{n,n^{\prime}}$
(B.3) $\displaystyle=$
$\displaystyle\langle\psi_{0}|\hat{a}_{\alpha_{2}}\hat{a}_{\alpha_{1}}^{\dagger}|\psi_{0}\rangle=\langle\psi_{0}|\hat{a}_{\alpha_{1}}^{\dagger}\hat{a}_{\alpha_{2}}|\psi_{0}\rangle~{}~{}\text{for
$\alpha_{1}\neq\alpha_{2}$}$ $\displaystyle=$
$\displaystyle\sum_{s}c_{s+\alpha_{2}}c_{s+\alpha_{1}}^{*}B^{(\alpha_{1})s+\alpha_{1},s}_{n,n-1}B^{(\alpha_{2})s+\alpha_{2},s}_{n,n-1}.$

Then the matrix $\mathbf{N}_{11}$ can be expanded as

$\displaystyle\mathbf{N}_{11}$ $\displaystyle=$
$\displaystyle(\Omega_{2}+\Omega_{3})\mathbf{M}_{1}^{\prime}+\sum_{s}(\Omega_{1}-\Omega_{2})\mathbf{M}_{2s}^{\prime}$
(B.4) $\displaystyle+\Omega_{2}\mathbf{I}_{3\times 3},$

where $\mathbf{M}_{1}^{\prime}$ and $\mathbf{M}_{2s}^{\prime}$ are defined as

$\displaystyle\mathbf{M}_{1}^{\prime}=\begin{pmatrix}\mathcal{M}_{0,0}&\mathcal{M}_{0,1}&\mathcal{M}_{0,1}\\\
\mathcal{M}_{1,0}&\mathcal{M}_{1,1}&\mathcal{M}_{1,-1}\\\
\mathcal{M}_{-1,0}&\mathcal{M}_{-1,1}&\mathcal{M}_{-1,-1}\end{pmatrix},$ (B.5)

with
$\mathcal{M}_{\alpha_{1},\alpha_{2}}=\sum_{s}c_{s+\alpha_{1}}c^{*}_{s+\alpha_{2}}B^{(\alpha_{1})s+\alpha_{1},s}_{n,n-1}B^{(\alpha_{2})s+\alpha_{2},s}_{n,n-1}$,
and

$\displaystyle\mathbf{M}_{2s}^{\prime}=\begin{pmatrix}\mathcal{M}^{s}_{0,0}&\mathcal{M}^{s}_{0,1}&\mathcal{M}^{s}_{0,-1}\\\
\mathcal{M}^{s}_{1,0}&\mathcal{M}^{s}_{1,1}&\mathcal{M}^{s}_{1,-1}\\\
\mathcal{M}^{s}_{-1,0}&\mathcal{M}^{s}_{-1,1}&\mathcal{M}^{s}_{-1,-1}\end{pmatrix},$
(B.6)

with
$\mathcal{M}^{s}_{\alpha_{1},\alpha_{2}}=c_{s-\alpha_{1}}c^{*}_{s-\alpha_{2}}A^{(\alpha_{1})s-\alpha_{1},s}_{n,n-1}A^{(\alpha_{2})s-\alpha_{2},s}_{n,n-1}$.

So $\lambda_{m}\mathbf{I}_{3\times 3}-\mathbf{N}_{11}$ can be expanded as the
summation of a series of hermite matrices

$\displaystyle\lambda_{m}\mathbf{I}_{3\times 3}-\mathbf{N}_{11}$
$\displaystyle=$ $\displaystyle(\Omega_{2}+\Omega_{3})(n\mathbf{I}_{3\times
3}-\mathbf{M}_{1}^{\prime})+\sum_{s}(\Omega_{2}-\Omega_{1})\mathbf{M}_{2s}^{\prime}.$

To prove the positive-semidefinite property of a matrix, we need to calculate
the $3$ upper left determinants. Before the calculation, we give some
inequality relations.

$\displaystyle\mathcal{M}_{\alpha_{1},\alpha_{1}}\mathcal{M}_{\alpha_{2},\alpha_{2}}-\mathcal{M}_{\alpha_{1},\alpha_{2}}\mathcal{M}_{\alpha_{2},\alpha_{1}}$
(B.8) $\displaystyle=$
$\displaystyle\sum_{s_{1}}|P_{B}^{(\alpha_{1})}(s_{1})|^{2}\sum_{s_{2}}|P_{B}^{(\alpha_{2})}(s_{2})|^{2}-\sum_{s_{1},s_{2}}P_{B}^{(\alpha_{1})}(s_{1})[P_{B}^{(\alpha_{2})}(s_{1})]^{*}[P_{B}^{(\alpha_{1})}(s_{2})]^{*}P_{B}^{(\alpha_{2})}(s_{2})$
$\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{s_{1},s_{2}}\left|P_{B}^{(\alpha_{1})}(s_{1})P_{B}^{(\alpha_{2})}(s_{2})-P_{B}^{(\alpha_{2})}(s_{1})P_{B}^{(\alpha_{1})}(s_{2})\right|^{2}$
$\displaystyle\geq$ $\displaystyle 0,$
$\displaystyle\mathcal{M}_{1,1}\mathcal{M}_{0,-1}\mathcal{M}_{-1,0}+\mathcal{M}_{-1,-1}\mathcal{M}_{0,1}\mathcal{M}_{1,0}-\mathcal{M}_{1,-1}\mathcal{M}_{0,1}\mathcal{M}_{-1,0}-\mathcal{M}_{-1,1}\mathcal{M}_{0,-1}\mathcal{M}_{1,0}$
(B.9) $\displaystyle=$
$\displaystyle\sum_{s}\left[|P^{(1)}_{B}(s)|^{2}|\mathcal{M}_{0,-1}|^{2}+|P^{(-1)}_{B}(s)|^{2}|\mathcal{M}_{0,1}|^{2}-P^{(1)}_{B}(s)[P^{(-1)}_{B}(s)]^{*}\mathcal{M}_{0,1}\mathcal{M}_{0,-1}^{*}-[P^{(1)}_{B}(s)]^{*}P^{(-1)}_{B}(s)\mathcal{M}_{0,-1}\mathcal{M}_{0,1}^{*}\right]$
$\displaystyle=$
$\displaystyle\sum_{s}\left|P^{(1)}_{B}(s)\mathcal{M}_{0,-1}^{*}-P^{(-1)}_{B}(s)\mathcal{M}_{0,1}^{*}\right|^{2}$
$\displaystyle\geq$ $\displaystyle 0,$

where $P_{B}^{(\alpha)}(s)$ is defined in Eq. (30).

For matrix $n\mathbf{I}_{3\times 3}-\mathbf{M}_{1}^{\prime}$, we have

$\displaystyle\text{UL}_{1}=n-\mathcal{M}_{0,0}=n-\langle\psi_{0}|\hat{a}_{0}^{\dagger}\hat{a}_{0}|\psi_{0}\rangle\geq
0,$ (B.10) $\displaystyle\text{UL}_{2}$ $\displaystyle=$
$\displaystyle(n-\mathcal{M}_{0,0})(n-\mathcal{M}_{1,1})-\mathcal{M}_{0,1}\mathcal{M}_{1,0}$
(B.11) $\displaystyle\geq$ $\displaystyle
n(n-\mathcal{M}_{0,0}-\mathcal{M}_{1,1})$ $\displaystyle\geq$ $\displaystyle
0,$ $\displaystyle\text{UL}_{3}$ $\displaystyle=$
$\displaystyle(n-\mathcal{M}_{0,0})[(n-\mathcal{M}_{1,1})(n-\mathcal{M}_{-1,-1})-\mathcal{M}_{1,-1}\mathcal{M}_{-1,1}]$
(B.12)
$\displaystyle+\mathcal{M}_{0,1}[-\mathcal{M}_{1,0}(n-\mathcal{M}_{-1,-1})-\mathcal{M}_{1,-1}\mathcal{M}_{-1,0}]$
$\displaystyle-\mathcal{M}_{0,-1}[\mathcal{M}_{1,0}\mathcal{M}_{-1,1}+\mathcal{M}_{-1,0}(n-\mathcal{M}_{1,1})]$
$\displaystyle\geq$ $\displaystyle
n(n-\mathcal{M}_{0,0})(n-\mathcal{M}_{1,1}-\mathcal{M}_{-1,-1})-n\mathcal{M}_{0,0}\mathcal{M}_{1,1}-n\mathcal{M}_{0,0}\mathcal{M}_{-1,-1}$
$\displaystyle=$ $\displaystyle 0.$

We have applied the relation
$\mathcal{M}_{0,0}+\mathcal{M}_{1,1}+\mathcal{M}_{-1,-1}=n$ in the calculation
of $\text{UL}_{l}$ ($l=1,2,3$). So matrix $n\mathbf{I}_{3\times
3}-\mathbf{M}_{1}^{\prime}$ is positive-semidefinite.

Also for matrix $\mathbf{M}_{2s}^{\prime}$, we calculate
$\text{UL}^{s}_{1}\geq 0$ and $\text{UL}^{s}_{2}=\text{UL}^{s}_{3}=0$. So
combined with the conditions $\Omega_{2}+\Omega_{3}>0$ and
$\Omega_{2}-\Omega_{1}>0$, we conclude from Eq. (VII.2.1) that
$\lambda_{m}\mathbf{I}_{3\times 3}-\mathbf{N}_{11}$ is positive-semidefinite,
which indicates that the eigenvalues of $\mathbf{N}_{11}$ satisfy the
condition $\lambda\leq\lambda_{m}$.

Then we need to consider the conditions for $\lambda=\lambda_{m}$. Based on
the Minkowski determinant theorem Marcus and Minc (1964), $\text{UL}_{3}$ in
Eq. (B.12) has to equal to zero, which indicates

$\displaystyle\frac{P^{(\alpha_{1})}_{B}(s)}{P^{(\alpha_{2})}_{B}(s)}=R_{\alpha_{1},\alpha_{2}},$
(B.13)

where $R_{\alpha_{1}\alpha_{2}}$ is an $s$-independent number. And we assume
that $P_{B}^{(\alpha_{2})}(s)\neq 0$.

Moreover, we note that the eigenvectors of matrices $n\mathbf{I}_{3\times
3}-\mathbf{M}_{1}^{\prime}$ and $\mathbf{M}_{2s}^{\prime}$ are the same, which
would induce

$\displaystyle 1+R_{1,0}r_{1,0}+R_{-1,0}r_{-1,0}=0,$ (B.14)

where we have also defined

$\displaystyle\frac{P^{(\alpha_{1})}_{A}(s)}{P^{(\alpha_{2})}_{A}(s)}=r_{\alpha_{1},\alpha_{2}}.$
(B.15)

Also we note that $R_{1,0}=-r_{-1,0}$ because

$\displaystyle\frac{A^{(-1)s+1,s}_{n,n-1}}{A^{(0)s,s}_{n,n-1}}=\frac{\langle
n-1,s;n-1|\hat{\Theta}\hat{a}_{-1}^{\dagger}|n,s+1;n\rangle}{\langle
n-1,s;n-1|\hat{\Theta}\hat{a}_{0}^{\dagger}|n,s;n\rangle}$ (B.16)
$\displaystyle=$ $\displaystyle\frac{\langle
n-1,s;n-1|\hat{a}_{-1}^{\dagger}\hat{\Theta}-2\hat{a}_{1}|n,s+1;n\rangle}{\langle
n-1,s;n-1|\hat{a}_{0}^{\dagger}\hat{\Theta}+2\hat{a}_{0}|n,s;n\rangle}$
$\displaystyle=$ $\displaystyle-\frac{\langle
n-1,s;n-1|\hat{a}_{1}|n,s+1;n\rangle}{\langle
n-1,s;n-1|\hat{a}_{0}|n,s;n\rangle}$ $\displaystyle=$
$\displaystyle-\frac{B^{(1)s+1,s}_{n,n-1}}{B^{(0)s,s}_{n,n-1}},$

where operator $\hat{\Theta}=\hat{a}_{0}^{2}-2\hat{a}_{1}\hat{a}_{-1}$
describes the destruction of a spin singlet pair.

Similarly, we have $R_{-1,0}=-r_{1,0}$. So Eq. (B.14) is equivalent to

$\displaystyle R_{1,0}R_{-1,0}=\frac{1}{2}.$ (B.17)

Then we can calculate $\langle\vec{S}\rangle^{2}$ in the system

$\displaystyle\langle\vec{S}\rangle^{2}=\langle\hat{S}^{+}\rangle\langle\hat{S}^{-}\rangle+\langle\hat{S}_{z}\rangle^{2}$
(B.18) $\displaystyle=$ $\displaystyle
2|R_{1,0}+R_{-1,0}^{*}|^{2}\Big{(}\sum_{s}|P^{(0)}_{B}(s)|^{2}\Big{)}^{2}+$
$\displaystyle(|R_{1,0}|^{2}-|R_{-1,0}|^{2})^{2}\Big{(}\sum_{s}|P^{(0)}_{B}(s)|^{2}\Big{)}^{2}$
$\displaystyle=$ $\displaystyle n^{2}.$

The corresponding eigenvector involved in the calculation, which can be
expressed as $(1,R_{1,0},R_{-1,0})$, determines the symmetry properties of the
superfluid phase. So we conclude from Eq. (B.17) that the superfluid phase in
the excitation is in ferromagnetic state.

####  VII.2.2 antiferromagnetic case ($U_{2}>0$)

For odd filling number case, the matrix $\mathbf{N}_{11}(0)$ can be expressed
as

$\displaystyle\mathbf{N}_{11}(0)=3(\Upsilon_{2}+\Upsilon_{4})\mathbf{I}_{3\times
3}+K_{1}\mathbf{M}_{1}+K_{2}\mathbf{M}_{2},$ (B.19)

where matrices $\mathbf{M}_{1}$ and $\mathbf{M}_{2}$ are

$\displaystyle\mathbf{M}_{1}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}|c_{0}|^{2}&c_{0}c_{1}^{*}&c_{0}c_{-1}^{*}\\\
c_{0}^{*}c_{1}&|c_{1}|^{2}&c_{1}c_{-1}^{*}\\\
c_{0}^{*}c_{-1}&c_{-1}c_{1}^{*}&|c_{-1}|^{2}\end{pmatrix},$ (B.20)

and

$\displaystyle\mathbf{M}_{2}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}|c_{0}|^{2}&-c_{0}^{*}c_{-1}&-c_{0}^{*}c_{1}\\\
-c_{0}c_{-1}^{*}&|c_{-1}|^{2}&c_{1}c_{-1}^{*}\\\
-c_{0}c_{1}^{*}&c_{-1}c_{1}^{*}&|c_{1}|^{2}\end{pmatrix}.$ (B.21)

###  VII.3 analysis of the spin phase boundary in MI phase with one atom per
site

In this section, we give detailed analysis of the minimum value of Eq. (59),
which would induce different phases in the MI phase with one atom per site.

For the excitation from $|1,0;1\rangle$ to $|1,1;1\rangle$, we start from Eq.
(59) and define

$\displaystyle f(\eta)$ $\displaystyle=$
$\displaystyle(q+zJ_{2}+J_{1}\eta)^{2}-(J_{2}-J_{1})^{2}\eta^{2},$ (C.1)

with $-z\leq\eta\leq z$. Then we calculate the minimum value of $f(\eta)$.

When $J_{2}\geq 2J_{1}$, the minimum is obtained by comparing $f(z)$ and
$f(-z)$, and we have

$\displaystyle
f(\eta)\Big{|}_{\text{min}}=(q+zJ_{2}-zJ_{1})^{2}-z^{2}(J_{2}-J_{1})^{2}.$
(C.2)

When $J_{2}<2J_{1}$, the case becomes more complicated, we need to first
decide whether the global minimum point is within the range $[-z,z]$.

If $q\geq zJ_{2}-z\frac{J_{2}^{2}}{J_{1}}$, the minimum is on the point
$\eta=-z$ with

$\displaystyle
f(\eta)\Big{|}_{\text{min}}=(q+zJ_{2}-zJ_{1})^{2}-z^{2}(J_{2}-J_{1})^{2},$
(C.3)

and if $q<zJ_{2}-z\frac{J_{2}^{2}}{J_{1}}$, the minimum is on the global
minimum point with

$\displaystyle
f(\eta)\Big{|}_{\text{min}}=-\frac{(J_{2}-J_{1})^{2}}{J_{1}^{2}-(J_{2}-J_{1})^{2}}(q+zJ_{2})^{2}<0.$
(C.4)

###  VII.4 SBO method for MI phase with two atoms per site

In this section, we give the self-consistent equations for the spin-excitation
in the MI phase with two atoms per site. Then by solving the equations, we
will derive self-consistent results of the phase transition process.

For MI phase with two atoms per site,
$H_{\alpha\alpha^{\prime},\beta\beta^{\prime}}$ in Eq. (67) can be expressed
as

$\displaystyle H_{ss,ss}$ $\displaystyle=$
$\displaystyle-\frac{20t^{2}}{3U_{0}},$ (D.1) $\displaystyle
H_{ms,sm^{\prime}}$ $\displaystyle=$
$\displaystyle-\frac{8t^{2}}{3U_{0}}\delta_{mm^{\prime}},$ (D.2)
$\displaystyle H_{mm^{\prime},ss}$ $\displaystyle=$
$\displaystyle\Big{(}-\frac{20t^{2}}{3U_{0}}+\frac{3U_{2}}{z}\Big{)}\delta_{mm^{\prime}},$
(D.3) $\displaystyle H_{ms,m^{\prime}s}$ $\displaystyle=$
$\displaystyle(-1)^{-m+1}\frac{8t^{2}}{3U_{0}}\delta_{m,-m^{\prime}},$ (D.4)
$\displaystyle H_{sm,lm^{\prime}}$ $\displaystyle=$
$\displaystyle\frac{4\sqrt{7}t^{2}}{3U_{0}}C_{0,0;2,l}^{2,l}C_{2,m;2,m^{\prime}}^{2,l}\delta_{l,m+m^{\prime}},$
(D.5) $\displaystyle H_{ml,m^{\prime}l^{\prime}}$ (D.6) $\displaystyle=$
$\displaystyle\delta_{m+m^{\prime},l+l^{\prime}}\Big{[}\left(\frac{6U_{2}}{z}-\frac{16t^{2}}{3U_{0}}\right)C_{2,m;2,m^{\prime}}^{0,m+m^{\prime}}C_{2,l;2,l^{\prime}}^{0,l+l^{\prime}}$
$\displaystyle+\left(\frac{6U_{2}}{z}-\frac{4t^{2}}{U_{0}}\right)C_{2,m;2,m^{\prime}}^{1,m+m^{\prime}}C_{2,l;2,l^{\prime}}^{1,l+l^{\prime}}$
$\displaystyle+\left(\frac{6U_{2}}{z}-\frac{8t^{2}}{3U_{0}}\right)C_{2,m;2,m^{\prime}}^{2,m+m^{\prime}}C_{2,l;2,l^{\prime}}^{2,l+l^{\prime}}$
$\displaystyle+\left(\frac{6U_{2}}{z}-\frac{4t^{2}}{U_{0}}\right)C_{2,m;2,m^{\prime}}^{3,m+m^{\prime}}C_{2,l;2,l^{\prime}}^{3,l+l^{\prime}}$
$\displaystyle+\left(\frac{6U_{2}}{z}-\frac{12t^{2}}{U_{0}}\right)C_{2,m;2,m^{\prime}}^{4,m+m^{\prime}}C_{2,l;2,l^{\prime}}^{4,l+l^{\prime}}\Big{]}.$

In the above expressions, we have assumed that
$|m\rangle,|l\rangle,|m^{\prime}\rangle,|l^{\prime}\rangle$ represent states
that belong to set $\\{S=2\\}$. And $C_{S_{1},m_{1};S_{2},m_{2}}^{S,m}$ are
the known Clebsch-Gordon coefficients.

For singlet ($|s\rangle$) to ferromagnetic ($|2\rangle$) excitation, we give
the equations connecting to $G_{s2,2s}(\mathbf{k},\omega)$

$\displaystyle\mathbf{M}_{s\rightarrow
2}\begin{pmatrix}G_{s2,2s}(\mathbf{k},\omega)\\\
G_{-2s,2s}(\mathbf{k},\omega)\\\
G_{02,2s}(\mathbf{k},\omega)\end{pmatrix}=\begin{pmatrix}D_{s2}\\\ 0\\\
0\end{pmatrix},$ (D.7)

with

$\displaystyle\mathbf{M}_{s\rightarrow 2}=\omega\mathbf{I}_{3\times
3}+\begin{pmatrix}E_{s}-E_{2}&0&\gamma\\\ 0&E_{-2}-E_{s}&0\\\
\gamma&0&E_{0}-E_{2}\end{pmatrix}+$
$\displaystyle\zeta(\mathbf{k})\begin{pmatrix}D_{s2}H_{s2,s2}&D_{s2}H_{2s,-2s}&D_{s2}H_{2s,02}\\\
D_{-2s}H_{s-2,s2}&D_{-2s}H_{s-2,-2s}&D_{-2s}H_{s-2,02}\\\
D_{02}H_{20,s2}&D_{02}H_{20,-2s}&D_{02}H_{20,02}\end{pmatrix},$

where $E_{\mu}=V_{\mu}+z\sum_{\nu}D_{\nu}H_{\mu\mu,\nu\nu}$ and
$\gamma=z\sum_{\nu}D_{\nu}H_{0s,\nu\nu}$. We have eliminated
$G_{-20,2s}(\mathbf{k},\omega)$ and $G_{-11,2s}(\mathbf{k},\omega)$ in the
equations because $E_{2}<E_{m}$ ($m=-2,\pm 1,0$) and we have $D_{\mu}=0$ for
$\mu\neq s,2$ in the calculation. Then we can solve $D_{2}$ with integrations
of the Green’s functions.

## References

  * Bloch et al. (2008) I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
  * Jaksch et al. (1998) D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, Phys. Rev. Lett. 81, 3108 (1998).
  * Greiner et al. (2002) M. Greiner, O. Mandel, T. Esslinger, T. W. Hänsch, and I. Bloch, Nature 415, 39 (2002).
  * Jördens et al. (2010) R. Jördens, L. Tarruell, D. Greif, T. Uehlinger, N. Strohmaier, H. Moritz, T. Esslinger, L. De Leo, C. Kollath, A. Georges, et al., Phys. Rev. Lett. 104, 180401 (2010).
  * Medley et al. (2011) P. Medley, D. M. Weld, H. Miyake, D. E. Pritchard, and W. Ketterle, Phys. Rev. Lett. 106, 195301 (2011).
  * Jiang et al. (2016) J. Jiang, L. Zhao, S.-T. Wang, Z. Chen, T. Tang, L.-M. Duan, and Y. Liu, Phys. Rev. A 93, 063607 (2016).
  * Stamper-Kurn and Ueda (2013) D. M. Stamper-Kurn and M. Ueda, Rev. Mod. Phys. 85, 1191 (2013).
  * Hall et al. (1998) D. S. Hall, M. R. Matthews, J. R. Ensher, C. E. Wieman, and E. A. Cornell, Phys. Rev. Lett. 81, 4531 (1998).
  * Matthews et al. (1998) M. R. Matthews, D. S. Hall, D. S. Jin, J. R. Ensher, C. E. Wieman, E. A. Cornell, F. Dalfovo, C. Minniti, and S. Stringari, Phys. Rev. Lett. 81, 243 (1998).
  * Widera et al. (2005) A. Widera, F. Gerbier, S. Fölling, T. Gericke, O. Mandel, and I. Bloch, Phys. Rev. Lett. 95, 190405 (2005).
  * Widera et al. (2006) A. Widera, F. Gerbier, S. Fölling, T. Gericke, O. Mandel, and I. Bloch, New Journal of Physics 8, 152 (2006).
  * Ho (1998) T.-L. Ho, Phys. Rev. Lett. 81, 742 (1998).
  * Ohmi and Machida (1998) T. Ohmi and K. Machida, Journal of the Physical Society of Japan 67, 1822 (1998).
  * Katsura and Tasaki (2013) H. Katsura and H. Tasaki, Phys. Rev. Lett. 110, 130405 (2013).
  * Snoek and Zhou (2004) M. Snoek and F. Zhou, Phys. Rev. B 69, 094410 (2004).
  * Demler and Zhou (2002) E. Demler and F. Zhou, Phys. Rev. Lett. 88, 163001 (2002).
  * Imambekov et al. (2003) A. Imambekov, M. Lukin, and E. Demler, Phys. Rev. A 68, 063602 (2003).
  * Tsuchiya et al. (2004) S. Tsuchiya, S. Kurihara, and T. Kimura, Phys. Rev. A 70, 043628 (2004).
  * Kimura et al. (2005) T. Kimura, S. Tsuchiya, and S. Kurihara, Phys. Rev. Lett. 94, 110403 (2005).
  * Pai et al. (2008) R. V. Pai, K. Sheshadri, and R. Pandit, Phys. Rev. B 77, 014503 (2008).
  * Zhou and Semenoff (2006) F. Zhou and G. W. Semenoff, Phys. Rev. Lett. 97, 180411 (2006).
  * Bernier et al. (2007) J.-S. Bernier, K. Sengupta, and Y. B. Kim, Phys. Rev. B 76, 014502 (2007).
  * Song and Zhou (2008) J. L. Song and F. Zhou, Phys. Rev. A 77, 033628 (2008).
  * Natu et al. (2015) S. S. Natu, J. H. Pixley, and S. Das Sarma, Phys. Rev. A 91, 043620 (2015).
  * So and Ueda (2016) K. H. Z. So and M. Ueda, arXiv:1604.01163 (2016).
  * de Forges de Parny et al. (2013) L. de Forges de Parny, F. Hébert, V. G. Rousseau, and G. G. Batrouni, Phys. Rev. B 88, 104509 (2013).
  * Batrouni et al. (2009) G. G. Batrouni, V. G. Rousseau, and R. T. Scalettar, Phys. Rev. Lett. 102, 140402 (2009).
  * Imambekov et al. (2004) A. Imambekov, M. Lukin, and E. Demler, Phys. Rev. Lett. 93, 120405 (2004).
  * Rizzi et al. (2005) M. Rizzi, D. Rossini, G. De Chiara, S. Montangero, and R. Fazio, Phys. Rev. Lett. 95, 240404 (2005).
  * Apaja and Syljuåsen (2006) V. Apaja and O. F. Syljuåsen, Phys. Rev. A 74, 035601 (2006).
  * Svidzinsky and Chui (2003) A. A. Svidzinsky and S. T. Chui, Phys. Rev. A 68, 043612 (2003).
  * Chung and Yip (2009) M.-C. Chung and S. Yip, Phys. Rev. A 80, 053615 (2009).
  * Rodríguez et al. (2011) K. Rodríguez, A. Argüelles, A. K. Kolezhuk, L. Santos, and T. Vekua, Phys. Rev. Lett. 106, 105302 (2011).
  * Haley and Erdös (1972) S. B. Haley and P. Erdös, Phys. Rev. B 5, 1106 (1972).
  * Sajna et al. (2015) A. S. Sajna, T. P. Polak, R. Micnas, and P. Rożek, Phys. Rev. A 92, 013602 (2015).
  * Law et al. (1998) C. K. Law, H. Pu, and N. P. Bigelow, Phys. Rev. Lett. 81, 5257 (1998).
  * Wu (1996) Y. Wu, Phys. Rev. A 54, 4534 (1996).
  * Sheshadri et al. (1993) K. Sheshadri, H. Krishnamurthy, R. Pandit, and T. Ramakrishnan, EPL (Europhysics Letters) 22, 257 (1993).
  * Kimura (2013) T. Kimura, Phys. Rev. A 87, 043624 (2013).
  * Ohliger (2008) M. Ohliger, Diploma Thesis (2008).
  * Chen and Levy (1973) H. H. Chen and P. M. Levy, Phys. Rev. B 7, 4267 (1973).
  * Yip (2003) S. K. Yip, Phys. Rev. Lett. 90, 250402 (2003).
  * Stenger et al. (1999) J. Stenger, S. Inouye, D. M. Stamper-Kurn, H.-J. Miesner, A. P. Chikkatur, and W. Ketterle, Nature 396, 345 (1999).
  * Micnas (1976) R. Micnas, Journal of Physics C Solid State Physics 9, 3307 (1976).
  * Papanicolaou (1988) N. Papanicolaou, Nuclear Physics B 305, 367 (1988).
  * Marcus and Minc (1964) M. Marcus and H. Minc, _A survey of matrix theory and matrix inequalities_ (Allyn and Bacon,, 1964).