Frustration induced incommensurate solids in the extended Bose-Hubbard model
FFrustration induced incommensurate solids in the extended Bose-Hubbard model
Kwai-Kong Ng ∗ and Min-Fong Yang † Department of Applied Physics, Tunghai University, Taichung 40704, Taiwan (R.O.C.) (Dated: March 1, 2021)The extended Bose-Hubbard model with nearest-neighbor and next-nearest-neighbor repulsiveinteractions on a square lattice is investigated by using the quantum Monte Carlo method. We findthat, for the cases of weak next-nearest-neighbor interactions and small hoppings, incommensuratesolids of fractional densities varying from 1/4 to 1/2 can be stabilized in the thermodynamic limit.We further show that the continuous changes of ordering wave vectors from ( π, π/
2) [or ( π/ , π )]to ( π, π ) within the incommensurate solid phase can be understood by the mechanism of domainwall formation. The related ground-state phase diagram and thermal phase transitions are alsodiscussed. I. INTRODUCTION
Polar molecules trapped in optical lattices provide anunique opportunity to study the dipole-dipole interac-tions in real experimental setups [1–4]. One primary in-vestigation is to identify possible quantum phases thatmay uniquely arise from these long-range interactions butcannot be observed in systems with only short-range in-teractions. Recent numerical studies [5–8] of hard-corebosons with infinite-range interactions on two dimen-sional (2D) square lattices have presented evidences ofMott insulators with checkerboard, stripe, and star or-dering at densities ρ = 1/2, 1/3, and 1/4, respectively.Furthermore, supersolids around the Mott lobes with ρ =1/2 and 1/4 are also found by doping the solids with par-ticles or vacancies. Besides, devil’s staircase (that is, asequence of commensurate phases separated by series ofjumps) is observed for finite-size systems, which signalsthe presence of the incommensurate phases in the ther-modynamic limit. A question is: are all these interestingphases uniquely stabilized by long-range interactions?Actually, the Mott insulators and the supersolids canbe realized in short-range models as well. In previousnumerical investigations on the extended Bose-Hubbardmodel that includes only the nearest-neighbor (nn) andthe next-nearest-neighbor (nnn) interactions, Mott in-sulators at 1/2 and 1/4 fillings and the associated su-persolids have indeed been observed [9–14]. The modelHamiltonian for the hard-core bosons on a 2D square lat-tice is, H = − t (cid:88) (cid:104) i,j (cid:105) ( b † i b j + h.c. ) + V (cid:88) (cid:104) i,j (cid:105) n i n j + V (cid:48) (cid:88) (cid:104) i,j (cid:105) (cid:48) n i n j − µ (cid:88) i n i . (1)Here b i ( b † i ) is the annihilation (creation) bosonic opera-tor on site i , t the hopping integral, n i the particle num-ber on site i , and V ( V (cid:48) ) the nn (nnn) repulsive interac-tion. µ denotes the chemical potential which controls the ∗ [email protected] † [email protected] particle density in the grand canonical ensemble. For thecases of strong V (cid:48) ( V (cid:48) > V / t and weak nnn interac-tions V (cid:48) ( V (cid:48) < V / V (cid:48) cases [9–13], there exists no supersolid phases in thepresent case of weak V (cid:48) . The incommensurate solids arefound to be characterized by a continuous change in or-dering wave vectors as model parameters are varied. Theformation of such a phase can be explained by the prolif-eration of domain walls from half-filled or quarter-filledsolids. Following the analysis in Ref. [16] in the context ofthe anisotropic triangular lattice, analytical calculationson the domain-wall excitation energies are made for thepresent system on a square lattice. Our analytical pre-dictions on the phase boundaries agree very well withthe QMC numerical results. Furthermore, the quantum a r X i v : . [ c ond - m a t . o t h e r] F e b phase transitions from the incommensurate solids to thehalf-filled or quarter-filled solids are found to be contin-uous in the thermodynamic limit, while the transitionto the superfluid phase is first-order. Upon increasingtemperatures, the incommensurate solids will melt intonormal fluids via a first-order thermal phase transition,which are demonstrated by hysteresis of the structurefactors and double peaks in histogram.This paper is organized as follows. The ground-statephase diagram at fixed V and V (cid:48) are identified in Sec. II.The general features of the incommensurate solids areexplained in Sec. III by the domain wall formation. Theenergetic perspective of the domain wall are elucidated inSec. IV. Analytic expressions of phase boundaries of theincommensurate solids and of the domain wall density arethen derived. More quantum and thermal phase transi-tions out of the incommensurate solids are discussed inSec. V. We summarize our work in Sec. VI. II. QUANTUM PHASE DIAGRAM
In this work we utilized the well established stochasticseries expansion algorithm [17, 18] in the QMC calcula-tions as the Hamiltonian H in Eq. (1) is sign-problemfree. Unless mentioned otherwise, we set V (cid:48) = 1 asthe energy scale. The inverse temperature is set to be β = 1 / (2 L ), where L is the linear size of lattice ( L = 12,24, 36, and 48 are used). To identify various quantumphases, besides the particle density ρ , two more quan-tities are measured in our simulations., The superfluid-ity ρ s is calculated to signal off-diagonal U(1) symmetrybreaking, which is defined by the fluctuations of windingnumbers W x and W y , ρ s = 12 β ( (cid:104) W x (cid:105) + (cid:104) W y (cid:105) ) . (2)On the other hand, the translational broken symmetry ischaracterized by the structure factor S ( Q ) = 1 L (cid:88) i,j (cid:104) n i n j (cid:105) e i Q · r ij , (3)where Q denotes the wave vector and r ij the displace-ment between two sites i and j .In the parameter regime we studied in our model, theoff-diagonal and the diagonal orderings do not coexist,that means no supersolid is found. In Fig. 1, the ground-state phase diagram µ vs t with V = 4 contains a super-fluid (SF) and two solid states of definite ordering wavevectors Q signaling the peak in structure factor. Theirphase boundaries are determined by the abrupt changesof ρ s or S ( Q ) (see Fig. 2) and then by extrapolating tothe thermodynamic limit with L → ∞ . The half-filledphase is a checkerboard solid (CBS) with Q = ( π, π ),instead of a striped solid, because the nn repulsion isdominant ( V > V (cid:48) ) so that nn occupation is avoided.Reducing the chemical potential such that ρ = 1 /
4, the
FIG. 1. Ground-state phase diagram with V = 4 obtainedfrom QMC calculations after extrapolating to thermodynamiclimit. The dotted lines are the phase boundaries predicted bythe domain wall analysis. The insets show the CBS and theQFS lattice structures. Here V (cid:48) = 1 as the energy unit. ground state is then a quarter-filled solid (QFS) of starordering [12, 13] with Q = ( π/ , π ) or ( π, π/ t = 0, the degeneracy is lifted for finite t . Thepreferred QFS ordering is shown in the inset of Fig. 1 asit prevents nn occupation up to second order of hopping t to reduce the cost of the nn potential energy. Interest-ingly, for finite but small t , the phase in between the CBSand the QFS is neither a superfluid nor supersolid, butan incommensurate solid (ICS) with ordering wave vec-tors Q continuously varying from ( π/ , π ) [or ( π, π/ π, π ) in the thermodynamic limit. As explained inthe next section, we found that ICS can be understoodby inserting linear domain walls into the solid phases. Byincreasing hopping t , all solid phases eventually melt intoa superfluid phase directly without passing an intermedi-ate supersolid phase. Further discussions on these phasetransitions will be provided in Section V. III. INCOMMENSURATE SOLIDS
For a given hopping t = 0 .
4, the particle density ρ isshown in Fig. 2 for L = 36. A series of density plateaus isobserved from ρ = 1 / m/L , with m being an integer from L/ L/
2. Thisis very similar to the staircase found in the model withinfinite-range interactions [3, 5]. In that case, although asmaller lattice size ( L = 12) is used, much more plateausof rational fillings are found between ρ =1/4 and 1/2.The reason for this difference is that the infinite-rangenature of the interactions allows more possible configu-rations of particle ordering to be stabilized, which, un-fortunately, also makes the physics more difficult to be ρ µ
36 9 36 10 36 11 36 12 36 13 36 14 36 1536 16 36 17 36 18SF QFS ICS CBS
FIG. 2. Order parameters as functions of µ with t = 0 . V = 4. Upper panel: the particle density ρ as a function of µ .The fractional numbers are the corresponding densities of theplateaus. Solid (Open) circles represent the data for L = 36( L = 24). Lower panel: superfluidity ρ s (filled squares) andstructure factors S ( Q ) (open circles) as functions of µ . Thecorresponding ordering wave vectors of the data are, fromleft to right, Q =(2 mπ/L , π ) with m = L/
4, ( L/
4) + 1, ..., L/
2. For simplicity, the structure factors for the orderingwave vectors Q =( π , 2 mπ/L ) are not shown. Note that thedata of S ( π, π ) is reduced by 5 times. analyzed. In contrast, for the present system with short-range interactions, we can understand the origin of thedensity plateaus by a simple picture of linear fluctuatingdomain walls. While both domain walls along the x andthe y directions are observed in our simulations, for sim-plicity, only the domain walls along the y direction willbe discussed hereafter.Under doping the half-filled CBS with holes, the po-tential energy gain is 4 V (cid:48) for each added isolated hole.Nevertheless, holes aligned together to form a domainwall [Fig. 3 (b)] can further gain kinetic energy throughparticle hoppings, which will be explained in details inthe next section when we discuss the domain wall dy-namics. For convenience, we define L x ( L y ) to be thelattice size in the x ( y ) direction. We note that evennumbers of domain walls are required to maintain theperiodic boundary condition imposed in our system. Be-cause L y / y direction out of the CBS, thetotal number of removed particle becomes nL y when 2 n domain walls appear in the lattice for n = 0, 1, 2, ..., L x /
4. Consequently, when there are L x / / ρ = 1 / − / / n domain walls,the particle densities become ρ = ( L x − n ) / (2 L x ), whichexactly correspond to those of the plateaus found in Fig. 2for both L = 24 and 36 (note that L = L x = L y in all our simulations). Our analysis clearly shows that the seriesof plateaus are formed by adding pairs of linear domainwalls from the CBS.Furthermore, the observed shift in the ordering wavevector Q in Fig. 2 implies the presence of linear domainwalls as well. We note that the insertion of domain wallsat equal distance will not immediately destroy the ordercompletely. Instead, due to producing long-period super-structures in the x direction, the ordering wave vectorwill be shifted by (2 π/L x ,
0) for each added pair of lineardomain walls along the y direction. Hence the plateaustates in ICS with 2 n domain walls along the y direc-tion can be characterized by the ordering wave vectors Q = ( π − nπ/L x , π ) or Q = (2 mπ/L x , π ) with m beingan integer from L x / L x /
2. (The same discussions ap-ply as well for the states with 2 n domain walls along the x direction.) This conclusion is consistent to the numericalresults shown in Fig 2.Another implication from the above analysis is thatas L → ∞ , the number of plateaus will become infinitysuch that ρ changes continuously from 1/4 to 1/2. Thedensity jump at the phase boundary between QFS (CBS)and ICS is then expected to be reduced to zero and thecorresponding commensurate-incommensurate transitionshould thus be continuous, rather than of first order. Thelocation of these phase boundaries in the ground-statephase diagram for small t can be determined analyticallyby considering the domain wall dynamics, as discussedin the next section. IV. DOMAIN WALL DYNAMICS
Here we follow the analysis in Ref. [16] for incommen-surate supersolids on anisotropic triangular lattices todiscuss the domain wall dynamics in the present model.Started from the half-filled CBS state, one can split theCBS into half by removing a column of bosons (i.e., re-moving L y / E p /L y =( µ − V (cid:48) ) /
2. After shifting half of the lattice upwards, adomain wall is created without further cost in potentialenergy. To be specific, we define the domain walls byconnecting the center of unoccupied bonds, shown as thered solid zigzag line in Fig. 3(b). Note that the domainwall can fluctuate to gain the kinetic energy, because thebosons around the shifted boundary are now free to movesideways. As illustrated in Fig. 3(b), after the hoppingof bosons, parts of the domain wall fluctuate in the op-posite direction of the particle movement (shown as thered dotted lines). To estimate the gain in kinetic energy,one can map the one-dimensional (1D) domain wall ontoa 1D spin-1/2 XY chain [16, 19, 20]. That is, by tracingthe zigzag chain along y direction, a up (down) spin isassigned if the x -coordinate of a given unoccupied bondis increased (decreased) by one unit. The kinetic energyof a domain wall thus corresponds to the ground-stateenergy of the 1D spin chain and has the value (per unit FIG. 3. Illustrations of domain wall motion. (a) Removinga column of bosons. (b) Shifting the right half of the latticecosts no extra potential energy but gains kinetic energy dueto the fluctuation of the domain wall (red line). (c) Multipledomain walls. (d) QFS state with a domain wall density ρ D =1 / length) ∆ E k /L y = − t/π [19]. Therefore, the total exci-tation energy (per unit length) of a domain wall becomes∆ E/L y = ∆ E p /L y + ∆ E k /L y = ( µ − V (cid:48) ) / − t/π .When its excitation energy is reduced to zero by decreas-ing µ , the phase transition out of the CBS induced byproliferation of domain walls will occur. This gives thecritical µ c for the CBS-ICS phase transition, i.e. µ c = 4 V (cid:48) + 4 tπ . (4)This quantum phase boundary is plotted in Fig. 1 (theupper dotted line), which is in good agreement with theQMC results.For the transition from QFS to ICS, the same argu-ment applies. Doping the QFS with a column of bosons,it actually generates two domain walls and the total costin potential energy per unit length equals to (4 V (cid:48) − µ ) / − t/π ). Therefore, the critical µ c for the QFS-ICSphase transition is given by µ c = 4 V (cid:48) − tπ , (5)which is again in good agreement with the numerical re-sults. Deviations in analytical predictions of µ c and µ c are expected for large t , since higher order correctionsand SF fluctuations become important then.Following the analysis in Ref. [16], we also determinethe domain wall density ρ D as a function of µ , which is ρ D µ ρ D t=0.4t=0.2 ρ D ρ D f( ρ D )f( ρ D ) FIG. 4. Domain wall density ρ D as a function of µ for t = 0 . t = 0 . f ( ρ D ) to a power-lawfunction Aρ αD for each case. Here V = 4 and L = 24 aretaken. defined as: ρ D = 1 L x L y (cid:88) x,y ¯ n x,y ¯ n x +1 ,y , (6)where ¯ n x,y = 1 − n x,y is the hole number at site ( x, y ).Here the summation of ¯ n x,y ¯ n x +1 ,y over x coordinatescounts the number of hole-hole bonds in the x direction,and then gives the number of linear domain walls alongthe y direction. The QMC results of ρ D are shown inFig. 4.Clearly, the domain walls will interact with each other.In combination with the discussions in the first paragraphof this section for isolated domain walls, the total energyfor interacting domain walls at a finite density ρ D willtake the following form: E ( ρ D ) = L x L y V (cid:48) ρ D (cid:20) − µ V (cid:48) − tπV (cid:48) + f ( ρ D ) (cid:21) . (7)Here the last term accounts for an effective repulsive in-teraction energy between two neighboring domain walls.This yet unknown function f ( ρ D ) can be determined byfitting to the numerical results.We note that the transition from the state with 2 M − M domain walls occurs when E [(2 M − /L x ] = E (2 M/L x ). This condition leads to arecursive relation of f (2 M/L x ): µ M = µ c + 2 V (cid:48) (cid:20) ( M − f ( 2 M − L x ) − M f ( 2 ML x ) (cid:21) , (8) ρ s S( t ρ s S( (a)(b) µ=4.0µ=3.2 FIG. 5. Superfluid density ρ s (open circles) and structurefactor S ( Q ) (open squares) as functions of t for (a) µ = 4 . µ = 3 .
2. Here V = 4 and L = 24. Both quantitiesshow the signatures of first-order phase transition from theICS to the SF states. which can be solved to have f ( 2 ML x ) = 1 V (cid:48) (cid:34) µ c − M (cid:88) i =1 µ i M (cid:35) . (9)As a result, the simulated values of µ M can be usedto evaluate the discrete values of f (2 M/L x ). As shownin the inset of Fig. 4, a power-law function Aρ αD can befitted to f ( ρ D ) with the exponent α = 3 . t = 0 . ρ D = (cid:20) µ c − µ V (cid:48) A ( α + 1) (cid:21) /α , (10)which is found to be in good agreement with the simu-lated data in Fig. 4. Our results thus justify the applica-tion of the domain wall dynamics in the present system.Interestingly, we repeated the same calculations on do-main wall density for t = 0 . α = 3 . V. INCOMMENSURATE-SUPERFLUID ANDTHERMAL PHASE TRANSITIONS
In our parameter regime, the quantum phase transi-tions from all solid states to the SF phase are observed tobe discontinuous and there appears no immediate super-solid phase. These discontinuous transitions are demon-strated by the abrupt changes of the superfluid density ρ s and the structure factor S ( Q ) at the phase boundaries.Here we focus on the transitions out of the ICS states. ρ ρ s ρ (a)(b) FIG. 6. Boson density ρ (circles) and superfluid density ρ s (squares) as functions of V (cid:48) with V = 4 for (a) µ = 4 . µ = 4 .
6. Here L = 36 and V = 4. In Fig. 5, we present our data for µ = 4 . t with ρ = 1 / t , wefind that both ρ s and S ( Q ) change abruptly across thephase boundaries. We note that the solid ordering isdestroyed at the same critical t where the superfluidityemerges. Hence no signs of coexistence of both order pa-rameters is found. The absence of supersolid phase couldbe accounted by the same argument for the CBS and theQFS that phase separation is more favorable than thesupersolid phase [9].In the previous discussions, we fix the nnn interaction V (cid:48) to be unity. We now turn to the effect from varying V (cid:48) .When V (cid:48) is increased, the enhanced nnn frustrations willreduce the cost in potential energy for the domain wallformation (see section IV). Therefore, the boson densityof the system is expected to follow a series of densitydrops as more domain walls of holes can be generated forlarger V (cid:48) . Our numerical results for V = 4 at µ = 4 . µ = 4 .
0, the system moves from the 1/3-filling stateall the way to the quarter-filling state as V (cid:48) increases.However, at µ = 4 .
6, when V (cid:48) > .
24, the frustrationbecomes too strong such that the diagonal long-rangeorder cannot be sustained and the system melts into asuperfluid phase instead. An abrupt jump in superfluiddensity is observed that signals a direct first-order phasetransition from the ICS to the SF and there exists nointermediate supersolid phase.We now consider the thermal phase transitions out ofthe ICS and the QFS phases. It is expected that thebroken translation symmetry in these two phases will berestored by thermal fluctuations at high temperatures.This symmetry restoration can be detected by the dis-appearance of the corresponding structure factors. Care-ful studies of the order-disorder transitions indicate thatthese transitions are of first order. In Fig. 7, we mea- S ( π / , π ) T S ( π / , π ) increasing Tdecreasing T0 0.010.020.030.040.050.06010002000 T=0.208T=0.1389 µ=2.0µ=4.0 (a)(b) FIG. 7. Hysteresis of the structure factors S( Q ) for the casesof (a) µ = 4 . µ = 2 .
0, which correspond the 1/3-filling plateau state and the 1/4-filling QFS state respectively.Here V = 4 and L = 24. The insets show the double-peakfeature of S( Q ) around the transitions. sure the structure factors S (2 π/ , π ) and S ( π/ , π ) ofthe ρ = 1 / ρ = 1 / µ in-creases from the QFS to the ICS, the thermal transitionbecomes a more weakly first-order one. This is consis-tent to the fact that the thermal transition of half-filledCBS is known to be continuous, so that approaching theCBS from the QFS via the ICS, the thermal transitionsare expected to evolve from strongly first-order to weaklyfirst-order and then to second-order. Since the orderingwave vector Q of the maximal structure factor varies con-tinuously in the ICS phase into the CBS, it is natural toexpect a smooth evolution of transition order instead ofan abrupt change as Q approaches to ( π , π ). VI. CONCLUSION
We have shown that the ICS states can emerge in abosonic model on a square lattice under short-range frus-trated interactions. This observation has been overlookedin previous studies of the extended Bose-Hubbard model.The ICS phase appears in between the half-filled and thequarter-filled solids and is characterized by a series ofdensity plateaus with fractional values of densities. Theshort-range character of the interactions in our modelsimplifies the analysis and allows us to address funda-mental questions about the nature of the plateau states.By following the analysis in Ref. [16], we show that thefractional particle densities and the ordering wave vec-tors of the plateau states can be explained by the domainwall formation. This indicates that the incommensuratephase originates from proliferation of domain walls. Fur-thermore, the measured domain wall densities ρ D agreewell with the predicted values from the theory of theinteracting domain walls. We find that the interactionterm behaves as a power-law function of ρ D with an ex-ponent likely to be independent of the system parapme-ters. In the thermodynamic limit, we expect the order-ing wave vector Q characterizing different broken sym-metries changes continuously from Q = ( π/ , π ) for thequarter-filled star order to Q = ( π, π ) for the checker-board order. In addition, the widths of the intermediateplateau states become diminished. There is no super-fluidity found in the ICS phase, differs from the case ofanisotropic triangular lattice where an incommensuratesupersolid is observed [16]. It is not surprising since, incontrast to the case of the square lattice, bosons havemore degrees of freedom to hop around on the triangularlattice. Based on this observation, incommensurate su-persolids may come out in our model if the nnn hoppingis turned on. It is also interesting to consider the effect ofstrong frustration (large V (cid:48) ) on the ICS phase. We showthat, by increasing V (cid:48) , either the ICS state is melted intoa SF or the system is transformed to be a QFS. Never-theless, it is known that for even stronger V (cid:48) anothertype of star order becomes favorable at quarter filling.Whether similar ICS phases can still appear around suchquarter-filled solids is not clear without extensive studieson a wide range of parameter regime. Furthermore, onemay wonder if the present mechanism of domain wall for-mation can be applied to explain the staircase observedin the infinite-range models. This is worth to investigatein a further study. ACKNOWLEDGMENTS
The numerical computations are performed in the Cen-ter for High Performance Computing of the THU. Thiswork is financial supported by the Ministry of Scienceand Technology (R.O.C.) Grant No. 108-2112-M-029 -005and 109-2112-M-029-006. M.F. Yang is also supported bythe Ministry of Science and Technology (R.O.C.) GrantNo. 109-2112-M-029-005. [1] T. Lahaye, C. Menotti, L. Santos, M. Lewenstein, andT. Pfau, The physics of dipolar bosonic quantum gases,Rep. Prog. Phys , 126401 (2009).[2] K.-K. Ni, S. Ospelkaus, M. H. G. de Miranda, A. Pe’er,B. Neyenhuis, J. J. Zirbel, S. Kotochigova, P. S. Julienne,D. S. Jin, and J. Ye, A high phase-space-density gas ofpolar molecules, Science , 231 (2008).[3] S. Ospelkaus, K.-K. N. A. Pe’er, J. J. Zirbel, B. Neyen-huis, S. Kotochigova, P. S. Julienne, J. Ye, and D. S.Jin, Efficient state transfer in an ultracold dense gas ofheteronuclear molecules, Nature Phys. , 622 (2008).[4] R. Landig, L. Hruby, N. Dogra, M. Landini, R. Mottl,T. Donner, and T. Esslinger, Quantum phases from com-peting short- and long-range interactions in an opticallattice, Nature , 476 (2016).[5] B. Capogrosso-Sansone, C. Trefzger, M. Lewenstein,P. Zoller, and G. Pupillo, Quantum phases of cold po-lar molecules in 2d optical lattices, Phys. Rev. Lett. ,125301 (2010).[6] T. Ohgoe, T. Suzuki, and N. Kawashima, Quantumphases of hard-core bosons on two-dimensional latticeswith anisotropic dipole-dipole interaction, Phys. Rev. A , 063635 (2012).[7] D. Yamamoto, A. Masaki, and I. Danshita, Quantumphases of hardcore bosons with long-range interactionson a square lattice, Phys. Rev. B , 054516 (2012).[8] P. Rosson, M. Kiffner, J. Mur-Petit, and D. Jaksch, Char-acterizing the phase diagram of finite-size dipolar bose-hubbard systems, Phys. Rev. A , 013616 (2020).[9] G. G. Batrouni and R. T. Scalettar, Phase separation insupersolids, Phys. Rev. Lett. , 1599 (2000).[10] F. H´ebert, G. G. Batrouni, R. T. Scalettar, G. Schmid,M. Troyer, and A. Dorneich, Quantum phase transitionsin the two-dimensional hardcore boson model, Phys. Rev. B , 014513 (2001).[11] K.-K. Ng, Y. C. Chen, and Y. C. Tzeng, Quarter-filledsupersolid and solid phases in the extended bose–hubbardmodel, Journal of Physics: Condensed Matter , 185601(2010).[12] K.-K. Ng and Y.-C. Chen, Supersolid phases in thebosonic extended hubbard model, Phys. Rev. B ,052506 (2008).[13] L. Dang, M. Boninsegni, and L. Pollet, Vacancy super-solid of hard-core bosons on the square lattice, Phys. Rev.B , 132512 (2008).[14] K.-K. Ng, Thermal phase transitions of supersolids in theextended bose-hubbard model, Phys. Rev. B , 184505(2010).[15] S. V. Isakov, H.-C. Chien, J.-J. Wu, Y.-C. Chen, C.-H. Chung, K. Sengupta, and Y. B. Kim, Commensuratelock-in and incommensurate supersolid phases of hard-core bosons on anisotropic triangular lattices, EPL (Eu-rophysics Letters) , 36002 (2009).[16] X.-F. Zhang, S. Hu, A. Pelster, and S. Eggert, Quantumdomain walls induce incommensurate supersolid phase onthe anisotropic triangular lattice, Phys. Rev. Lett. ,193201 (2016).[17] A. W. Sandvik, Stochastic series expansion method withoperator-loop update, Phys. Rev. B , R14157 (1999).[18] O. F. Sylju˚asen and A. W. Sandvik, Quantum montecarlo with directed loops, Phys. Rev. E , 046701(2002).[19] F. Mila, Exact result on the mott transition in a two-dimensional model of strongly correlated electrons, Phys.Rev. B , 14047 (1994).[20] N. G. Zhang and C. L. Henley, Stripes and holes in atwo-dimensional model of spinless fermions or hardcorebosons, Phys. Rev. B68