Two supersolid phases in hard-core extended Bose-Hubbard model

The effect of the next-nearest-neighbor (nnn) tunneling on the hard-core extended Bose-Hubbard model on square lattices is investigated. By means of the cluster mean-field theory, the ground-state phase diagrams are determined. When a modest nnn tunneling is introduced, depending on its sign, two distinct supersolid states with checkerboard crystal structures are found away from half-filing. The characters of various phase transitions out of these two supersolid states are discussed. In particular, for the case with kinetic frustration, the existence of a half supersolid phase possessing both solid and unconventional superfluid orders is established. Our work hence sheds light on the search of this interesting supersolid phase in real ultracold lattice gases with frustrated tunnelings.


I. INTRODUCTION
Ultracold atomic gases in optical lattices provide an extremely versatile and well-controlled experimental platform for exploring complex phenomena in quantum matters [1][2][3]. Owing to remarkable control over microscopic model parameters, they can simulate the physics of strongly correlated systems in regimes which are not easily accessible in condensed matters. The rapid advances in experimental techniques open the door to the search for exotic phases and in probing their quantum critical behaviors. One of the most striking experiments is the observation of the superfluid-Mott transition of ultracold bosons in an optical lattice [4], which can be described by the seminal Bose-Hubbard model [5,6]. This quantum phase transition is driven by the competition between the kinetic energy of particles hopping between lattice sites and the on-site repulsive interaction.
Recently, by using the lattice shaking technique, the effective hopping parameters can be tuned to be negative or even complex [3,[7][8][9][10]. As a result, lattice bosons with frustrated hoppings and/or coupled to artificial gauge fields have now been realized, which can give rise to unconventional superfluid (SF) phases [11,12]. Besides, intriguing supersolid (SS) states have been predicted to occur as well in the Bose-Hubbard model with hopping frustration induced by staggered flux [13].
The SS phases are fascinating quantum states that have superfluid and solid long-range orders simultaneously. Such phases usually come out around the phase boundary of commensurate solid phases and behave as intermediate states between the superfluid and the crystalline ones. In addition to hopping frustration induced by the gauge fields [13], ultracold lattice gases with dipole-dipole interaction [14][15][16] or the spin-orbit coupling [17][18][19] also provide possible candidates in simulating such quantum phases. In particular, the existence of the SS phases has been established in extended Bose-Hubbard models (i.e., Bose-Hubbard models with additional intersite interactions) [20][21][22][23][24][25], which can be naturally realized with dipolar bosons. Due to intersite repulsions among lattice bosons, insulating solid states at commensurate filling factors can be stabilized. Under doping particles or holes into these insulating solid states, the SS states arises when dopants delocalize and undergo Bose-Einstein condensation rather than a phase separation. The late realization of the extended Bose-Hubbard model with nearest-neighbour repulsion for magnetic erbium atoms lays the groundwork in searching for the SS states [26]. Very recently, the SS phases are observed experimentally in atomic quantum gases trapped in an optical lattice inside a high-finesse optical cavity [27,28]. They are stabilized by the long-range interaction mediated by a vacuum mode of the cavity.
More exotic states of matter are expected when one incorporates the hopping frustration into dipolar lattice gases. To uncover the interesting physics of such systems, an extended Bose-Hubbard model with a frustrated nextnearest-neighbor tunneling on square lattices is investigated recently [29]. In the hard-core limit, the grandcanonical Hamiltonian reads where i, j indicates the nearest-and i, j nextnearest neighbors. The destruction (creation) operators of hard-core bosons on site i are denoted by b i (b † i ). The parameters t and t stand for the nearest-neighbor (nn) and the next-nearest-neighbor (nnn) hopping integrals, respectively. The repulsion between the nn sites is represented by V , n i = b † i b i is the number operator of bosons, and µ is the chemical potential. Here both t and V > 0 are taken to be positive, while t can be either positive or negative. When t < 0, the hopping frustration appears.
In Ref. [29] where the frustrated case with t < 0 is considered, by using the tensor network state method [30][31][32], a peculiar half supersolid (HSS) phase is found to appear away from half-filling for large V and modest |t |. Its stabilization results from the interplay of the frustrated nnn hopping and the nn repulsion. This HSS state is characterized by the preferential superfluid flow along the arXiv:1703.05866v2 [cond-mat.quant-gas] 7 Nov 2017 nnn bonds and nonzero Bose-Einstein condensate mainly on one sublattice. In addition, the particle density shows a checkerboard long-range order. The lattice's translational symmetry is thus spontaneously broken. A phase separation regime consisting of the uniform SF and the HSS phases is observed in their fixed-density models [29]. This implies a first-order phase transition between these two phases.
Actually, the non-frustrated nnn hopping (t > 0) can support as well a similar supersolid phase possessing both the SF and the checkerboard crystalline long-range orders [33]. The appearance of such a checkerboard supersolid (CSS) state for the postive-t case has been demonstrated by means of quantum Monte Carlo (QMC) simulations. In contrast to the the negative-t case, the SF-CSS phase transitions are however found to be of secondorder [33], instead of being discontinuous.
In this paper, we focus on the detailed ground-state phase diagrams of the model in Eq. (1) for both signs of t and elucidate the nature of the phase transitions therein. The cluster mean-field theory is employed here, which has been applied to various spin [34][35][36][37][38] and boson [39][40][41][42][43][44][45][46][47] systems with success. We note that the main difference between two supersolid states, the HSS and the CSS states, comes from the distinct momentum states at which bosons condenses. Depending on the sign of t , Bose condensate in the HSS phase occurs at the state with momenta k = (0, π) and (π, 0), but at k = (π, π) in the CSS phase. Schematic pictures of these two phases are illustrated in Figs. 1(a) and (b). By measuring the types of condensation, both supersolid states can thus be identified. Besides the uniform SF and the supersolid phases, there exists as well an incompressible checkerboard solid (CBS) state at half-filling [see Fig. 1(c)]. This solid state breaks as well the lattice's translational symmetry and is signaled by a plateau in the density profile versus chemical potential. Various transitions between these quantum phases are summarized in Fig. 2, which is the main result of this work.
For weak nnn hopping t , direct first-order transitions between the CBS and the SF phases are observed. Upon increasing |t |, an intermediate supersolid state (either CSS or HSS) emerges. In accordance with previous results [29,33], we find that the SF-HSS transitions are of first order, while the SF-CSS transitions are continuous. In contrast, both transitions from each supersolid phase to the CBS phase are shown to be continuous and belong to the superfluid-insulator universality class [5]. According to our CMF results, the stability of the two supersolid phases is enhanced when the nn repulsion V and the magnitude of the nnn hopping t are increased. Since the model under consideration lies within the reach of present experimental setup, hopefully our predictions may be realized in the near future. This paper is organized as follows: In Sec. II, the cluster mean-field theory is briefly described and the measured order parameters are introduced. Our numerical results are presented and discussed in Sec. III. We summarize our work in Sec. IV.
In the first subsection, we explain how the ground states are obtained under the CMF approach. In the second subsection, the ways in determining the order parameters and thus the ground-state phase diagrams are described. The method of the cluster-size scaling is discussed at the end.

A. ground states in the CMF theory
In the CMF theory, the whole system is decomposed into a cluster of lattice sites and its surroundings, as illustrated in Fig. 3. The interaction between these two parts is approximated by the mean fields at the borders of the cluster. The effective Hamiltonian of the cluster is then given by where H C is the exact Hamiltonian of the cluster C and H ∂C contains the effects of the mean fields acting only on sites at the boundary ∂C. By using the standard meanfield decoupling, the latter takes the form for our model in Eq. (1), where the prime over indicates that the site i locates on the boundary ∂C and the site j is outside the cluster. Such mean-field decoupling amounts to the approximation of the many-body wave function being written as a product of the cluster and the environment wave functions.
After diagonalizing the effective Hamiltonian H eff C , the mean fields of the surroundings, b j and n j , can be obtained self-consistently from the corresponding expectation values inside the cluster. We note that the strategy in determining the values of mean fields is not unique. The implementation suggested in Refs. [40] and [41] is adopted here. A sublattice structure expected to emerge in the parameter range is first imposed, into which the cluster of N C sites is embedded. To be compatible with all the expected phases in our system (cf. Fig. 1), states with a four-sublattice structure are considered in the present work. Assuming the mean fields b j to be real, there are thus eight mean-field parameters, b j and n j , in total. Under this four-sublattice pattern, a cluster of N C = 4 × 4 sites is illustrated in Fig. 3. The mean fields to which the cluster is coupled are then calculated by averaging over all sites belonging to the same sublattice. In some cases, this prescription can prevent possible numerical instability happened in other proposals [38].
Given the system parameters t, t , V , and µ, the selfconsistent procedure for determining the ground state includes the following key steps: If the net deviation is smaller than the given tolerance, that is, , then output |GS as the ground state, b j and n j as the self-consistent mean fields. Otherwise, set b j = b j , n j = n j and return to step (i).

B. order parameters and phase diagram
To characterize different phases, several order parameters are calculated. Nonzero Bose condensate indicates the presence of the SF long-range order. The condensate density ρ 0 at a momentum k can be evaluated through the Fourier transform of the mean fields b j , where N C is the total number of sites in the cluster C. In a uniform SF phase, we have Bose condensate at k = (0, 0) with ρ 0 (0, 0) = 0 and ρ 0 (k) = 0 for all other k's. Instead, unconventional SF orders are found in the CSS and the HSS phases. From the inspection of their condensation patterns depicted in Fig. 1, one expects that both ρ 0 (0, 0) and ρ 0 (π, π) are nonzero in the CSS phase, while only ρ 0 (π, 0) = ρ 0 (0, π) can be finite in the HSS phase.
On the other hand, long-range solid order appears if there exists density modulation with a commensurate wave vector k. This can be detected by nonzero values of the Fourier transform of the local densities n j , Its value at k = (0, 0) is nothing but the averaged density. In the thermodynamic limit, this quantity has a direct relation, |ñ(k)| 2 = S(k)/N C , to the static structure factor S(k) ≡ (1/N C ) i,j∈C n i n j e ik·(ri−rj ) . For systems with the CBS order,ñ(π, π) = 0 is expected. In particular, the classical picture of the incompressible CBS state at half-filling predictsñ(π, π) = 1/2. Besides, in order to determine possible first-order phase transitions signaled by energy level crossing [48], we need to measure the ground-state energies of our systems. In the CMF theory, the expression of the energy per site is where E G is the lowest energy eigenvalue of the effective Hamiltonian H eff C . We note that double counting terms coming from the mean-field decoupling of two-cluster interactions must be subtracted. This explains the occurrence of the second term in Eq. (6).
The phase diagrams presented in Fig. 2 are obtained by means of the CMF approach using 4×4 clusters. We note that the phase boundaries will be slightly affected upon increasing the cluster sizes, since longer-ranged manyparticle correlations can be incorporated. By applying a cluster-size scaling, the phase transition points for an infinite lattice can be extrapolated [40,41,43]. In this work, we take the scaling parameter λ = N B /(N C × z/2) introduced in Refs. [40] and [41]. Here, N B denotes the number of the nn bonds within the cluster and z is the coordination number for the lattice (z = 4 for the present case of square lattice). Because the denominator represents the number of the nn bonds of the original lattice per N C sites, the parameter λ provides an indication of how much the correlation effects between particles are considered in the cluster. The infinite-size limit thus corresponds to λ = 1. In the next section, linear extrapolations toward λ = 1 for some typical phase transition points are performed to illustrate the cluster-size dependence of our results.

III. NUMERICAL RESULTS AND DISCUSSIONS
A. non-frustrated case of t > 0 In order to provide benchmarks on the performance of the CMF method, we begin with the non-frustrated case of t > 0, where some QMC data are available in the literature [33]. Fig. 4(a) shows our CMF results of various order parameters as functions of the chemical potential µ for V /t = 3.5 and t/t = 0.01 obtained by using 4 × 4 clusters. For this special case of t t , upon decreasing µ, the incompressible CBS state withñ(π, π) = 1/2 and ρ 0 (0, 0) = ρ 0 (π, π) = 0 transits continuously into the CSS state with coexistence of the CBS and the SF orders. This continuous melting transition occurs at µ c /t 4.0. Our findings are in accordance with those calculated by the QMC approach (see Fig. 3 in Ref. [33]), in which the critical value µ c /t 4.2 is observed. As seen from the inset, the values of this transition point will be shifted to larger ones when the cluster size is increased. Excellent agreement can be reached after a linear extrapolation toward the scaling parameter λ = 1. As a further benchmark, the CMF results using 4 × 4 clusters for V /t = 3.2 and µ/t = 2.9 are presented in Fig. 4(b). By increasing the nn hopping t, the uniform SF state will be eventually stabilized, such thatñ(π, π) vanishes and ρ 0 (k) is nonzero for k = (0, 0) only. The CSS-SF transition is found to be as well continuous and it occurs at t c /t 0.34. They are consistent with the QMC results (see the lowest panel of Fig. 4 in Ref. [33]), where t c /t 0.27 is found. Again, excellent agreement can be obtained after a linear extrapolation toward the scaling parameter λ = 1, as shown in the inset of Fig. 4(b). All these results validate the CMF approach to the present hard-core boson problem. Here t = 1 is taken as energy unit. In both cases, ρ0(π, 0) = ρ0(0, π) = 0. Insets in the panels show the cluster-size scalings of the CMF data for the phase transition points µc. The scaling procedure is discussed in Fig. 4. Lower panels show the details of the order parameters around the SF-CSS and the CSS-CBS transition points in panel (b) with enlarged scales. Lines in these panels are (c) the square-root and (d) the linear fits of the data.
The requirement of t < t studied in Ref. [33] is hard to be satisfied in most cold atom experiments. One may wonder if similar physics remains in the more practical situation such that the nn hopping is larger than the nnn one (i.e., t > t ). In the rest of this paper, only weak t cases are considered, and we take t = 1 as the energy unit.
In Fig. 5, the CMF results using 4 × 4 clusters for two typical values of t = 0.1 and 0.5 with fixed V = 8 are presented. For the t = 0.1 case shown in panel (a), a direct SF-CBS transition is observed, at which the SF and the CBS order parameters, ρ 0 (0, 0) andñ(π, π), change discontinuously. This first-order transition occurs at µ c 1.14. As seen from the inset, the transition point approaches to µ c 1.49 in the infinite-size limit (i.e., λ → 1). The absence of the supersolid phase for this weak t case is expected. In the t = 0 limit, it has been known that the hard-core boson model [49][50][51] (or the equivalent spin-1/2 XXZ model [52,53]) does not stabilize the supersolid phase on a square lattice and the direct SF-CBS transition is of first order. Turning on a small t would not change the whole story, as shown by our results.
According to the Ginzburg-Landau-Wilson paradigm, a direct transition between the SF and the CBS states, which break different symmetries, must be of first order. Otherwise, an intermediate supersolid phase with coexistence of SF and solid long-range orders will emerge in between. The occurrence of the CSS phase between the SF and the CBS phases at modest nnn hoppings has been established in Fig. 5(b). For t = 0.5, within our CMF results using 4 × 4 clusters, the CSS phase is stabilized when 1.46 µ 2.80. The transitions at these two phase boundaries are both continuous. As seen from the inset, the stability region of the CSS phase shrinks but remains finite in the infinite-size limit (i.e., λ → 1).
Two transitions out of the CSS phase have distinct characters. The SF-CSS transition is expected to be of Ising type, since the two-fold sublattice symmetry broken in the CSS state due to the CBS order will be restored after this transition. Around this critical point, our findings of the CBS order parameterñ(π, π) show a square-root dependence on µ,ñ(π, π) ∝ (µ − µ c ) 1/2 , as presented in Fig. 5(c). That is, the critical exponent has a mean-field value β = 1/2. This reflects the meanfield character of our approach. On the other hand, the CSS-CBS transition should belong to the SF-insulator universality class [5], such that the condensate fraction vanishes linearly around the transition point. Our results of both condensate fractions ρ 0 (0, 0) and ρ 0 (π, π) shown in Fig. 5(d) are in agreement with this prediction.
The phase diagrams for the extended Bose-Hubbard model in Eq. (1) with t > 0 is depicted in Figs. 2(a) and (c). We find that the CSS phase exists in extended regions of the phase diagrams. From our CMF calculations, besides the t t limit investigated in Ref. [33], this phase can appear as well for weak t cases as long as strong nn repulsion V is present.

B. frustrated case of t < 0
Now we turn to the frustrated case of t < 0. Due to the negative-sign problem caused by the hopping frustration, QMC simulations are no longer applicable. Therefore, to obtain quantitative predictions of the present frustrated systems, other theoretical approaches, such as the tensor network state method employed in Ref. [29], are necessary. With success in the t > 0 case presented above, we continue our CMF studies on the t < 0 side.
Our CMF results using 4×4 clusters for two typical values of t = −0.1 and −0.3 with fixed V = 8 are presented in Fig. 6. Similar to the cases of weak positive t , a direct first-order transition between the SF and the CBS phases is observed for the t = −0.1 case [see Fig. 6(a)]. This SF-CBS transition occurs at µ c 0.82 and the value of the transition point approaches to µ c 1.09 in the infinite-size limit (i.e., λ → 1).
Upon increasing the magnitude of the negative t , an intermediate HSS phase with coexistence of SF and solid long-range orders can emerge in between. Similar to the CSS phase in the non-frustrated case, the HSS state is characterized as well by the preferential superfluid flow along the nnn bonds. Nevertheless, Bose condensate in the HSS phase occurs at the state with momenta k = (0, π) and (π, 0), instead of k = (π, π) in the CSS phase. Here t = 1 is taken as energy unit. In both cases, ρ0(π, π) = 0 and ρ0(0, π) = ρ0(π, 0). Insets in the panels show the cluster-size scalings of the CMF data for the phase transition points µc. The scaling procedure is discussed in Fig. 4. Panel (c) shows the groundstate energies e as a function of µ calculated by Eq. (6)  The case of t = −0.3 is shown in Fig. 6(b), in which the HSS phase is stabilized when 0.41 µ 1.05 within our CMF results using 4 × 4 clusters. As seen from the inset, the stability region of the HSS phase shrinks but remains finite in the infinite-size limit (i.e., λ → 1).
In agreement with the SF-insulator universality class [5], we find that the condensate fraction ρ 0 (π, 0) vanishes linearly around the continuous HSS-CBS transition [see Fig. 6(d)]. However, the phase transition between the uniform SF and the present supersolid phases is found to be of first order, instead of being continuous in the non-frustrated case of t > 0. The discontinuous jumps in the order parameters at the SF-HSS transition point is a clear indication that phase separation would occur in a canonical system, as observed in Ref. [29].
First-order quantum phase transitions come from energy level crossing in the ground state and are thus signaled by a cusp in the ground-state energy curve [48]. In Fig. 6(c), the energy curve in the ground state around the discontinuous SF-HSS transition is shown, where the location of the cusp gives the transition point we want. In the present work, all the first-order transition points, including those of the direct SF-CBS transitions shown in Figs. 5(a) and 6(a), are determined by this way.
Distinct characters of the SF-CSS transition for the t > 0 case and the SF-HSS transition for the t < 0 case can be understood from different condensation patterns in these two supersolid states. Bose condensate at the state with momenta k = (π, 0) and (0, π) in the HSS state gives b j ∝ e i(π,0)·rj +e i(0,π)·rj = cos(πx)+cos(πy). This leads to alternating phase change in π among sites within the sublattice on which bosons condense, as displayed in Fig. 1(b). Thus a smooth evolution from the HSS state to the uniform SF state should be unlikely. On the contrary, there is no phase difference in the condensate for the CSS state [see Fig. 1(a)]. A continuous transition from the CSS state to the uniform SF state is thus expected. The phase diagrams for the extended Bose-Hubbard model in Eq. (1) with t < 0 is presented in Figs. 2(a) and (b). The existence of the HSS phase in extended regions of the phase diagrams is demonstrated. Indicated by our CMF calculations, the stability of the HSS phase is enhanced when the nn repulsion V and the magnitude of the nnn hopping t are increased. As seen from the phase diagrams in Fig. 2, the HSS phase for t < 0 exists in larger parameter regions than the CSS phase for t > 0 does. This implies that, in stabilizing the supersolid phase, the effect of frustrated nnn hopping is more significant. In particular, when t /t < −1/2, the uniform SF phase becomes absent and the HSS phase can exist for all incommensurate fillings. This can be understood from the single-particle dispersion relation, which has local minima at k = (0, π) and k = (π, 0) when t /t < −1/2. This thus favors the Bose condensation consistent with the SF order in the HSS state, rather than the uniform SF one.

IV. SUMMARY AND CONCLUSIONS
Using the cluster mean-field theory, we explored in detail the formation of two supersolid phases, which arise in the hard-core extended boson Hubbard model of Eq. (1). We find that, in the presence of a large nn repulsion, by introducing a modest nnn hopping t , either the CSS or the HSS phases can be stabilized away from halffilling. Both supersolid phases come out around the phase boundary of the CBS phase and behave as intermediate states between the SF and the CBS states. Despite carrying the same checkerboard crystalline long-range order, bosons condense at different momentum states in these two supersolid phases, as shown in our calculations.
The characters of various phase transitions out of the CSS or the HSS states are discussed. Their transitions to the CBS state are found to be continuous and belong to the SF-insulator universality class [5]. However, for the frustrated case with t < 0, the SF-HSS transition is shown to be of first order, instead of being continuous observed at the SF-CSS transition for the t > 0 case. The distinct nature of these two phase transitions can be understood from their different condensation patterns. While only the model of hard-core bosons is considered in this work, similar phenomena are expected as well when the hard-core constraint is replaced by a finite but strong on-site repulsion. The latter condition can be accomplished, for example, by tuning the s-wave scattering length via a Feshbach resonance [54,55]. As mentioned in the Introduction, the realization of the extended Bose-Hubbard model has been achieved recently [26] and frustrated tunnelings can be produced by, say, lattice shaking techniques [3,[7][8][9][10]. That is, systems related to the model discussed here will be hopefully soon available experimentally. Our results thus provide a useful guide to the experimental search of these supersolid phases, especially the interesting HSS phase in real ultracold lattice gases with frustrated tunnelings.