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

TL;DR

This paper analyzes the phase transitions and spin excitations of spin-1 bosons in optical lattices, providing analytical and numerical insights into the MI-SF transition and various magnetic phases, with results aligning with Monte-Carlo simulations.

Contribution

It introduces an extended SBO method to study spin-1 bosons, deriving phase diagrams and spin phase boundaries considering quantum and thermal fluctuations.

Findings

Quantitative phase diagram of MI-SF transition agrees with Monte-Carlo results.

Finite-temperature effects are underestimated by mean-field approximation.

Boundaries of various spin phases in Mott lobes are analytically and numerically determined.

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.

Equations

$$\hat{H}$$

$${\vec{S}}_i$$

$$E^{(0)}(S_i,n_i)$$

$$\hat{H}$$

$$G^{i\alpha j\beta }(t-t^{\prime })$$

$$(\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 }$$

$$G^{\alpha \beta }(\mathbf{k},\omega )$$

$$\mathbf{G}(\mathbf{k},\omega )$$

$$\mathbf{\Pi }(\mathbf{k},\omega )$$

$$\mathbf{\Gamma }(\mathbf{k},\omega )$$

$$\left(\begin{array}{cc}{\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{array}\right)$$

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

$$A_{n,n\pm 1}^{(\alpha )s,s+\alpha }$$

$$B_{n,n-1}^{(\alpha )s,s-\alpha }$$

$${\Omega }_1$$

$${\Omega }_2$$

$${\Omega }_3$$

$${\mathbf{N}}_{11}^{0,0}$$

$${\mathbf{N}}_{11}^{1,1}$$

$${\mathbf{N}}_{11}^{-1,-1}$$

$$\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{]},$$

$${\mathbf{N}}_{11}^{1,0}$$

$${\mathbf{N}}_{11}^{-1,0}$$

$${\mathbf{N}}_{11}^{-1,1}$$

$$\frac{P_B^{({\alpha }_1)}(s)}{P_B^{({\alpha }_2)}(s)}$$

$$R_{1,0}{R}_{-1,0}=\frac{1}{2},$$

$$P_A^{(\alpha )}(s)$$

$$P_B^{(\alpha )}(s)$$

$$z{t}_c[(n+1){\Omega }_2+n{\Omega }_3]=1.$$

$${\lambda }_0=3({{\rm Y}}_2+{{\rm Y}}_4),$$

$${\lambda }_{\pm }=\frac{1}{2}({{\rm Y}}_1+{{\rm Y}}_3+7{{\rm Y}}_2+7{{\rm Y}}_4)$$

Full text

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

[TABLE]

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

[TABLE]

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

[TABLE]

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)

[TABLE]

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$

[TABLE]

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)

[TABLE]

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)

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

and

[TABLE]

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

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

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

[TABLE]

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

where

[TABLE]

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

[TABLE]

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

[TABLE]

where

[TABLE]

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)

[TABLE]

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

[TABLE]

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.

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).

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.

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)

[TABLE]

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,

[TABLE]

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

[TABLE]

and

[TABLE]

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

[TABLE]

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)$

[TABLE]

with

[TABLE]

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

[TABLE]

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

[TABLE]

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).

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 [math] 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

[TABLE]

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

[TABLE]

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

[TABLE]

with

[TABLE]

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

[TABLE]

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.

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

[TABLE]

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

[TABLE]

with

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

then combined with the equation

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

[TABLE]

and

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.

[TABLE]

[TABLE]

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

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

[TABLE]

[TABLE]

[TABLE]

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

[TABLE]

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

[TABLE]

where we have also defined

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

and

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

[TABLE]

[TABLE]

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)$

[TABLE]

with

[TABLE]

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.