Mixed-order transition and tricritical point associated with checkerboard supersolidity in the two-dimensional $t_2-V_1$ model
Amrita Ghosh Satyaki Kar Sudhakar Yarlagadda





^1 Department of Physics, Ben-Gurion University of the Negev, Beer-Sheva 8410501, Israel.
^2 Institute of Physics, Bhubaneswar, India.
^3 CMP Division, Saha Institute of Nuclear Physics, HBNI, India


We use Quantum Monte Carlo method employing stochastic-series-expansion technique to study the ground state properties of the t_2-V_1 model on a square lattice. We find that, away from half-fillings, the minimal combination of nearest-neighbor  repulsion V_1 and next-nearest-neighbor hopping t_2 may give rise to checkerboard supersolidity.
  The nature of the quantum phase transition, where  the superfluid  changes to a checkerboard supersolid,
  depends on the relative strength of V_1/t_2 and the average site occupancy. Interestingly, the model exhibits a mixed-order transition near half filling; 
  at a higher (lower) filling,  tricriticality  is witnessed followed by a second-order transition at densities even further away from half filling.
  Close to half filling, the model displays the extreme Thouless effect and transits from a superfluid to a checkerboard solid.   

Mixed-order transition  and tricritical point associated with checkerboard supersolidity  in
the two-dimensional t_2-V_1 model
    Sudhakar Yarlagadda^3
    Received: date / Accepted: date
==============================================================================================================================



§ INTRODUCTION

Ultracold atoms in optical lattices offer a novel platform
where generic many-body phenomena of
condensed matter physics can be explored<cit.>. 
One can can study exotic quantum phases such as  lattice supersolidity, which is the homogeneous
coexistence of superfluidity/superconductivity and crystal order in discrete lattices.
In fact, one such example is the recent experimental realization of supersolid phases due to competing
short-range and infinite-range interactions for bosonic atoms in
optical lattices<cit.>.  The
infinity-range interactions are generated by a vacuum mode of
the cavity and can be independently manipulated<cit.>.
The optical-lattice platform has been quite successful in simulating  models  for strong correlations such
as the Bose-Hubbard models<cit.> which captures the 
essential physics of lattice supersolidity as shown by various theoretical 
works on two-dimensional square lattices<cit.>.

Besides in artificially engineered optical lattices, lattice supersolidity can also
occur in naturally formed systems. Lattice supersolidity has been observed in a variety
of systems such as  three-dimensional doped bismuthates <cit.>,  quasi-two-dimensional doped
dichaclcogenides<cit.> and molecular crystals<cit.>, and 
quasi-one-dimensional doped trichalcogenides<cit.> and doped spin-ladder systems<cit.>. 


In an earlier communication<cit.>, it was shown that the t_1-t_2-t_3-V_1 model on a square lattice can produce (π,π) (or checkerboard)
supersolidity
when the same-sublattice tunneling is sizeable and  nearest-neighbor repulsion is large.
In fact, the purpose of the present paper is to demonstrate that the t_2-V_1 model is the minimum model for checkerboard supersolidity 
and to elucidate the rich physics manifested by this model.

The study of t_2-V_1 model dates back to constructing effective Hamiltonians
arising due to the presence of  cooperative breathing modes observed in
oxide systems. 
Many oxides such as  manganites<cit.>, cuprates<cit.>, and bismuthates<cit.> show evidence of cooperative
strong electron-phonon interactions.
By including  cooperative strong electron-phonon couplings in a Holstein model, an effective Hamiltonian was obtained where the dominant transport comes from double hopping
and the dominant repulsion is between nearest neighbors both in one dimension <cit.> and in two dimensions <cit.>.

Phase transitions 
classified by the Ehrenfest scheme are first order, second order, 
and higher order.
In a first-order transition, the order parameter jumps, whereas  its fluctuations are not large on either side of the transition.
Furthermore, coexistence and hysteresis are some of the
usual features associated with this  transition.
On the other hand, a second-order transition is characterized by a lack of discontinuity and 
anomalously large fluctuations of the order parameter.
Contrastingly, Thouless<cit.> 
found mixed-order transition indicated by  order parameter
jumping and displaying large fluctuations. Subsequently,
 several systems with such mixed-order transition have been 
 reported <cit.>. Interestingly the term “extreme Thouless effect”
 (i.e., an extreme version of this
 mixed-order)   ,
was coined  to denote a 
transition where both the jump and the fluctuations are maximal<cit.>.
In this paper, we report another instance of the extreme Thouless effect
 in the context of a minimal model (i.e., the t_2-V_1 model) for checkerboard supersolidity.



The outline of the paper is as follows.
In section II, we discuss our numerical techniques and details of the calculations. Section III describes the results obtained and the following discussions.
Lastly in section IV, we summarize our work and discuss briefly its novelty.




§ FORMULATION


To study the various phases of the two-dimensional t_2-V_1 model for HCBs, we employ stochastic-series-expansion (SSE) technique<cit.>, a quantum Monte Carlo (QMC) method, involving directed loop updates<cit.>. By identifying b_i,j^†=S_i,j^+, b_i,j=S_i,j^- and n_i,j=S_i,j^z+1/2, to employ SSE, we recast the HCB Hamiltonian in terms of spin-1/2 operators. The converted Hamiltonian, in terms of 2t_2, takes the form of an extended XXZ Hamiltonian given by

    H=   -∑_i,j1/2(S_i+1,j+1^+ S_i,j^-+S_i-1,j+1^+ S_i,j^-+H.c.)
       +∑_i,jΔ_1 (S_i,j^z S_i+1,j^z+S_i,j^z S_i,j+1^z)-h∑_i,j S_i,j^z

with Δ_1=V_1/(2t_2). In the above equation we have introduced the term -h∑_i,j S_i,j^z, where h is a variable and can be thought of as a dimensionless external magnetic field acting on the system. By tuning this variable we can access different magnetizations or filling-fractions of the system.

In Ref. Ghosh1, the one-dimensional t_2-V_1 model was shown to undergo a discontinuous transition from a superfluid phase, with equally populated sublattices, to a checkerboard supersolid state (with coexisting CDW and superfluidity), where all the particles occupy a single sublattice. Hence, to study the competition or coexistence of these two diagonal (CDW) and off-diagonal (superfluidity) long-range orders, in the two-dimensional version of the t_2-V_1 model, we use two order parameters : structure factor S(Q⃗) and superfluid density ρ_s. The structure factor per site is expressed as

    S(Q⃗)=1/N_s^2∑_i,j∑_l,me^iQ⃗·(R⃗_i,j-R⃗_l,m)⟨ S^z_i,j S^z_l,m⟩,

where ⟨⋯⟩ denote the ensemble average and N_s is the total number of sites of the system. We study S(Q⃗) at all possible values of wavevector Q⃗ and identify the ones for which the structure factor shows peaks.
For Q⃗=(π,π), the structure factor takes the form

    S(π,π)=1/N_s^2∑_i,j∑_l,m(-1)^(i-l)(-1)^(j-m)⟨ S^z_i,j S^z_l,m⟩,

which, in terms of number operator n_i,j, can be re-expressed as

    S(π,π)=1/N_s^2∑_i,j∑_l,m   (-1)^(i-l)(-1)^(j-m)
       ⟨(n_i,j-1/2) (n_l,m-1/2)⟩,

If for a site (i,j), the sum (i+j) is even then we call it an even site, otherwise it is called an odd site. A square lattice with even number of sites can always be divided into two equal sublattices: even sublattice containing all the even sites and odd sublattices which contains the odd sites. We define the number operators giving the total number of HCBs at even and odd sites as N̂_̂ê=∑_i+j=even n_i,j and N̂_̂ô=∑_i+j=odd n_i,j respectively. Now, the summation in Eq. <ref> can be divided into two parts based on the fact that there are two possible scenarios; either both the sites (i,j) and (l,m) belong to the same sublattice or they belong to two different sublattices. Noting that (-1)^(i-l)(-1)^(j-m) takes the value +1(-1) when (i,j) and (l,m) belong to the same sublattice (two different sublattices), Eq. <ref> can be written as


    S(π,π)=1/N_s^2   ⟨ [N̂_̂ê^2+N̂_̂ô^2-N_s/2(N̂_̂ê+N̂_̂ô)+N_s^2/8 ] ⟩
    
     -1/N_s^2   ⟨ [ 2N̂_̂êN̂_̂ô-N_s/2(N̂_̂ê+N̂_̂ô)+N_s^2/8 ] ⟩ ,

which reduces to

    S(π,π)=1/N_s^2⟨ (N̂_̂ê-N̂_̂ô)^2 ⟩.


Since the Hamiltonian consists only NNN hopping, both N̂_̂ê and N̂_̂ô commute with the Hamiltonian H; hence we obtain

    S(π,π)=1/N_s^2(N_e-N_o)^2,

where N_e (N_o) denotes the total number of HCBs at even (odd) sites. Thus, when both the sublattices are equally occupied, i.e., N_e=N_o, the structure factor at wavevector (π,π) attains its minimum value, S(π,π)_ min=0. On the other hand, when all the particles occupy only one sublattice, i.e., either N_e=N_p (N_p being the total number of particles in the system) and N_o=0 or vice-versa, the structure factor maximizes to 

    S(π,π)_ max=(N_p/N_s)^2=ρ^2,

where ρ=N_p/N_s is the filling-fraction of the system. In terms of magnetization m of the system, the maximum value of the structure factor reduces to

    S(π,π)_ max=(m+1/2)^2.

In terms of winding numbers along x and y directions, W_x and W_y, the superfluid density can be expressed as

    ρ_s=1/2β⟨ W_x^2+W_y^2⟩,

where β denotes the inverse temperature. One can calculate the winding number in the x direction as W_x=1/L_x(N_x^+-N_x^-), where N_x^+ (N_x^-) represents the total number of operators transporting spin in the positive (negative) x direction with L_x being the length of the lattice along the x direction.

As discussed in Ref. Scalettar, for a L× L square lattice, SSE simulations done at low enough temperatures, such that the inverse temperature β∼ L, can capture the ground-state properties of the system. Since numerical calculations for both β=3L/2 and β=2L produce same results (within the error bars of our calculations), in this paper we report the results for β=3L/2 only.
Now, in SSE, a small parameter ϵ is introduced to make sure that all the two-spin matrix elements are positive. However, for large values of anisotropy (i.e., large values of Δ_1) the autocorrelation time can be affected by the value of ϵ. Therefore, it is important to study the autocorrelation time at various regions of the phase diagram to make sure that the bin size used in the calculation of the observables are much larger than the autocorrelation time in all the cases.
To calculate the autocorrelation time we use the formula

    τ_ int[m]=1/2+∑_t=1^∞ A_m(t),

with

    A_m(t)=⟨ m(i+t) m(i)⟩-⟨ m(i)⟩^2/⟨ m(i)^2⟩-⟨ m(i)⟩^2,

where i and t represents the Monte Carlo steps and ⟨⋯⟩ indicates the average over time. Since, in the case of two-dimensional t_2-V_1 Hamiltonian, all the phase transitions occur at apparently lower values of anisotropy, we find that we can safely use ϵ=Δ_1/4, keeping consecutive
bins used in the QMC simulation as independent of the other (with bin-size > the autocorrelation times). We choose the magnetic fields in the vicinity of the phase transitions, where we expect the autocorrelation times to be larger, as well as far away from them. The bin size we have used for our calculations is 7,00,000, which is sufficient to keep the autocorrelation times well within the bin size.


§ RESULTS AND DISCUSSIONS

To construct the phase diagram, we study the system at different anisotropy values by varying the magnetization m from 0 to 0.5. In terms of particle filling-fraction, this corresponds to the variation of particle density ρ from 1/2 to 1. The particle-hole symmetric Hamiltonian forces the physics at any filling fraction for particles to be identical to the one for holes at the same filling. One should note that, since we vary the particle density ρ from 1/2 to 1, the relevant physics should be determined by the holes. Therefore, the maximum value of the structure factor, i.e. S(π,π)_ max reduces to

    S(π,π)_ max=(1-ρ)^2=(1/2-m)^2.





Fig. <ref> displays the variation of the superfluid density ρ_s and structure factor S(π,π) for four different values of Δ_1, as the magnetization m of the system is tuned from 0 to 0.5. For a small value of NN anisotropy, Δ_1=1.10 (see Fig. <ref>(a)), the system manifests superfluidity over the whole range of magnetization, where both the sublattices are equally occupied by the particles, giving rise to a zero structure factor. At a slightly higher value of Δ_1=1.214 (see Fig. <ref>(b)), we find
superfluid (SF) density to undergo a downward kink at around m≈ 0.26. Beyond this point, a S(π,π) order develops in the system
which continuously increases to soon mimic the S(π,π)_ max curve, as given by Eq. <ref> (accompanied by a continuous decrease in the superfluid density). 
The system undergoes a continuous transition from a SF to a checkerboard supersolid (cSS) phase with both the S(π,π) order and superfluidity coexisting homogeneously. As soon as S(π,π) becomes equal to S(π,π)_ max,
we get single sublattice occupation in the ground state. As we further go to higher values of Δ_1, surprisingly the nature of the phase transition changes. For example, for Δ_1=1.22 one can see from Fig. <ref>(c) that initially the system manifests a SF phase, but around m≈ 0.188 the system undergoes striking jumps in order parameters, similar to a first-order phase transition from a SF phase to a cSS phase, pushing all the particles to occupy a single sublattice. 
On the other hand, similar to a second-order transition, there is a large fluctuation in the particle number in a sublattice. In fact, at the
transition there is a large  degeneracy (i.e., equal to approximately the  total number of particles) and the entropy shoots up. 
The degenerate states at this strange phase transition can have 
single-sublattice occupancy ranging from 0 particles to the total number of particles; thus, the fluctuations of  single-sublattice occupancy
is spread all over the
permitted range of variation of the parameter. This transition is the exotic extreme Thouless effect <cit.>. In the context of  the 
one-dimensional t_2-V_1 model dealt with in Ref. Ghosh1, a similar
extreme Thouless effect was detected as can be seen from Figs. 1, 4, and 5 of this work.

At even higher Δ_1 values, such as Δ_1=1.50, the particles form a checkerboard solid (cS) at half-filling (i.e., m=0) with S(π,π)=S(π,π)_ max, where one sublattice is completely filled and the other one is completely empty. As we move slightly away from the half-filling, superfluidity develops and the system retains a cSS order which continues all the way to m=0.5. Now, based on the competition between the hopping and the repulsion term in the Hamiltonian, one could explain Fig. <ref>. It is important to note that, in general, larger NN repulsion V_1 assists the formation of checkerboard solid by forcing the particles to be in a single sublattice, whereas the NNN hopping t_2 helps a particle to hop in the same sublattice. In the case of Fig. <ref>(a), the NN repulsion is not large enough to restrict the particles in one sublattice, rather it is energetically favorable for the system if the particles lower their energy by hopping to different sites thereby giving rise to superfluidity for all filling-fractions. On the other hand, for large repulsion values, such as Δ_1=1.50 in Fig. <ref>(d), at half-filling the particles arrange themselves in alternate sites to avoid NN occupation and thus form a checkerboard solid as depicted in Fig. <ref>. If one extra particle is now added in the system the particle can occupy any one of the empty sites, because no matter which site it resides on it will feel the same amount of extra repulsion 4V_1. Now by the virtue of NNN hopping t_2, this extra particle can hop to its NNN sites with a checkerboard solid in the background without costing any additional energy, thereby resulting in the coexistence of superfluidity and CDW, i.e. a checkerboard supersolid phase. This cSS phase persists all the way up to filling-fraction 1. Now, for the intermediate values of repulsion, such as Δ_1=1.214 and  Δ_1=1.22 in Figs. <ref>(b) and <ref>(c), the repulsion is not strong enough to form a cS at half-filling. Since there are a substantial amount of holes present in the system (for filling-fraction 1/2 and in its vicinity), initially it is energetically favorable for the system to be in the superfluid phase, where the particles can lower their energy by hopping to other sites. In this phase although the particles experience NN



repulsion, the energy lowered by the hopping process is large enough to overcome this energy cost.






The results shown in Fig. <ref> imply that, in between Δ_1=1.214 and Δ_1=1.22, the system must have passed through a tricritical point at which the nature of the phase transition (between SF and cSS) changes from continuous to mixed  order. 
To locate the tricritical point, in Fig. <ref>, we plot the structure factor S(π,π) and superfluid density ρ_s as a function of magnetization m for a number of Δ_1 values.
Let us here describe only the results from the 32×32 lattice.
We observe that for Δ_1=1.17 the system manifests superfluid phase only. As we increase the NN repulsion, the system passes through a continuous phase transition from SF to cSS. For lower values of Δ_1, the system never reaches a state where all the particles occupy a single sublattice, but as the  Δ_1 value goes up the system gets closer to the single sublattice occupancy state. For some NN anisotropy between Δ_1=1.216 and Δ_1=1.217, the nature of the phase transition changes from continuous to mixed order, which corresponds to a tricritical point. Up to Δ_1 = 1.2205 the system passes through a mixed-order transition from SF to cSS. As the value of Δ_1 is further increased, at Δ_1=1.221 the half-filled system manifests a cS phase and goes into a cSS phase in a continuous manner as the density is varied. Although the value of S(π,π) agrees with the maximum possible value of structure factor S(π,π)_ max for densities close to half-filling, it starts to deviate from the S(π,π)_ max curve as we go towards the fully-filled system indicating a deviation from the single sublattice occupancy state. For even larger values of Δ_1 the particles always occupy a single sublattice for all filling-fractions and thus the S(π,π) and S(π,π)_ max curves merge together.
It should be noted that in Fig. <ref>, for the system sizes considered,
we have restricted the magnetization axis from 0.15 to 0.5 only because this is the window where the phase transition takes place. For magnetization values between m=0 and m=0.15, the structure factor S(π,π) is essentially zero for all anisotropy values up to Δ_1=1.2205, whereas it matches with the theoretical S(π,π)_ max curve for Δ_1=1.221 and above.

Fig.<ref> displays the variation of the structure factor S(π,π) for different values of Δ_1 and 
as a function of magnetization m measured on a 24× 24 , 28× 28 and
32× 32 size systems, respectively. Comparing these results one can find that, as we increase the system size the number of mixed-order lines increases.
For a smaller system size one has to go for  smaller repulsion window to identify the tricritical point which makes it  harder to detect.




The complete ground-state phase diagram is displayed in Fig. <ref> for HCBs on 24× 24, 28× 28 and 32× 32 square lattices. In the phase diagram, the dashed (solid) line represents a discontinuous (continuous) transition from superfluid (SF) to checkerboard supersolid (cSS) region as the magnetization of the system is varied, where the filled circles denote the tricritical points. Note that the three different colors, i.e. blue, magenta and green, represent the phase boundaries for the three different system sizes 32× 32, 28× 28 and 24× 24, respectively. Apart from the slight shift of the tricritical point, the phase diagram appears to be by and large independent of the system size. We now discuss the phase diagram for the 32× 32 square lattice.  For anisotropy values Δ_1⩽1.17, the system manifests superfluidity for all values of magnetization. Beyond this point (i.e. Δ_1=1.17)  a small region of cSS phase starts to develop, through continuous phase transition, close to m≈0.5. As we increase the Δ_1 value, the magnetization value at which the system goes into the cSS region shifts towards m=0. 





In the region between Δ_1=1.216 and Δ_1=1.217, which is represented as the filled blue circle in Fig. <ref>, the nature of the transition changes from continuous to mixed order, thereby giving rise to a tricritical point. 
We observe that for Δ_1⩾ 1.221, the system manifests a checkerboard solid (cS) at half-filling
and beyond m=0 a cSS region, persisting all the way to m=0.5, is developed.



Figs. <ref>(a) and <ref>(b)  show the snapshots of the spin configurations barely before and after the discontinuous jump seen in the ground state
corresponding to Fig. <ref>(b). We can clearly see building of checkerboard order (i.e., occupancy along NNN sites or in a single sublattice) beyond this jump. 
Fig. <ref>(c) also shows a configuration [corresponding to m=0 in Fig. <ref>(c)] with full checkerboard order where we obtain S(π,π)=S(π,π)_ max. 

Now, we want to study the nature of the phase transitions as we vary the magnetization of the system for a fixed value of NN anisotropy. As already seen in Fig. <ref> (b),
the order parameters show a continuous variation as a function of the magnetization for Δ_1=1.214 
which signifies continuous phase transition between different phases. On the other hand, the discontinuous jumps in the order parameters in Fig. <ref> (c) for Δ_1=1.22, indicates the existence of a mixed-order phase transition. However, to bring out the nature of the phase transitions along the m axis of the phase diagram, a more reliable process is to study the order parameters, i.e. magnetization, superfluid density and structure factor, as a function of the magnetic field h. Fig. <ref> (a) shows that as we vary the magnetic field h, keeping the NN anisotropy value fixed at Δ_1=1.214, the order parameters change in a continuous fashion. This evidently rules out the possibility of a first-order or mixed-order phase transition and establishes the fact that at Δ_1=1.214, as we move along the m-axis of the phase diagram, the superfluid (SF) and the checkerboard supersolid (cSS) phases are separated by a continuous phase transition. At a higher value of NN anisotropy, Δ_1=1.22, a discontinuous jump in the structure factor accompanied by a sudden drop in the superfluid density is observed as a function of the magnetic filed h (see Fig. <ref> (b)); however, there seems to be no visible jump associated with the magnetization curve. 
Now, zooming into the m-h curve in the region, where the discontinuous jumps in S(π,π) and ρ_s are observed, we see that the magnetization curve shows a very small but sharp drop as the magnetic field is varied. Usually, whenever a first-order phase transition is encountered by the variation of the magnetic field, a sudden upward jump in the magnetization curve is observed which signifies the existence of a phase-separated region. In contrast to the usual scenario, the mixed-order phase transition encountered in the case of Δ_1=1.22, is not associated with any phase-separated region. The sharp drop in the magnetization curve simply indicates a mixed-order transition from a superfluid phase with equal sublattice-occupancy to a checkerboard supersolid phase where only one sublattice is occupied. Next, for Δ_1=1.50. Fig. <ref> (c) depicts continuous variations of the order parameters as the magnetization of the system is tuned. This means that in the phase diagram shown in Fig. <ref>, as we move along the m-axis keeping value of NN anisotropy fixed at Δ_1=1.50, a continuous phase transition, from a checkerboard solid (cS) to a checkerboard supersolid (cSS) phase, is encountered.
The following subsection gives some detailed analysis of the various transitions encountered in the present system.



 §.§ Nature of transitions

To study the nature of the phase transitions encountered while moving along the Δ_1-axis of the phase diagram at a fixed value of magnetization, first it should be noted that, in our simulations, we can not tune the magnetization of the system directly. Instead, we introduce a magnetic field in the system, by tuning which we can access different magnetization values of the system. Since the resulting magnetization for a particular value of magnetic field generally fluctuates during the simulation, it becomes almost impossible to study the nature of the phase transition by varying the Δ_1 value at a fixed value of magnetization. Nevertheless, while in a charge-density-wave (CDW) state the system always shows plateau in the m-h curve, i.e., the magnetization of the system remains unchanged over a range of magnetic field values. Therefore, we can choose any magnetic field lying in the plateau and obtain the desired magnetization value for different values of NN anisotropy.

Fig. <ref>(a) shows that, at 1/2-filling (corresponding to m=0), as the NN anisotropy Δ_1 is varied from 1.0 to 1.5,  the structure factor S(π,π) sharply jumps from 0 to its maximum value 0.25 at Δ_1≈1.221. At the same time the superfluid density ρ_s dramatically drops down to zero. In the phase diagram depicted in Fig. <ref>, as we move along the Δ_1-axis at m=0,
this indicates a mixed-order phase transition from a U(1) symmetry breaking superfluid (SF) to a translational symmetry breaking checkerboard solid (cS) state at Δ_1≈1.221.  
It is important to point out that only when the system is in the CDW state after the phase transition, the magnetization can be fixed at m=0; but before the transition, in the superfluid phase, the magnetization is given by m=0± 0.00000063.

As we have already discussed, the dashed (solid) line in the phase diagram in Fig. <ref>, signifies a  
mixed-order (continuous) transition between the superfluid and checkerboard supersolid phase. Since the nature of the transition between any two phases should be independent of the driving parameters (Δ_1 or m in this case), irrespective of whether we cross the dashed (solid) line horizontally (by varying the magnetization at a fixed Δ_1 value) or vertically (i.e., moving along the Δ_1-axis at a fixed magnetization value) in the phase diagram, the nature of the transition should remain mixed-order (discontinuous). To demonstrate this point, we concentrate on the phase diagram for the 32× 32 square lattice around m≈0.204 and observe that as the Δ_1 value is increased from 1.20 to 1.25, the system goes through a phase transition from SF to cSS. To determine the nature of this phase transition, we vary the magnetic field in very small steps so that we can obtain the magnetization as close as possible to 0.204 for a number of Δ_1 values between 1.20 and 1.25. Fig. <ref>(b) depicts the variation of S(π,π) and ρ_s in terms of Δ_1 at magnetization m=0.20404 ± 0.00001. The sharp jumps in the order parameters clearly indicate that in the phase diagram, as we move along the Δ_1-axis keeping the magnetization fixed at m≈0.204, the phase transition, encountered between the superfluid (SF) and checkerboard supersolid (cSS) phase, is mixed order in nature.

Next, in order to determine the nature of the SF-cSS transition while crossing the solid line vertically in the phase diagram of a 32× 32 lattice, we plot, in Figs. <ref>(c) and <ref>(d), the order parameters S(π,π) and ρ_s as a function of Δ_1 while keeping the magnetization fixed at m=0.33840 ± 0.00003. The continuous variation of the order parameters rules out the possibility of a mixed-order transition. This establishes the point that, irrespective of whether we cross the solid line horizontally or vertically, the nature of the SF-cSS transition remains the same, i.e., continuous.





§ SUMMARY

To summarize briefly, the present work deals with the t_2-V_1 model on a square lattice which turns out to be the minimum model exhibiting checkerboard supersolidity.
This model has a well defined physical origin that goes back to dominant-particle-transport mechanism of double hopping realized in a system with
cooperative normal mode at strong electron-phonon interaction. Alternately, it can always be realized in an optical lattice platform using cold atoms.
The fascinating feature of  this model is its rich ground state phase diagram characterized by various exotic phases and unusual 
quantum phase transitions. It shows a tricritical point,
at an optimum strength V_1/t_2, that separates mixed-order  and continuous transitions involving  (π,π)-checkerboard orders. Importantly, the system displays extreme Thouless effect close to half-filling.



  Authors thank Pinaki Sengupta and Deepak Dhar for valuable comments and suggestions on the work. Though all the numerics is done using the computational cluster at SINP, India, the final communication of the work occurred after AG joined Ben-Gurion University, Israel, under QuantERA project InterPol.






1
ultracold
I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885
(2008).
landig1
R. Landig, L. Hruby, N. Dogra, M. Landini, R. Mottl, T. Donner,
and T. Esslinger, Nature (London) 532, 476 (2016).
landig2
R. Mottl, F. Brennecke, K. Baumann, R. Landig, T. Donner, and
T. Esslinger, Science 336, 1570 (2012).
esslinger
 H. Ritsch, P. Domokos, F. Brennecke, and T. Esslinger, Rev.
Mod. Phys. 85, 553 (2013).
zoller1
D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller,
Phys. Rev. Lett. 81, 3108 (1998).
scalettar1 G. G. Batrouni and R. T. Scalettar, Phys. Rev. Lett. 84, 1599
(2000).
scalettar2
F. Hébert, G. G. Batrouni, R. T. Scalettar, G. Schmid,
M. Troyer, and A. Dorneich, Phys. Rev. B 65, 014513 (2001).
boninsegni1
Long Dang, Massimo Boninsegni and Lode Pollet, Phys. Rev. B 78, 132512 (2008).
zoller2
B. Capogrosso-Sansone, C. Trefzger, M. Lewenstein, P. Zoller, and G. Pupillo
Phys. Rev. Lett. 104, 125301 (2010).
wessel1
Y.-C. Chen, R. G. Melko, S. Wessel, and Y.-J. Kao, Phys. Rev.
B 77, 014524 (2008).
pinaki
P. Sengupta, L. P. Pryadko, F. Alet, M. Troyer, and G. Schmid
Phys. Rev. Lett. 94, 207202 (2005).
troyer1
G. Schmid and M. Troyer, Phys. Rev. Lett. 93, 067003
(2004).
kar
S. Kar and S. Yarlagadda, Ann. Phys. 375, 322 (2016).
Datta
S. Datta, S. Yarlagadda, Solid State Commun. 150,  2040 (2010).
lv2
Xiao Huo, Yong-Yong Cui, Dali Wang, and Jian-Ping Lv, Phys. Rev. A 95, 023613 (2017).
bismuthate3
A. M. Gabovich , A. I. Voitenko and M. Ausloos, Phys. Rep. 367,  583 (2002).
bismuthate2
S. H. Blanton, R. T. Collins, K. H. Kelleher, L. D. Rotter, Z. Schlesinger, D. G. Hinks, and Y. Zheng,
Phys. Rev. B 47, 996 (1993).
quasi_2d_1
R. L. Withers, J.A. Wilson, J. Phys. C 19, 4809 (1986).
quasi_2d_2
J. Merino, R. H. McKenzie, Phys. Rev. Lett. 87, 237002 (2001) 
quasi_1d_1
W. W. Fuller, P. M. Chaikin, N.P. Ong, Phys. Rev. B 24, 1333 (1981).
quasi_1d_2
A. Rusydi, W. Ku, B. Schulz, R. Rauer, I. Mahns, D. Qi, X. Gao, A.T.S. Wee, P. Abbamonte, H. Eisaki, Y. Fujimaki,
S. Uchida, M. Rübhausen, Phys. Rev. Lett. 105, 026402 (2010).
quasi_1d_3
P. Abbamonte, G. Blumberg, A. Rusydi, A. Gozar, P. G. Evans, T. Siegrist, L. Venema, H. Eisaki, E. D. Isaacs,
G. A. Sawatzky, Nature 431, 1078 (2004).
manganites1
A. Lanzara, N. L. Saini, M. Brunelli, F. Natali, A. Bianconi, P. G. Radaelli, and S.-W. Cheong, Phys. Rev. Lett. 81, 878 (1998).
shen
A. Damascelli, Z. Hussain, and Z. X. Shen, Rev. Mod. Phys. 75, 473 (2003).
tvr
A. Taraphder, R. Pandit, H. R. Krishnamurthy, and T. V.
Ramakrishnan, Int. J. Mod. Phys. B 10, 863 (1996).
ravindra
R. Pankaj and S. Yarlagadda, Phys. Rev. B 86, 035453 (2012).
Ghosh1
A. Ghosh and S. Yarlagadda, Phys. Rev. B 90, 045140 (2014).
ghosh2
A. Ghosh, S. Yarlagadda, Phys. Rev. B 96, 125108 (2017).

thouless
D. J. Thouless, Phys. Rev. 187, 732 (1969).
gross
D. J. Gross, I. Kanter, and H. Sompolinsky, Phys. Rev. Lett. 55,
304 (1985)
blossey
R. Blossey and J. O. Indekeu, Phys. Rev. E 52, 1223
(1995).
kafri
Y. Kafri, D. Mukamel, and L. Peliti, Phys. Rev. Lett. 85, 4988
(2000).
schwarz
J. M. Schwarz, A. J. Liu, and L. Q. Chayes, Europhys. Lett. 73,
560 (2006).
fisher
C. Toninelli, G. Biroli, and D. S. Fisher, Phys. Rev. Lett. 96,
035702 (2006).
mukamel1
 A. Bar and D. Mukamel,  Phys. Rev. Lett., 112 
015701 (2014).
mukamel2 
A. Bar and D. Mukamel, J. Stat. Mech. (2014) P11001.
dhar
K. E. Bassler, D. Dhar, and R. K. P. Zia , J. Stat. Mech. (2015) P07013.



sandvik1
A. W. Sandvik, Phys. Rev. B 56, 11678 (1997).
sandvik2
A. W. Sandvik, AIP Conf. Proc. 1297, 135 (2010).
sandvik3
O. F. Syljuåsen and A. W. Sandvik, Phys. Rev. E 66, 046701 (2002).
Syljuasen
O. F. Syljuåsen, Phys. Rev. E 67, 046701 (2003).


Scalettar
G. G. Batrouni, R. T. Scalettar, G. T. Zimanyi, and A. P. Kampf, Phys. Rev. Lett. 74, 2527 (1995).











sahinur1
S. Reja, S. Yarlagadda, and P. B. Littlewood, Phys. Rev. B 84, 085127 (2011).
sahinur2
S. Reja, S. Yarlagadda, and P. B. Littlewood, Phys. Rev. B 86, 045116 (2012).
bismuthate3
A. M. Gabovich , A. I. Voitenko and M. Ausloos, Phys. Rep. 367,  583 (2002).
bismuthate2
S. H. Blanton, R. T. Collins, K. H. Kelleher, L. D. Rotter, Z. Schlesinger, D. G. Hinks, and Y. Zheng,
Phys. Rev. B 47, 996 (1993).
quasi_2d_1
R. L. Withers, J.A. Wilson, J. Phys. C 19, 4809 (1986).
quasi_2d_2
J. Merino, R. H. McKenzie, Phys. Rev. Lett. 87, 237002 (2001) 

quasi_1d_1
W. W. Fuller, P. M. Chaikin, N.P. Ong, Phys. Rev. B 24, 1333 (1981).
quasi_1d_2
A. Rusydi, W. Ku, B. Schulz, R. Rauer, I. Mahns, D. Qi, X. Gao, A.T.S. Wee, P. Abbamonte, H. Eisaki, Y. Fujimaki,
S. Uchida, M. Rübhausen, Phys. Rev. Lett. 105, 026402 (2010).

quasi_1d_3
P. Abbamonte, G. Blumberg, A. Rusydi, A. Gozar, P. G. Evans, T. Siegrist, L. Venema, H. Eisaki, E. D. Isaacs,
G. A. Sawatzky, Nature 431, 1078 (2004).
coldatom1
Matthew P. A. Fisher, Peter B. Weichman, G. Grinstein, and Daniel S. Fisher, Phys. Rev. B 40, 546 (1989).
coldatom2
D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, Phys. Rev. Lett. 81, 3108 (1998).
coldatom3
C. Orzel, A. K. Tuchman, M. L. Fenselau, M. Yasuda, and M. A. Kasevich, Science 291, 2386 (2001).
coldatom4
M. Greiner, O. Mandel, T. Rom, A. Altmeyer, A. Widera, T. W. Hnsch and I. Bloch, Nature (London) 415, 39 (2002).
coldatom5
L. M. Duan, E. Demler, and M. D. Lukin, Phys. Rev. Lett. 91, 090402 (2003).

scalettar1 G. G. Batrouni and R. T. Scalettar, Phys. Rev. Lett. 84, 1599
(2000).
scalettar2
F. Hébert, G. G. Batrouni, R. T. Scalettar, G. Schmid,
M. Troyer, and A. Dorneich, Phys. Rev. B 65, 014513 (2001).
boninsegni1
Long Dang, Massimo Boninsegni and Lode Pollet, Phys. Rev. B 78, 132512 (2008).
zoller
B. Capogrosso-Sansone, C. Trefzger, M. Lewenstein, P. Zoller, and G. Pupillo
Phys. Rev. Lett. 104, 125301 (2010).
wessel1
Y.-C. Chen, R. G. Melko, S. Wessel, and Y.-J. Kao, Phys. Rev.
B 77, 014524 (2008).
pinaki
P. Sengupta, L. P. Pryadko, F. Alet, M. Troyer, and G. Schmid
Phys. Rev. Lett. 94, 207202 (2005).
troyer1
G. Schmid and M. Troyer, Phys. Rev. Lett. 93, 067003
(2004).
kar
S. Kar and S. Yarlagadda, Ann. Phys. 375, 322 (2016).
Datta
S. Datta, S. Yarlagadda, Solid State Commun. 150,  2040 (2010).
lv2
Xiao Huo, Yong-Yong Cui, Dali Wang, and Jian-Ping Lv, Phys. Rev. A 95, 023613 (2017).
wessel2
S. Wessel and M. Troyer, Phys. Rev. Lett. 95, 127205 (2005).
kedar1
D. Heidarian and K. Damle, Phys. Rev. Lett. 95, 127206 (2005).
kedar2
A. Sen, P. Dutt, K. Damle, R. Moessner, Phys. Rev. Lett. 100,  147204 (2008).
melko1 R. G. Melko, A. Paramekanti, A. A. Burkov, A. Vishwanath,
D. N. Sheng, and L. Balents, Phys. Rev. Lett. 95, 127207 (2005).
boninsegni2
 M. Boninsegni and N. Prokof'ev, Phys. Rev. Lett. 95, 237204
(2005).
melko2
R. G. Melko, A. Del Maestro, and A. A. Burkov, Phys. Rev. B
74, 214517 (2006).
troyer2
L. Pollet, J. D. Picon, H. P. Büchler, and M. Troyer, Phys. Rev.
Lett. 104, 125302 (2010).

wessel3
S. Wessel, Phys. Rev. B 75, 174301 (2007).
ye
 J. Y. Gan, Y. C. Wen, J. Ye, T. Li, S.-J. Yang, and Y. Yu,
 Phys. Rev. B 75, 214509 (2007).

mishra 
T. Mishra, R. V. Pai, and S. Mukerjee, Phys. Rev. A 89, 013615
(2014).
ghosh
A. Ghosh and S. Yarlagadda, Phys. Rev. B 90, 045140 (2014).
nigel
T. Bilitewski and N. R. Cooper
Phys. Rev. A 94, 023630 (2016).
expt_supersolid
R. Landig, L. Hruby, N. Dogra, M. Landini, R. Mottl, T. Donner,
and T. Esslinger, Nature (London) 532, 476 (2016).

cuprates1
A. Lanzara, P. V. Bogdanov, X. J. Zhou, S. A. Kellar, D. L. Feng, E. D. Lu, T. Yoshida, H. Eisaki, A. Fujimori, K. Kishio, J.-I. Shimoyama, T. Noda, S. Uchida, Z. Hussain, and Z. X. Shen, Nature (London) 412, 510 (2001).
cuprates2
G.-H. Gweon, T. Sasagawa, S. Y. Zhou, J. Graf, H. Takagi, D.-H. Lee, and A. Lanzara, Nature (London) 430, 187 (2004).
manganites1
A. Lanzara, N. L. Saini, M. Brunelli, F. Natali, A. Bianconi, P. G. Radaelli, and S.-W. Cheong, Phys. Rev. Lett. 81, 878 (1998).
manganites2
A. J. Millis, P. B. Littlewood, and B. I. Shraiman, Phys. Rev. Lett. 74, 5144 (1995).
manganites3
F. Massee, S. de Jong, Y. Huang, W. K. Siu, I. Santoso, A. Mans, A. T. Boothroyd, D. Prabhakaran, R. Follath, A. Varykhalov, L. Patthey, M. Shi, J. B. Goedkoop, and M. S. Golden, Nat. Phys. 7, 978 (2011).
fullerides
O. Gunnarsson, Rev. Mod. Phys. 69, 575 (1997).

QMC1
J. E. Hirsch and E. Fradkin, Phys. Rev. B 27, 4302 (1983).
QMC2
J. E. Hirsch, Phys. Rev. B 31, 6022 (1985).
QMC3
E. Berger, P. Valášek, and W. von der Linden, Phys. Rev. B 52, 4806 (1995).
QMC4
Z. B. Huang, W. Hanke, E. Arrigoni, and D. J. Scalapino, Phys. Rev. B 68, 220507(R) (2003).
QMC5
R. P. Hardikar and R. T. Clay, Phys. Rev. B 75, 245103 (2007).
QMC6
A. Macridin, G. A. Sawatzky, and M. Jarrell, Phys. Rev. B 69, 245111 (2004).
exact1
A. Dobry, A. Greco, J. Lorenzana, and J. Riera, Phys. Rev. B 49, 505 (1994).
exact2
A. Dobry, A. Greco, J. Lorenzana, J. Riera, and H. T. Diep, Europhys. Lett. 27, 617 (1994).
exact3
B. Bäuml, G. Wellein, and H. Fehske, Phys. Rev. B 58, 3663 (1998).
DMRG1
M. Tezuka, R. Arita, and H. Aoki, Phys. Rev. B 76, 155114 (2007).
DMRG2
Shigetoshi Sota and Takami Tohyama, Phys. Rev. B 82, 195130 (2010).
DMFT1
J. K. Freericks and M. Jarrell, Phys. Rev. Lett. 75, 2570 (1995).
DMFT2
M. Capone, G. Sangiovanni, C. Castellani, C. Di Castro, and M. Grilli, Phys. Rev. Lett. 92, 106401 (2004).
DMFT3
W. Koller, D. Meyer, Y. Ono, and A. C. Hewson, Europhys. Lett. 66, 559 (2004).
DMFT4
W. Koller, D. Meyer, and A. C. Hewson, Phys. Rev. B 70, 155103 (2004).
DMFT5
G. S. Jeon, T.-H. Park, J. H. Han, H. C. Lee, and H.-Y. Choi, Phys. Rev. B 70, 125114 (2004).
DMFT6
G. Sangiovanni, M. Capone, C. Castellani, and M. Grilli, Phys. Rev. Lett. 94, 026401 (2005).
DMFT7
G. Sangiovanni, M. Capone, and C. Castellani, Phys. Rev. B 73, 165123 (2006).
DMFT8
J. Bauer and A. C. Hewson, Phys. Rev. B 81, 235113 (2010).
DMFT9
Johannes Bauer and Giorgio Sangiovanni, Phys. Rev. B 82, 184535 (2010).
slaveboson1
M. Grilli and C. Castellani, Phys. Rev. B 50, 16880 (1994).
slaveboson2
J. Keller, C. E. Leal, and F. Forsthofer, Physica B 206-207, 739 (1995).
slaveboson3
E. Koch and R. Zeyher, Phys. Rev. B 70, 094510 (2004).
slaveboson4
U. Trapper, H. Fehske, M. Deeg, and H. Buttner, Z. Phys. B: Condens. Matter 93, 465 (1994).
slaveboson5
C. A. Perroni, V. Cataudella, G. De Filippis, and V. Marigliano Ramaglia, Phys. Rev. B 71, 113107 (2005).
variational1
Y. Takada and A. Chatterjee, Phys. Rev. B 67, 081102 (2003).
variational2
H. Fehske, D. Ihle, J. Loos, U. Trapper, and H. Buttner, Z. Phys. B: Condens. Matter 94, 91 (1994).
largeN
R. Zeyher and M. L. Kuli ́c, Phys. Rev. B 53, 2850 (1996).
Gutzwiller1
A. Di Ciolo, J. Lorenzana, M. Grilli, and G. Seibold, Phys. Rev. B 79, 085101 (2009).
Gutzwiller2
P. Barone, R. Raimondi, M. Capone, C. Castellani, and M. Fabrizio, Phys. Rev. B 77, 235115 (2008).
cluster
Alexandre Payeur and David S ́en ́echal, Phys. Rev. B 83, 033104 (2011).


sandvik_paper
A. W. Sandvik, Phys. Rev. B 56, 11678 (1997).
sandvik_review
A. W. Sandvik, AIP Conf. Proc. 1297, 135 (2010).
dir_loop1
O.F. Syljuasen, A.W. Sandvik, Phys. Rev. E 66, 046701 (2002).
dir_loop2
O.F. Syljuasen, Phys. Rev. E 67, 046701 (2003).
Batrouni
G.G. Batrouni, R.T. Scalettar, G.T. Zimanyi, A.P. Kampf, Phys. Rev. Lett. 74, 2527 (1995).
Damle
A. Sen, K. Damle, and T. Senthil, Phys. Rev. B 76, 235107 (2007).
Wessel1
S. Wessel, Phys. Rev. B 78, 075112 (2008).
Wessel2
S. Wessel, Phys. Rev. B 86, 140501(R) (2012).
