Quantum phase diagrams and time-of-flight pictures of spin-1 Bose systems in honeycomb optical lattices

By treating the hopping parameter as a perturbation, with the help of cumulant expansion and the re-summing technique, the one-particle Green’s function of a spin-1 Bose system in a honeycomb optical lattice is calculated analytically. By the use of the re-summed Green’s function, the quantum phase diagrams of the system in ferromagnetic cases as well as in antiferromagnetic cases are determined. It is found that in antiferromagnetic cases the Mott insulating states with even filling factor are more robust against the hopping parameter than that with odd filling factor, in agreement with results via other different approaches. Moreover, in order to illustrate the effectiveness of the re-summed Green’s function method in calculating time-of-flight pictures, the momentum distribution function of a honeycomb lattice spin-1 Bose system in the antiferromagnetic case is also calculated analytically and the corresponding time-of-flight absorption pictures are plotted.


Introduction
Since it was proposed theoretically [1,2] and realized exper imentally [3], spinor Bose-Einstein condensate has opened a new avenue to investigate magnetic systems by utilizing ultra cold atomic gases [4]. Immediately after the realization of optical lattice [5], confining spinor bosons in optical lattices becomes one of the hottest research topics in the field of cold atoms and has been extensively investigated [6][7][8][9][10], although the experimental realization of the optical lattice spinor Bose system has been achieved only very recently [11].
As is known, for atoms trapped in a magnetooptical trap, the spin degree of freedom is frozen and the atoms become effectively spinless. However, the quantum gases will keep their spin degree of freedom if they are trapped by purely optical means. Such an additional spin degrees of freedom, together with other tunable parameters, makes the spinor Bose lattice systems unique for studying quantum magnetic prop erties of matter (for review, see [12][13][14], even for quantum information processing [15]. Among spinor Bose systems in optical lattices, spin1 lattice Bose gases should be the simplest one and have attracted lots of attention. Such systems can be depicted by the extended Bose-Hubbard Hamiltonians [1,2] with spindependent terms which will dramatically change the properties of the system [12]. In fact, the nature of onsite spin-spin interaction can be either ferromagnetic ( 87 Rb) [16] or antiferromagnetic ( 23 Na) [3], depending on the relative magnitudes of the swave scattering lengths in spin0 and spin2 channels [17]. The possibility of appearance of dif ferent quantum phases in such systems has been predicted both analytically and numerically [6,7,[18][19][20][21][22]. To deter mine the quantum phase boundaries, a variety of methods have also been applied, such as meanfield methods [23][24][25], strong coupling expansion [27], and Monte Carlo simulations [26,28] et al. However, most of these efforts have been only devoted to simple Bravais lattice systems.
With the progress of experimental techniques, scalar Bose gases trapped in a nonBravais honeycomb optical lattice has been accomplished [29]. Due to its complexity, trapping a spinor Bose gas in the honeycomb lattice has been only real ized most recently [11]. However, theoretical studies on the quantum phase transitions of spinor Bose systems in hon eycomb lattice are less explored. In this paper, by utilizing a generalized Green's function method [30], together with cumulants expansion [31] and resummed Green's function technique [32,33], we will investigate the quantum phase transitions of the spin1 Bose system in a honeycomb optical lattice. As a demonstration of the usefulness of the general ized Green's function method, the timeofflight absorption pictures of the system in antiferromagnetic cases will also be presented. Our analytical results may serve as a reference object for the upcoming experiments.
The paper is arranged as follows. In section 2, the spinor Bose-Hubbard model is introduced, and the possible ground states of the pure local part of the spinor Bose-Hubbard Hamiltonian are discussed. In section 3, by the use of the cumulants expansion as well as the resummed Green's func tion technique, the oneparticle Green's function of a spinor Bose-Hubbard system in a honeycomb lattice is constructed. The phase boundary equations for both antiferromagnetic and ferromagnetic cases are presented in section 4, where the phase diagrams for both cases are also plotted. To cross check the validity of our analytical method, we also calculate the phase boundaries of spin1 Bose-Hubbard systems in triangu lar and square lattices and compare them with previous results [28,47], an excellent agreement is found. In section 5, as an example, the momentum distribution function of antiferro magnetic case with the filling factor n = 3 is calculated ana lytically, the corresponding timeofflight pictures are plotted. Section 6 is devoted to summary.

The model
A system of spin1 bosons trapped in a honeycomb optical lat tice can be portraited by the Hamiltonian [1,6,7,34] ˆˆ( ) ( )ˆ( ) where ˆ( ) ψ α r and ˆ( ) † ψ α r are the field operators for a spin1 boson in hyperfine states M is the mass of atom, μ is the chemical potential, ( ) V r represents the optical lattice potential and F is the vector of the three spin1 generators of the rotation group The coefficients c 0 and c 2 are defined as with a 0 and a 2 being the swave scattering lengths in singlet and quintet spin channels respectively. To investigate the ground state properties of the spinor Bose system as well as its quantum phase transitions, a convenient way is to artifi cially introduce a spin symmetry breaking term from the start and set it back to zero at the end of calculations. Thus, an additional Zeeman term with infinitesimal external magnetic field η has to be included in the Hamiltonian (1).
Due to the modulation of the optical lattice and the independ ence of the lattice potential on the spin of bosons, at extremely low temperature, the field operators ˆ( ) † ψ α r and ˆ( ) ψ α r can be expanded on the basis of spinindependent, orthonor mal Wannier functions ( ) − w r r i of the lowest band [35,36] aŝ a w a w r r r r r r , , where ˆ † α a i (ˆα a i ) is the bosonic creation (annihilation) opera tor on site i with α = m F . In the tightbinding limit, the Hamiltonian of the system (1) reduces to the extended Bose-Hubbard model for spin1 lattice bosons [6,7,23,34,37] where the nearestneighboring hopping parameter is and the interactions U 0, 2 are Note that ⟨ ⟩ i j , in ( As we can see from equation (5), the Hamiltonian consists of two parts, the hopping termˆ〈 .

This leads to
for odd filling factor n.
As for ferromagnetic interaction with U 2 < 0, to minimize the energy, the ground state of the system has to be the state of ⟩ |n n n , ; . Hence the constraint on chemical potential for different filling factor becomes With the configurations of the localized states in hand, we can now tune on the hopping parameter J gradually. It can be imagined that the localized states will be softened by quantum fluctuations and be transformed into some delocalized states while J goes beyond some critical value, and thus a quantum phase transition is expected.

Cumulants expansion of the Green's function
As is well known, for secondorder phase transition, the one particle Green's function of the system diverges at the critical point [39,40]. Meanwhile, the Green's function can reveal the momentum distribution of the system, which has a direct connection with the timeofflight measurements [13,[41][42][43] in real experiments. Based on these facts, in order to investi gate the quantum phase transition of the spinor boson opti cal lattice systems as well as the corresponding timeofflight absorption pictures, we are going to calculate the oneparticle Green's function of the system, which is defined as , , .
j j 1 (14) Here τ is imaginary time, ⟨ ⟩ ⋅ denotes the statistical average with respect to the full Hamiltonian (5), and ˆτ T is the imagi nary time ordering operator.
Since the hopping part (8) and the diagonal one (9) does not commute with each other, the exact eigenstates of the full Hamiltonian can hardly be determined. Thus we calculate the oneparticle Green's function (14) by treating the hopping term as perturbation in Dirac representation, i.e.
Tr e ,0 , and the time evolution (setting = 1) [44]. Using the expansion of ˆ( ) β U , 0 , it is easy to find that the perturbative expansion of the Green's function con sists of terms like with ⟨ ⟩ ⋅ 0 denoting the average with respect to Ĥ 0 . We set J ij = J if i and j are the nearest neighbors of each other, oth erwise J ij = 0. Now, the problem of calculating the oneparticle Green's function with respect to the full ground state is reduced to the calculation of the nparticle Green function with respect to the unperturbed state which is defined as It turns out that the calculation of ( ) G n 0 can be further simplified by means of the socalled cumulant expansions [31]. Since the eigenstates of the unperturbed part Ĥ 0 is local, each annihila tion (creation) operator in the expression of ( ) G n 0 has to be paired with a creation (annihilation) operator on the same site with the same quantum numbers. Therefore, ( ) G n 0 can be decomposed in in which all the operators are at the same lattice site [31,32,45]. In this case the decomposition of ( ) G n 0 is straightforward.
As an example, we decompose the oneparticle and two particle Green's functions in terms of the cumulants as and In order to simplify the discussion, the mth order cumu lant can be represented as an arrowline diagram composed of m incoming arrow lines (representing the creation opera tors) and m outgoing arrow lines (representing the annihi lation operators) at the same lattice site. For example, the first order and the second order cumulants can be repre Correspondingly, the hopping parameter J ij between nearest neighbor sites is also sketched as = J ij . Thus, with the help of cumulants expansion, can be easily expressed diagrammatically [31]. When the cancel lation effect of the numerator and the denominator in equa tion (15) is taken into account, it is found that the oneparticle of secondorder phase transition. This property can not be ful filled by just counting finite orders of perturbation parameter J/U 0 . On the other hand, J/U 0 is no longer a small quantity near the phase boundaries. Hence, we regroup the diagrams based on the number of loops rather than on the order of J. This is the socalled resummed Green's function technique [32,33]. At present we only take line diagrams into account, as the diagram containing n loops will be smaller than the corresponding line diagram at least by factor of 1/d n with d being the dimension of the system. This choice is going to be accurate in the limit → ∞ d . In this sense, our choice of the resummed Green's function may be looked upon as a type of meanfield approximation. Nevertheless, from the above discussion, we argue that in principle the resummed Green's function method can go beyond mean field by including loop diagrams constituted by higher order cumulants and related hoppings.
Since the honeycomb lattice is a nonBravais one, we need to divide the honeycomb lattice into two sublattices, denoted by sublattice A and sublattice B respectively [30]. In this case, the hopping amplitude from sublattice A to sublattice B J AB and the hopping from sublattice B to sublattice A J BA can be expressed in momentum space as where a is the lattice constant. The cumulants on different sublattices in terms of Matsubara frequency are denoted dia Here the conservation of angular momentum has already been taken into account, i.e. α α = ′ and ω ω = ′ m m has been adopted. Thus, the resummed oneparticle Green's function of the system is further expressed diagrammatically as only consists of con nected diagrams [31,32]. From the above diagrammatic rep resentation, it is not difficult to recognize that diagrams which only consists of the first order cumulants are line diagrams, while the higher order cumulants together with hoppings will form loops in corresponding diagrams.
In principle, the oneparticle Green's function (15) may be accurately calculated out by counting all the possible dia grams. However, in practice one can only calculate the Green's function perturbatively by selecting a specific group of dia grams. Actually, this selection is very subtle. On one hand, the Green's function should diverge at the phase boundaries (22) which corresponds to the following analytical expression Note that since the cumulants are local quantities on each lat tice site and there is no bias on these two sublattices in our system, to this end the cumulants on different sublattices will make the same contribution, i.e.  , ; . To do this, we first check how the creation and annihilation operators act on the eigenstate ⟩ |S m n , ; of Ĥ 0 . Since the crea tion operator ˆ † α a (ˆα a ) creates a spin1 particle with its spin ori entation specified by α, we can easily have where α M Smn , α N Smn , α O Smn and α P Smn are superposition para meters, and their detail expressions are presented in the appen dix. Combining equations (17), (18), (24) and (25), we obtain the following explicit expression of the first order cumulant in Matsubara space At zerotemperature limit, the above expression reduces to

Phase boundaries of spin-1 Bose systems on honeycomb optical lattices
For ferromagnetic interaction, both theoretical analysis [21] and quantum Mont Carlo results [28] indicate that the super fluidMott insulator phase transitions are always of second order. However, for antiferromagnetic interaction, the situa tion becomes complicated. The mean field analysis [21,46] as well as some numerical results [22,28] [22]. In the present work, we will concentrate on the cases with strong U 2 so that a sec ondorder phase transition can occur.
By summarizing the geometric sequences in equation (23) The divergency of the Green function (28) at the phase trans ition point [39] requires that Solving the above quadratic equation (29) and plugging in the culmulant result (27), finally we obtain the following phase boundary equation  [39,40]. In the fol lowing, we are going to discuss the phase boundaries case by case.
Let us first consider an antiferromagnetic system with even filling factor. In this case the unperturbed ground state is ⟩ | n 00 , i.e. S = 0 and m = 0 in the above equation (30). With the help of the recursion relations of the superposition parameters, it is easy to find that = α N 0 with ˜/ = U U U 2 2 0 and ˜/ µ µ = U 0 . The situation for odd filling cases is a little bit complicated, where the unperturbed ground state is threefold degenerated when the external field is absent. Hence, after taking the limit → η 0 in equation (30), we have to evaluate the following equation for all possible α and m We plot the phase boundary of the spin1 Bose system in a honeycomb lattice in figure 1 . When crossing the phase boundary, the system will go from the MI state to a superfluid phase with m = 0, the socalled polar state [34]. As shown in figure 1, the Mott insulator phase with even filling factor is considerably more stable against the super fluid phase than that with odd filling factor. Similar behav ior has also been found in square optical lattice systems via Monte Carlo simulation [28] and in triangular lattice systems via strong coupling expansion [47]. Actually, this remarkable feature is understandable as in the antiferromagnetic case with even filling factor, all particles in MI states form spinsinglet pairs, which presents an additional spin gap [6,34].
Considering that there is no existing results for the phase diagram of a spinor Bose system in a honeycomb lattice, in order to make cross check for our analytical result, we apply our method to triangular lattice system. We compare our ana lytical result with the result taken from [47] in figure 1 (right) and a good agreement is found.
For ferromagnetic cases, the phase boundaries can be obtained with the same procedure. From the recursion relation of the superposition para meters, we know that . By looking for the minimum of (30) with respect to the spin index α, the analytical expression of the phase boundary for ferromagnetic case of a spinor Bose system in a honeycomb lattice can be expressed as The phase diagram of a ferromagnetic system in a honey comb lattice with / = − U U 0.1 2 0 is plotted in figure 2 (left). The phase boundary separates the Mott insulator states ⟩ |n n , from superfluid phase. Our analytical result for ferromagnetic cases is also cross checked by being applied to square lattice system. The corresponding result is plotted in figure 2 (right) together with numerical result obtained via the diagonalized mean field approach [28], and a good agreement is also found.

Time-of-flight pictures
Different from quantum phase diagram, timeofflight absorp tion spectra [5,13] reveals the momentum distribution prop erties of a ultracold optical lattice system and can be detected directly in experiment. To offer a direct comparations with possible experimental observations, as an example, we inves tigate the timeofflight picture of spinor1 Bose gases in a honeycomb optical lattice in antiferromagnetic cases. This can be realized by trapping ultracold 23 Na [3] atoms in opti cal lattices, with the scattering lengths for different channels being ≈ a a 50 0 B and ≈ a a 56 2 B respectively [48,49] (a B being the Bohr radius).
As is known, the momentum distribution function of the system is [33,36] k r r i (36) By making use of (4) and the Fourier transform of the operatorŝˆˆ † † the momentum distribution function can be expressed as (38) where ( ) w k is the Fourier transformation of the Wannier func tion ( ) w r . With the help of the oneparticle Green's function , it can be rewritten as Since we have approximated , now the problem of calculating the density distribution (39) is transformed to calculate the corresponding Green's function . The dots represent the perturbative mean field results [47].
, which can be obained via the following inverse Fourier transformation of (23) According to equation (23), one needs to determine the explicit form of J and U 0, 2 first. Experimentally, a honeycomb optical lattice can be created [29] by three laser beams with wavelength λ ≈ 830 nm intersecting with angles of / π 2 3, lead ing to the lattice constant / λ = a 2 3 3. The optical trapping potential of the lattice is Under the harmonic approximation, the potential ( ) V r can be expanded around the minimal point ( )  (42) and the corresponding Wannier function takes the form To create a 2D system in real experiment, people first cre ate a 3D stacked optical lattice system and then set trapping potential in zdirection extremely stronger than that in x and ydirection. In this case the influence of the third direc tion can be avoided and the system shows twodimensional properties [29]. Plugging (43) into equation (7), and taking = V E 60 z R (which would be deep enough in practice), we have˜( and˜( a a a V E 4 6 9 60 .
According to the preceding discussion, while increas ing the hopping amplitude or equivalently decreasing the lattice depth in experiment, in an antiferromagnetic sys tem the component with α = 0 and m = 0 will go from Mott insulating phase to superfluid state first. Therefore, in this paper we will concentrate on the case with α = 0 and m = 0 for antiferromagnetic system. Without causing misleading, the index of α will be dropped in abbreviation in the following.
For antiferromagnetic interaction, from equations (45) and (46), it is not difficult to verify that ≈ U U 0.04 2 0 for 23 Na. In order to keep the system to stay in the region of second order phase transition, according to [22], we will keep the occupa tion number n = 3 in the following calculation.
For the timeofflight processes released from the ground state ⟩ | n 1, , the explicit form of the cumulant is With the help of the residue theorem and general Leibniz rule, the above equation (48) can be calculated out analytically, yet very tediously. Combining the momentum distribution func tion (39) with (20), (21), (43), (44), (45), in figure 3 we plot the timeofflight absorption pictures of the system released from state ⟩ |1, 3 for different lattice depth Ṽ 0 . From the figures we see clearly that in both cases the momentum distribution ( ) α n k exhibits sharp peaks at the reciprocal lattice sites of the honeycomb lattice for suffi ciently shallow lattice depth, this reflects the coherence of the condensate, indicating that the system is in a superfluid state. For deeper lattice depth, the heights of the peaks decay and the widths of the peaks broaden. The sharp peaks disappear entirely for sufficient deep lattice depth, which is the case in Mott insulating phase. The graphs in figure 3 qualitatively show that the system undergoes a phase trans ition as the opti cal lattice potential varies, which is in agreement with our calcul ation of the phase diagram shown in figure 1.

Summary
In this work, by treating the hopping parameter in the spinor Bose-Hubbard model as a perturbation, with the help of cumulant expansion and the resumming technique, we have calculated the oneparticle Green's function of a spin1 Bose system in a honeycomb optical lattice analytically. By use of the resummed Green's function, we determined the quantum phase diagrams of the system in both ferromagnetic cases and antiferromagnetic cases. We found that in antiferromagnetic cases the Mott insulating states with even filling factor are more robust against the hopping parameter than cases with odd filling factor. This is in agreement with results via dif ferent approaches [28,47]. It should be pointed out that although the lowest order resummed Green's function is a sort of meanfield treatment, it is shown in section 3 that the developed method here is ready to go beyond mean field by adding loop diagrams which consists of higher order cumu lants. Moreover, according to equation (26), it is also straight forward to investigate the finite temperature properties of the system via the resummed Green's function.
The momentum distribution function of the system, meas ured directly in experiment via timeofflight technique by changing the depth of the optical lattice trapping potential, can also be exposed by the Green's function. We therefore have also calculated analytically the timeofflight absorption pictures of the honeycomb optical lattice spin1 Bose system in antiferromagnetic cases. It is clearly shown that the time offlight absorption pictures dramatically change when cross ing the phase boundaries by changing the depth of the optical lattice trapping potential. Our results may serve as a reference object for further studies. 1 [34]. With these in hand, by making use of the commutation relations among ˆ † α a and spin operators, after a laborious yet straightforward calculation, we have the matrix elements of creation and annihilation operators