Frustration-Induced Supersolid Phases of Extended Bose-Hubbard Model in the Hard-Core Limit
FFrustration-induced supersolid phases of extended Bose-Hubbard model in thehard-core limit
Wei-Lin Tu , Huan-Kuang Wu , and Takafumi Suzuki Institute for Solid State Physics, University of Tokyo, Kashiwa, Chiba 277-8581, Japan Department of Physics, Condensed Matter Theory Center and Joint Quantum Institute,University of Maryland, College Park, MD 20742, USA Graduate School of Engineering, University of Hyogo, Hyogo, Himeji 670-2280, Japan (Dated: March 10, 2020)We investigate exotic supersolid phases in the extended Bose-Hubbard model with infinite pro-jected entangled-pair state, numerical exact diagonalization, and mean-field theory. We demonstratethat many different supersolid phases can be generated by changing signs of hopping terms, and theinteractions along with the frustration of hopping terms are important to stabilize those supersolidstates. We argue the effect of frustration introduced by the competition of hopping terms in thesupersolid phases from the mean-field point of view. This helps to give a clearer picture of thebackground mechanism for underlying superfluid/supersolid states to be formed. With this knowl-edge, we predict and realize the d -wave superfluid, which shares the same pairing symmetry withhigh- T c materials, and its extended phases. We believe that our results contribute to preliminaryunderstanding for desired target phases in the real-world experimental systems. I. INTRODUCTION
Supersolid (SS) phase, after its first proposition byPenrose and Onsager [1], has attracted much attentionfrom both experimental and theoretical aspects [2]. Thenature of the SS state is characterized by the coexistenceof crystal order and superfluidity, namely the coexistenceof diagonal and off-diagonal long-range order. The SSphase is formed by adding dopants at a commensuratefilling, where perfect crystal exists. These dopants thencondensate and contribute to the superflow [3–5] whilekeeping the crystal order. This scenario is called the“defect-condensation” [6].Although great efforts have been paid searching forthe SS phase experimentally, the observation of its ex-otic property has been never reported yet in natural ma-terials. In 2004, Kim and Chan have reported the co-existence of solidity and superfluidity in helium-4 [7].They have observed a sudden drop of the non-classicalrotational moment of inertia in the torsional oscillatorexperiment. After their work, many experiments havebeen performed in order to argue more of this fascinatingphenomenon [8–10]. However, it has later been demon-strated that the supersolidity is lost when the effect ofelastic properties of helium is removed by careful exper-iments [11, 12].With the observation of supersolidity from heliumseeming to be unlikely, our attention turns to anotherplatform, the ultracold atomic gases [13–16]. Many pos-sible scenarios to realize the SS state have been proposedfor the ultracold atomic gases by adding additional in-teractions. Among all, the dipolar gases [17, 18] seemto become a more reliable platform for supersolidity. Inquite recent experiments, strong evidences of supersolid-ity have been found in dipolar atoms made of erbium anddysprosium [19–24].With the advanced techniques introduced by ultracoldexperiments, theorists have also been urged to investigate this issue more deeply. One of the most commonly usedHamiltonians for tackling this issue is the extended Bose-Hubbard (EBH) model, which can be experimentally re-alized by the dipolar gases [17, 18, 25]. In earlier works,it has been revealed that by the next-nearest-neighborinteraction, the SS phase can be formed and appears viaa second-order phase transition from the superfluid (SF)phase on a square lattice [26, 27]. When the furtherlong-range interaction, namely the dipole-dipole interac-tion in the cold atom experiment, is included into theEBH model, the SS phase still survives [28]. These re-sults have highlighted the importance of a longer-rangeinteraction for the stabilization of the SS state.In the effect of the long-range interaction, the frustra-tion between interactions, which seems to be essentialto realize the quantum spin liquid in many modern con-densed matter systems[29, 30], can be crucial in formingthe SS phase. A direct mapping from hard-core bosons toS=1/2 spins enables us to discuss ordering phases fromone side complementary to the other side [31]. On theother hand, the abundant underlying physics of fermionicHubbard models concerning high- T c superconductivity[32–36] also implies that there could be some interest-ing features for the bosonic side.Early efforts upon the EBH model for searching theSS state were made by mainly the quantum Monte Carlo(QMC) method and contributed to fruitful results [6, 37–39]. In these studies, the systems with the frustrated hop-ping interactions were ignored owing to the negative signproblem in QMC calculations. However, recent worksconcerning the same model have shown that with thefrustration induced by a negative next-nearest-neighborhopping t (cid:48) , a peculiar SS phase, named after half super-solid (HSS), can be generated [40, 41]. To further inves-tigate the details of HSS, one might need other numericalapproaches.In this paper, our central method for solving the EBHmodel in the hard-core limit is through one of the ten- a r X i v : . [ c ond - m a t . o t h e r] M a r sor network (TN) ansatz, known as infinite projectedentangled-pair state (iPEPS) [42]. For this TN ansatzin a square lattice, quantum state on each lattice site isrepresented by a rank-5 tensor with one physical indexof dimension d and four auxiliary indices of dimension D in a two-dimensional plane. Because of the hard-corelimit, each site can have two different states and therefore d = 2, while D can be chosen freely but it helps enhancethe accuracy with larger numbers. To attain the ther-modynamic limit, here we use the corner-transfer-matrixalgorithm [33, 43, 44]. After the iterative convergence,we obtain a series of environment tensors surroundingthe unit cell, which we choose to be equal to the size of2 ×
2. For the corner and edge tensors, composing theenvironment tensors, their bond dimension, χ , is chosento be large enough ( χ ( D ) > D ) for minimizing the errorcaused by the usage of finite χ .To obtain a proper ground state of the target Hamilto-nian in iPEPS, we can achieve it by projecting the initialstates into the imaginary-time evolution or through thevariational update [45, 46]. In the following calculation,we apply the imaginary-time evolution in which the en-vironment tensors are not included while all tensors ofthe unit cell being updated. Thus, this imaginary-timeevolution method is also called the simple update [47].Compared with another more accurate but computation-ally expensive update algorithm, namely the full-updatemethod [42, 48], the simple-update method would loseits accuracy, once the system of interest is highly cor-related. Nevertheless, the simple-update method withsmall D has provided the accurate phase diagrams witha dimer model [49], which is in fact more entangled thanour model owing to the longer-range interdimer couplingterm. Therefore, we apply simple-update iPEPS as one ofour numerical methodologies. Besides iPEPS, it has beenknown that the mean-field theory works well for the EBHmodel [41]. Thus, we employ the mean-field theory pro-posed by Matsubara and Matsuda [50] to compare withthe results by iPEPS. In addition to the above two meth-ods, we execute exact diagonalization (ED) numerically.The combination of ED and the finite-size scaling alsoprovides further evidences of the phases we have found.This paper is organized in the following way. In SectionII, we will present our results. We firstly define the orderparameters for distinguishing different phases and applythe iPEPS and mean-field approaches for constructingphase diagrams of given hopping and interacting ampli-tudes in Section II A and B. We then re-examine the exis-tence of these phases with ED in Section II C. In SectionII D, we analyze the causes of found phases in SectionII B through a mean-field interpretation and constructgeneral rules in pursuit of desired SF/SS states. We fol-low these rules and predict a d -wave SF along with itsrelated phases. Our conclusion is carried in Section IIIand the justification of finite- D iPEPS is included in theAppendix. FIG. 1: Configurations for different orders: (a)-(c) arefor real-space particle density modulation. Theycorrespond to (a) checkerboard, (b) stripe, and (c)quarter-filled modulations. Sites with darker color arethe occupied (largely filled) sites and the others denoteempty (lightly filled) sites. (d)-(g) stand for possiblepatterns of arg[ (cid:104) b i (cid:105) ] in real space, and we call them ashomogeneous, checkerboard, collinear, and star signpatterns, respectively. Notice that each pattern isinvariant under a global Z transformation. II. RESULTSA. Hamiltonian and order parameters
We study the hard-core EBH model on a square lattice.The Hamiltonian is written by H = − t (cid:88) (cid:104) i,j (cid:105) ( b † i b j + H.C. ) − t (cid:48) (cid:88) (cid:104)(cid:104) i,j (cid:105)(cid:105) ( b † i b j + H.C. )+ V (cid:88) (cid:104) i,j (cid:105) n i n j + V (cid:88) (cid:104)(cid:104) i,j (cid:105)(cid:105) n i n j − µ (cid:88) i n i , (1)where b † i ( b i ) stands for the creation (annihilation) op-erator of the hard-core boson and n i = b † i b i . (cid:104) i, j (cid:105) and (cid:104)(cid:104) i, j (cid:105)(cid:105) denote the summation for the nearest-neighbor(nn) and next-nearest-neighbor (nnn) pairs, respectively. V and V represent the inter-site interactions and aregenerally taken to be positive (repulsive). Since we con-sider the hard-core limit, there can be only two possiblestates, the empty and occupied states, for each latticesite.In recent two papers [40, 41], the authors have studiedthe EBH model and found HSS with negative t (cid:48) , whichcauses the frustration of Hamiltonian. This motivatesus to wonder if there are other exotic SS or even SFphases once our model is frustrated in a different way,and if we are able to understand/predict the existence ofdistinct SF/SS phases. Therefore, in this work, we striveto provide a larger scope for underlying phases carried bythe EBH model (1) in the hard-core limit. To properlydistinguish phases, we define our order parameters in thissubsection.Since we focus on phases for the filling in 0 ≤ (cid:104) n i (cid:105) ≤ n ≤ .
5, we can onlyFIG. 2: (a) Phase diagram by iPEPS for ( t, V , V ) = (1 . , . , D = 8 results near the boundaries for better estimating the transition points.Rules above apply for every phase diagram in the following content. Upper and lower inset are diagrams from ref.41 and our mean-field theory under the same conditions. In the upper inset, solid (dotted) curves are boundaries forsecond-(first-)order transition. We demonstrate the schematic real-space structures of our states next to each phaseand their definitions for each symbol are the same as those in Fig. 1. Note that the sign in the circles corresponds toarg[ (cid:104) b i (cid:105) ] and a site with no sign indication means its (cid:104) b i (cid:105) =0. (b) Variations of order parameters for t (cid:48) = − . t = 1 is taken to be the energy unit.have three different crystal configurations: checkerboard(CB), stripe, and quarter-filled (QF) modulations. Here, (cid:104) n i (cid:105) means the particle density on site i . The configura-tion of the particles are drawn in Figs. 1(a)–1(c). For thesolid phases, CB and stripe contain two filled sites whileQF only has one. These solids emerge at commensuratefillings for the 2 × n ( k ) = 1 N C (cid:88) i ∈ C (cid:104) n i (cid:105) e i k · r i , (2)where C means unit cell and therefore N C is equal to 4for iPEPS. r i is the coordinate of location for each site.Depending on the patterns of these three solid phases,˜ n ( k ) can be non-zero at k = (0 , , π ), ( π, π, π ).The finite value of ˜ n ( π, π ) or ˜ n (0 , π ) (˜ n ( π, n (0 , n (0 ,
0) isalways non-zero in the following calculations. In order to distinguish those three configurations in Figs. 1(a)–1(c),we apply the same definition encoded in ref. 6. Thatis, for a state to be of CB pattern, we demand that its˜ n ( π, π ) (cid:54) = 0 and ˜ n ( π, n (0 , π ) = 0. For the rest twoorders, the stripe phase fulfills the condition | ˜ n ( π, − ˜ n (0 , π ) | (cid:54) = 0 and ˜ n ( π, π ) = 0, while the QF phase satisfies | ˜ n ( π, − ˜ n (0 , π ) | = 0 under that ˜ n ( π,
0) and ˜ n (0 , π ) (cid:54) = 0.Besides the structural order of particle density, to seeif a phase is SF/SS, one needs to examine its condensatedensity ρ , which is characterized by the non-zero valueof the off-diagonal components. In this paper, we simplydefine this order parameter as the following: ρ = 1 N C (cid:88) i ∈ C |(cid:104) b i (cid:105)| . (3)This reflects only the net SF density. In addition, weevaluate the sign of (cid:104) b i (cid:105) , namely arg[ (cid:104) b i (cid:105) ], for each site,because we expect the frustration between t and t (cid:48) to gen-erate various patterns of arg[ (cid:104) b i (cid:105) ] . Notice that the Fouriertransform for SF density is not performed because usuallythe modulation of |(cid:104) b i (cid:105)| is related to the order of particledensity. Moreover, because of the frustration, (cid:104) b i (cid:105) maysuffer from a sign difference in contrast to nearby sitesand therefore, the Fourier transform for a specific k inEq. (2) might fail to reflect the true situation when adetailed analysis is insufficient. As a result, our strategyis to see whether or not the condensate density (Eq. (3))can coexist with the solid order defined in Eq. (2) andcalculate (cid:104) b i (cid:105) for each site in the unit cell for further cat-egorizing each SS phase. Modulations of the SF densitywith different momenta are revealed by ED in Sec. II C.Configurations of possible patterns for arg[ b i ] are shownin Figs. 1(d)–1(g), named after homogeneous, diagonal,collinear, and star sign patterns, respectively. They willbe further discussed in Sec. II D. B. Phase diagram construction
In this section, we strive to construct a couple of phasediagrams in two extremely different component sets ofthe Hamiltonian (1), using iPEPS and mean-field theory.In the following phase diagrams obtained by iPEPS, weadopt a fixed bond dimension to be equal to 8. Thischoice already provides a good estimate for the phase di-agrams. We will discuss the influence of bond dimension D in the Appendix [49]. Also, we employ the mean-fieldtheory [50] for the Hamiltonian (1) and construct thephase diagrams for the 4 × T = 0for comparing with the iPEPS results. t = 1 , V = 8 , V = 0 First we consider the case for ( t, V , V ) = (1 . , . , t (cid:48) and µ . This is the same scenario dis-cussed in ref. 41, where the cluster mean-field approachis applied. We show our results in Fig. 2(a) togetherwith theirs.In Fig. 2(a), the phase boundaries denoted by filled(empty) blue circles indicate the second-(first-)orderphase transition. iPEPS is able to capture the first-order phase transitions quantitatively. When the com-putations are started near the first-order transition withdifferent initial states, each energy converges at a differ-ent value. By detecting the energy cross, we get to de-cide the transition points, as discussed in ref. 51. On theother hand, for the second-order phase transitions, or-ders can be formed continuously near the boundaries. Todetermine the first-order transition points and evaluaterelative errors, we perform the calculation along a verti-cal cut with fixed t (cid:48) and slice the value of µ in the unitof 0.1. If we find a energy level crossing between consec-utive µ = x and µ = x + 0 .
1, we then adopt µ = x + 0 . (cid:104) O (cid:105) (cid:54) = 0 if its valueis larger than 10 − . We define the second-order transi-tion point as the midpoint of two consecutive µ , where the order parameter becomes nonzero. Thus, the errorof the second-order transition point is also equal to 0.1.The upper inset in Fig. 2(a) is the result presentedin ref. 41, where the solid (dotted) curves represent thesecond-(first-)order transition. Our mean-field phase dia-gram is also shown in the lower inset. One can clearly seethat these phase diagrams coming from different meth-ods qualitatively agree with each other. We notice thatthe mean-field phase boundaries deviate from those ofthe iPEPS calculation. However, in Ref. 41, they haveshown that after scaling, phase boundaries move upwardand become closer to the iPEPS results. Thus, the mean-field calculation also captures the right phases. In fact,the mean-field analysis can be helpful in understandingthe background causes of each phase, which will be dis-cussed in Sec. II D.Two largest phases in Fig. 2(a) are the SF and CBSphases, separated by a first-order boundary and two SSphases. If a direct transition between SF and CBS takesplace, the transition is of first order due to the breaking ofdifferent symmetries. In contrast, if CSS stays in betweenthese two phases, the second-order transition is allowedat the phase boundary between SF and CSS, because itis characterized by the emergence of CBS that breaksthe Z symmetry. Thus, we expect that the criticality ofthis phase transition is explained by the Ising universal-ity. The CSS-CBS transition belongs to the SF-insulatortransition class [52]. As a result, it is also continuous.Among all the phases in Fig. 2(a), HSS is considered tobe one of the exotic phases induced by frustrated hoppingterm [40, 41]. This frustration is unavoidable because ofthe negative t (cid:48) . One of the characteristic feature of HSSis that the staggered pattern of arg[ (cid:104) b i (cid:105) ] appears alongthe nnn bonds. We demonstrate the phase transitionsfrom SF to HSS and from HSS to CBS for t (cid:48) = − . (cid:104) b i (cid:105) from ++ to + – in nearby sites.Therefore, the phase transition must be of first order.Since the superfluidity disappears at the HSS-CBS tran-sition, its universality belongs to the SF-insulator classand therefore the second-order transition occurs. t (cid:48) = 1 , V = 0 , V = 8 Now let us consider another completely different sce-nario. We set ( t (cid:48) , V , V ) = (1 . , , .
0) while varying t and µ . The phase diagram is shown in Fig. 3(a).As shown in the inset of Fig. 3(a), the mean-field the-ory qualitatively explains the phase diagram obtained byiPEPS. The solid state now has a stripe-like modulation,which is due to the nnn repulsive interaction ( V ). Forsuperfluid states, we now have two different states, thenormal superfluid (SF) and staggered superfluid (SSF)states. Although we can remove the frustration amongFIG. 3: (a) Phase diagram by iPEPS for ( t (cid:48) , V , V ) = (1 . , , . t = − . t = 0 . t (cid:48) = 1 is taken to be the energy unit.the hopping terms thanks to the positive t (cid:48) , the staggeredpattern of arg[ (cid:104) b i (cid:105) ] appears. In | t | (cid:38) .
4, they will evolveinto two separate SS states continuously by increasing µ ,and end up into the stripe solid phase.In Figs. 3(b) and 3(c), we demonstrate the variationsof order parameters for superfluid → supersolid → solidphases. The phase transition from SF(SSF) to SS1(SS2)is explained by the emergence of the stripe solid order.This criticality is expected to associate with the three-dimensional Ashkin-Teller (AT) model, where the weakfirst-order transition or second-order transition explainedby the Ising universality class takes place [53]. In thepresent iPEPS calculation, the transition seems to becontinuous. However, we need to investigate more com-putations near the phase boundaries for making a deci-sive claim, probably with full-update technique or iPEPSwith variational optimization. Further researches will beconsidered in future works. The phase boundary betweenSS1(SS2) and the stripe solid is of continuous transitionbecause the SF-insulator [52] type transition is expected.Combining with the observation upon Fig. 2(a), we con-clude the following two rules for phase transition: (I)Direct SF → solid transition is of first order, becauseof the breaking of different symmetries, according to the Ginzburg-Landau-Wilson paradigm. However, if we haveSS sandwiched in between SF and solid, then continuoussymmetry breaking can happen. And (II) if there is alocal phase change of (cid:104) b i (cid:105) , then the phase transition fallsinto the first order owing to the absence of a local Z gauge in our model. One last point to be noted is thatFigs. 3(b) and 3(c) look nearly identical. This reflectsthe symmetry of phase diagram under positive and neg-ative t . The reason will be discussed in the next section. C. Exact diagonalization
In this section, we present results obtained from thestudy of our complementary ED calculation. Due tothe fact that the finite-size effect becomes influentialnear the phase boundaries, instead of demonstrating thephase transitions, we focus on providing further evidencesof the existence for each of the nine phases in Figs.2(a) and 3(a). We pick one set of system parameters, p ≡ ( t, t (cid:48) , V , V ), and particle density n for each phase.To evaluate the values in the thermodynamic limit, weperform the finite-size analysis by using the scaling for-mula O ( N C ) = O ∞ + α/ √ N C . Cluster sizes used here (a) CSS, p = (1, 0.7, 8, 0), n = 0.375 (b) HSS, p = (1, 0.7, 8, 0), n = 0.25 N C (c) CBS, p = (1, 0, 8, 0), n = 0.5 n ( , )( n ( , 0) + n (0, ))/2 N C (d) SF, p = (1, 0, 8, 0), n = 0.125(0, 0)( , )( ( , 0) + (0, ))/2 FIG. 4: Finite-size scaling from ED calculation on fourunderlying states: (a)CSS, (b)HSS, (c)CBS, and (d)SF,within the phase diagram in Fig. 2(a). Parameterset-up p = ( t, t (cid:48) , V , V ) and averaged particle filling n are indicated as the figure titles. Interpolation isapplied between two particle numbers closest to ourdesired filling when such particle number isincommensurate to the lattice size. Dashed lines arefittings to the scaling formula O ( N C ) = O ∞ + α/ √ N C for the non-vanishing orders. (a) SS1, p = (0.7, 1, 0, 8), n = 0.375 (b) SS2, p = ( 0.7, 1, 0, 8), n = 0.375 N C (c) Stripe Solid, p = (0, 1, 0, 8), n = 0.5 n ( , )( n ( , 0) + n (0, ))/2 N C (d) SF, p = (0.5, 1, 0, 8), n = 0.125(0, 0)( , )( ( , 0) + (0, ))/2 FIG. 5: Same figures as those in Fig. 4 for (a)SS1,(b)SS2, (c)Stripe solid, and (d)SF within the phasediagram in Fig. 3(a). Note that in (c), yellow and greysquares overlap completely. This is due to the relationof orders at t = 0 discussed in the main text.are mainly N C = 16, 20, and 24. When no stripe-likeorder is expected, we include the results up to the 26-sitecluster in the finite-size analysis. Note that the 26-site cluster used here [54] is not compatible with the stripeorder under the periodic boundary condition.Since translational symmetry breaking can not be seendirectly in the ED ground states, we examine the pres-ence of the long-range order from the correlation func-tions. For the particle and SF density, we calculate theirstructural orders defined by S n ( k ) = 1( N C ) (cid:88) i,j ∈ C (cid:104) n i n j (cid:105) e i k · ( r j − r i ) (4)and S ρ ( k ) = 1( N C ) (cid:88) i,j ∈ C (cid:104) b † i b j (cid:105) e i k · ( r j − r i ) , (5)respectively. We then define the order parameters bytaking the square roots: ˜ n ( k ) = (cid:112) | S n ( k ) | and ˜ ρ ( k ) = (cid:112) | S ρ ( k ) | . Therefore, CB and stripe orders are repre-sented by ˜ n ( π, π ) and (˜ n ( π,
0) + ˜ n (0 , π )) /
2, respectively.The results for the finite-size scaling plot for( t, V , V ) = (1 , ,
0) and ( t (cid:48) , V , V ) = (1 , ,
8) are shownin Figs. 4 and 5, respectively. For ( t (cid:48) , V , V ) = (1 , , H ( ± t, t (cid:48) , V , V ) are related by the unitarytransformation of T b i T − = ( − r ix + r iy b i , which leavesevery parameter unchanged but inverts the sign of t inthe Hamiltonian. Notice that this transformation cannotremove the frustration induced by the kinetic terms andwe will discuss this issue in the next section. Using theabove unitary transformation, we can infer the results ofdesired orders as t → − t , by mapping ˜ ρ (0 , → ˜ ρ ( π, π )and ˜ ρ ( π, π ) → ˜ ρ (0 , ρ (0 ,
0) and˜ ρ ( π, π ) in Fig. 5(c), where t = 0 is chosen. Thus, we onlyshow the results for the four states except SSF in Fig. 5,because the outcomes for SSF are equivalent to those ofSF by the above-mentioned transformation.Figure 4 shows the cluster-size-dependence of thestructural order parameters for four sets of p and n . InFigs. 4(a) and 4(b), we demonstrate two SS states whoseorder parameters for the particle density show non-zeroCB order (˜ n ( π, π )). In Fig. 4(a), the SF density possessesthe same modulation as that of the particle density inCSS. In contrast, the SF density for (˜ ρ ( π,
0) + ˜ ρ (0 , π )) / (cid:104) b i (cid:105) ] in HSSis classified into the collinear type in Fig. 1(f). In theCB solid phase, all SF density orders are extrapolated tozero and only ˜ n ( π, π ) remains, as shown in Fig. 4(c). An-other state with only one existing order is the SF state inFig. 4(d), where all but ˜ ρ (0 ,
0) are extrapolated to valueseither close to or below zero.As mentioned above, Figs. 5(a) and 5(b) are symmet-rical after interchanging ˜ ρ (0 ,
0) and ˜ ρ ( π, π ) orders. Non-zero (˜ n ( π,
0) + ˜ n (0 , π )) / (cid:13) - 4 (cid:13) Hopping configurations within 2 by 2 unit cell, corresponding to homogeneous, diagonal, collinear,and star sign patterns, respectively. The sign in the red circle represents arg[ (cid:104) b i (cid:105) ]. Dashed squares enclose the areaof unit cells. The total nn (nnn) hopping energies for each configuration are shown in the table. For t >
0, the moststable configuration is 1 (cid:13) , while it becomes configuration 2 (cid:13) for t <
0. As for t (cid:48) , if it stays to be positive thenconfigurations 1 (cid:13) and 2 (cid:13) are equally favored but as t (cid:48) <
0, configuration 3 (cid:13) becomes the most stable.where the only remaining order is (˜ n ( π,
0) + ˜ n (0 , π )) / ρ (0 ,
0) survives in the thermodynamiclimit, which means that the SF phase appears. We con-clude that the results by ED are consistent with those byiPEPS and mean-field theory.
D. Discussion
1. Formation of SF/SS states
So far, we have demonstrated several SF/SS phases,besides the CSS and HSS which were found earlier[40, 41]. In this section, we try to clarify the mecha-nism of the emergence of CSS and HSS, and then makea simple categorization. To begin with, we have learnedthat introducing the diagonal nnn terms into the Hamil-tonian (1) is necessary for the formation of SS [26, 27].The CB (stripe) pattern is caused by the dominant V ( V ) interaction, because bosons tend to occupy the nnn(nn) sites in avoidance of repulsive interaction. When V competes with V , the third structural order, namely QForder, can be formed [6]. We will demonstrate the QForder in the latter section.What interests us the most is the phase modulationof SF density on each site, leading to the exotic SF/SSphases. In order to understand the underlying causes, wefirst decouple our hopping observable through a mean- field treatment: b i = (cid:104) b i (cid:105) + ( b i − (cid:104) b i (cid:105) ) = (cid:104) b i (cid:105) + δb i , (6)and then our hopping value can be approximated as: χ ij ≡ (cid:104) b † i b j (cid:105) = (cid:104) b † i (cid:105)(cid:104) b j (cid:105) + (cid:104) δb † i δb j (cid:105) ≈ (cid:104) b † i (cid:105)(cid:104) b j (cid:105) . (7)The approximation in Eq. (7) holds as away from thephase boundaries. By decoupling the hopping term inthe mean-field way, we are now more able to look insidewhat happens for different t and t (cid:48) .Because the SS state always starts from SF by enhanc-ing the chemical potential (filling), we shall first ana-lyze various SF states. In Figs. 1(d)–1(g), we presentfour possible SF states with different patterns of on-sitearg[ (cid:104) b i (cid:105) ]. Now we start our iPEPS calculation with an ini-tial state of certain filling and a hopping strength. For thesimplicity, we focus on the sign of |(cid:104) b † i b j + H.C. (cid:105)| ≡ χ i,j for ( i, j ) ∈ (cid:104) ij (cid:105) and (cid:104)(cid:104) ij (cid:105)(cid:105) and suppose that all χ i,j arehomogeneous. From Eq. (7), we can compose the pat-tern of χ i,j on each bond in the unit cell, as shown inFigs. 6 1 (cid:13) –6 4 (cid:13) , corresponding to homogeneous, diagonal,collinear, and star sign patterns, as indicated in Figs.1(d)–1(g). The overall nn (nnn) kinetic energies aftermultiplied by t ( t (cid:48) ) are shown in the table of Fig. 6.One can easily find out that in t > t (cid:48) >
0, theconfiguration with the lowest kinetic energy is homoge-neous pattern, meaning that no frustration happens. In t > t (cid:48) <
0, the total nn hopping energy E nnkin stillFIG. 7: (a) Phase diagram by iPEPS for ( t (cid:48) , V , V ) = ( − . , . , . d -wavesymmetry (dSF). When | t | is small, dSF enters the QF half supersolid (QFHSS) phase through a first-order phasetransition; while when | t | is large enough, QF solid (QFS) is sandwiched between QFHSS and dSF. All phases endup in the stripe solid phase when µ is large enough. We also have four narrow phases in between stripesolid/QFHSS and QFS/dSF when | t | is large. These four phases contain two stripe half supersolid (SHSS1/SHSS2)and two d -wave supersolid (dSS1/dSS2), indicated in the enlarged phase diagram sectors. (b) Phase diagram fromthe mean-field theory under the same conditions. (c) Variations of order parameters for t = 0 . | t (cid:48) | = 1 istaken to be the energy unit.prefers the homogeneous pattern; E nnnkin , however, is in fa-vor of the collinear pattern. In this way, the frustrationbetween E nnkin and E nnnkin realizes the HSS phase in Fig. 2.As shown in Fig. 3, when t (cid:48) >
0, homogeneous and di-agonal patterns are equally favored for t = 0, where thephase boundary between SF and SSF exists. In fact, theconfiguration with the lowest E nnkin is selected by the signof t . The condition of t > (cid:13) , while t < (cid:13) . This results in the reduction of frustrationbetween the kinetic terms. In the same manner, two SSstates stabilized by V can be formed with the differentpatterns of on-site arg[ (cid:104) b i (cid:105) ].The rules for constructing desired SF/SS phases aresummarized as follows: (I) Next-nearest-neighbor terms are needed for SS.(II) Dominant V ( V ) favors CB (stripe) pattern. QForder can be formed under competing V and V .(III) In t >
0, homogeneous pattern of arg[ (cid:104) b i (cid:105) ] is favoredas t (cid:48) >
0, while it competes with collinear pattern ofarg[ (cid:104) b i (cid:105) ] as t (cid:48) < t <
0, diagonal pattern of arg[ (cid:104) b i (cid:105) ] is favoredas t (cid:48) >
0, while it competes with collinear pattern ofarg[ (cid:104) b i (cid:105) ] as t (cid:48) < V / V for deciding theleading real-space order (Rule(II)). As for the patternof arg[ (cid:104) b i (cid:105) ], we can have the homogeneous/diagonal signpattern in t (cid:48) > t > t < t (cid:48) <
0, homo-geneous sign pattern competes with collinear sign patternin t > t < t . Here, according to Rule (III) and (IV),we can see that such transformation will leave our modelwith t (cid:48) < t (cid:48) >
0) to be still (non-)frustrated.So far, our results fit quite well with the above rulesfrom the mean-field theory. But more importantly, wewould like to see if these rules can be used to predict theexistence of certain SF or SS states. Therefore, in thenext section, we examine the scenario with the set-up of( t (cid:48) , V , V ) = ( − . , . , . t (cid:48) = − , V = 4 , V = 4 For t = −
1, the system includes the frustration be-tween the hopping terms, independent of the sign of t .The frustration cannot be removed by the unitary trans-formation discussed in Sec. II C. Before we look at the re-sults, we conjecture what kinds of states can appear fromour rules. Because V and V are competing, we expectthat QF order can be generated (Fig. 1(c)) [6, 37, 38].Also, we expect that the collinear pattern of arg[ (cid:104) b i (cid:105) ] isfavored in our SF state when | t (cid:48) | is dominant. This pro-vides the so-called d -wave superfluid (dSF) [55], whichshares the same d -wave pairing symmetry as the high- T c superconductivity [32]. By increasing µ , we expect theQF solid and QF SS phases before finally entering thestripe solid phase [6]. However, our SS is presumablydifferent from the conventional ones where no frustrationis presented.The resulting phase diagram is shown in Fig. 7(a).First, iPEPS phase diagram again qualitatively agreeswith the one by mean-field theory in Fig. 7(b). Here weisolate the mean-field phase diagram for better demon-stration. In our phase diagram, there appear one su-perfluid state (dSF), five supersolid states (dSS1, dSS2,QFHSS, SHSS1, SHSS2), and two solid states (QFS,stripe solid) in µ >
0. All our phases start with dSF,which was named after d -wave Bose liquid in the pre-vious works [55, 56]. But unlike their model, which in-cluded a four-site ring interaction, our dSF is realizedcompletely from the two-site terms. More importantly,the frustrated hopping terms play a crucial role in bothscenarios. When | t | (cid:46) .
3, dSF enters the QF super-solid phase, which is named after the QF half supersolid(QFHSS), through a first-order transition. For | t | (cid:38) . µ increases, QFHSS is replaced by the stripe solid after the first order transition. This isbecause the QF solid order is not simply composed oftwo perpendicular stripe orders. The first-order transi-tion from the QF phase to the stripe phase has also beenreported in Ref. 57.Our complex phase diagram also contains the stripe SSphases in several narrow regions, indicated by those fourpanels in Fig. 7(a). The difference between dSS1/dSS2and SHSS1/SHSS2 lies on the following point: the orien-tation of stripe solid order stays parallel/perpendicularto the orientation of sign modulation of arg[ (cid:104) b i (cid:105) ]. Thiscan be understood with again the mean-field analysis.Since we have two sublattices, we assign (cid:104) b sub1 (cid:105) = a ∈ C and (cid:104) b sub2 (cid:105) = b ∈ C . Then because of the inequality ofarithmetic and geometric means: a + b (cid:62) | a || b | , (8)in order to lower the total energy, when t >
0, theytend to possess the same sign among one sublattice sites,leading to dSS1 and SHSS1, and vice versa in t < ± t due to the same cause explained in Sec.II C. In Fig. 7(c), we demonstrate the variation of orderparameters for t = 0 .
9. The sudden changes of the orderparameter, indicating first-order transition, are very clearon the boundaries of dSS/QFS and QFHSS/SHSS. Dueto the emergence of the stripe order, the dSF-dSS1(dSS2)transition is again expected to be in association with thetype of 3D AT model [53], followed by the first-orderdSS1(dSS2)-QFS transition resulted from the breaking ofdifferent symmetries [6, 57]. The QFS-QFHSS transitionis continuous because it belongs to the SF-insulator class[52]. The QFHSS-SHSS1(SHSS2) transition is of first-order since it is again the QF-to-stripe transition. At last,SHSS1(SHSS2) transits to the stripe solid continuously,once again due to the SF-insulator class.
III. CONCLUSION
In this work, we have presented the exotic SF and SSstates out of the EBH model in the hard-core limit andhave discussed the background mechanism of generatingdesired SF/SS states, by manipulating the hopping andinteracting strengths of the EBH model. From the de-tailed mean-field interpretation, we have concluded withseveral rules of the SF/SS formation. Thanks to the ad-vance of cold atom experiment, which is now able to tunethe effective hopping parameters to be negative or evencomplex [58], our states and unveiled rules can be appliedexperimentally. By a simple mapping, the EBH modelcan be mapped into XXZ or Heisenberg model, wherethe frustration of interactive terms becomes importantfor deciding whether a system should be ferro- or an-tiferromagnetic. Therefore, from the side of hard-coreboson, there could be some interesting physics that mayhelp interpret the intriguing quantum magnetic states,0FIG. 8: Extrapolations of energies to infinite D for (a)SF to CBS and (b) SF to HSS phase transitions. Wechoose t (cid:48) = 0 . t (cid:48) = − . µ Extc . We show the transition points, µ D =8 c and µ Extc , in the lower panel of each figure.including the quantum spin liquid, more deeply. Thiswould be one of our future subjects.
IV. ACKNOWLEDGEMENT
W.-L.T. is supported by Postdoctoral ResearchAbroad Program, Project No. 108-2917-I-564-007, fromMinistry of Science and Technology (MOST) of Taiwan.Special thanks to Japan-Taiwan Exchange Associationfor its sponsorship of a short-term research activity in2018. H.-K.W. is supported by JQI-NSF-PFC (sup-ported by NSF grant PHY-1607611). T.S. is supportedby the Creation of new functional devices and high-performance materials to support next-generation indus-tries (CDMSI) and Challenge of Basic Science - Explor-ing Extremes through Multi-scale Simulations (CBSM2)from MEXT, Japan.
V. APPENDIX: IPEPS EXTRAPOLATION
Our phase diagrams shown in the main text come fromthe simple-update iPEPS calculation with fixed bond di-mension, D = 8. But since the precision of iPEPS ansatzis highly dependent on the bond dimension, it is requiredto carefully check its influence to our present model. Inthis Appendix, we extrapolate D to infinity and argue thephase boundaries obtained from the iPEPS calculations.Since the phase diagram in Fig. 2(a) contains both frus-trated and non-frustrated cases, we will focus our studieson those phase boundaries. FIG. 9: Extrapolation of order parameters to infinite D for (a) HSS to CBS, (b) CSS to CBS, and (c) SF toCSS phase transitions. We choose t (cid:48) = − . t (cid:48) = 0 . µ D =8 c and µ Extc are shown in the table forthese three phase transitions separately.
A. First-order boundaries
For the phase diagrams with fixed D , to determine thefirst-order phase boundaries, we searched along the verti-cal cut for consecutive µ where the energy level crossingtakes place for the adjoint phases. Therefore, we needto extrapolate the phase energies to infinite D and seewhether or not the choices of transition points for ourfixed- D phase boundaries are reasonable. We then plotthe energies for CBS and SF at t (cid:48) = 0 . t (cid:48) = − . D more largely. Therefore, we extrapolate the SF energieswith a polynomial fit up to the third order to infinite D and obtain E D →∞ SF . On the other hand, energies for HSSand CBS seem to be less affected by the choice of D ,so we simply choose the values of linear extrapolation as E D →∞ HSS and E D →∞ CBS . We compare the extrapolated val-ues and determine the phase transition point from thecondition, E D →∞ SF = E D →∞ HSS or E D →∞ SF = E D →∞ CBS .As shown in Figs. 8(a) and 8(b), the transition pointsare located at µ Extc ≈ .
79 and 0 .
3, respectively, whichare around 0.05 larger than the transition points ob-tained from the D = 8 calculation. This is not surpris-ing, since E D =8 SF is only a upper bound of energy and E D =8 SF > E D →∞ SF . Also, E CBS and E HSS are more stableupon varying D . Thus, the determined transition pointsafter the extrapolation becomes larger. However, since E D =8 SF is already a good estimate, the variation of transi-1FIG. 10: Phase diagram with transition pointsdetermined by fixed- D criteria (blue) and extrapolation(red). Solid (dotted) lines stand for second-(first-)ordertransition.tion points is not obvious. B. Second-order boundaries
Similar to the case for the first-order boundaries, weneed to extrapolate to infinite D for determining thesecond-order transition points. However, the quantitywe need here depends on the related order parameters.As shown in Fig. 9, we demonstrate the behavior of orderparameters as a function of 1 /D at the phase boundariesof (a) HSS to CBS, (b) CSS to CBS, and (c) SF to CSS. We choose t (cid:48) to be equal to -0.5 for Fig. 9(a) and 0.5for (b) and (c). The values of transition points are puttogether with previous ones for D = 8 in the table of Fig.9. One can see that the variation between these two tran-sition points becomes larger for the CSS-CBS transition.It is due to the fact that the related order parameter, ρ ,is less sensitive to the different bond dimension. C. Extrapolated phase diagram
Since we have already checked the extrapolated transi-tion points along certain cuts for both first- and second-order phase transitions, we would like to reconstructagain the phase diagram and compare it with that ofa fixed bond dimension. Based on Ref. [49], we assumethat the variation of the transition point is almost inde-pendent along the same boundary. Since we have alreadysampled one cut for each phase boundary, we shift theboundaries by using the same values obtained in Sec. IVA and Sec. IV B. The result is shown in Fig. 10. Onecan see that the variation is not very obvious except forthe second-order CSS-CBS transition. The reason is, asmentioned above, due to the fact that the related or-der parameter ( ρ ) is not very sensitive to the variationof D . Recall that for fixed- D phase diagram, we deter-mine the second-order transition points from the con-dition (cid:104) O (cid:105) > − for consecutive µ . If we apply thesame standard for the extrapolated phase diagram, then µ Extc ∼ .
9, which is much closer to µ D =8 c . Therefore, thevariation at the CSS-CBS transition point reflects the de-fect of our previous criteria for determining the transitionpoints. Nevertheless, it does not cause any of qualitativechange and thus is still acceptable. [1] M. Penrose and L. Onsager, “Bose-Einstein Condensa-tion and Liquid Helium,” Phys. Rev. , 576 (1956).[2] M. Boninsegni and N. V. Prokof’ev, “Colloquium: Su-persolids: What and where are they?” Rev. Mod. Phys. , 759 (2012).[3] A. Andreev and I. Lifshits, “Quantum Theory of Defectsin Crystals,” Sov. Phys. JETP , 1107 (1969).[4] G. V. Chester, “Speculations on Bose-Einstein Condensa-tion and Quantum Crystals,” Phys. Rev. A , 256 (1970).[5] A. J. Leggett, “Can a Solid be ”Superfluid”?” Phys. Rev.Lett. , 1543 (1970).[6] Y. C. Chen, R. G. Melko, S. Wessel, and Y. J. Kao,“Supersolidity from defect condensation in the extendedboson Hubbard model,” Phys. Rev. B , 014524 (2008).[7] E. Kim and M. H. W. Chan, “Probable observation of asupersolid helium phase,” Nature , 225 (2004).[8] E. Kim and M. H. W. Chan, “Supersolid Helium at HighPressure,” Phys. Rev. Lett. , 115302 (2006).[9] A. S. C. Rittner and J. D. Reppy, “Observation of Clas-sical Rotational Inertia and Nonclassical Supersolid Sig-nals in Solid He below 250 mk,” Phys. Rev.Lett. ,165301 (2006). [10] J. Day and J. Beamish, “Pressure-Driven Flow of SolidHelium,” Phys. Rev. Lett. , 105304 (2006).[11] J. Day and J. Beamish, “Low-temperature shear moduluschanges in solid He and connection to supersolidity,”Nature , 863 (2007).[12] D. Y. Kim and M. H. W. Chan, “Absence of Supersolidityin Solid Helium in Porous Vycor Glass,” Phys. Rev. Lett. , 155301 (2012).[13] I. Bloch, J. Dalibard, and W. Zwerger, “Many-bodyphysics with ultracold gases,” Rev. Mod. Phys. , 885(2008).[14] I. Bloch, J. Dalibard, and S. Nascimbne, “Quantum sim-ulations with ultracold quantum gases,” Nature Phys. ,267 (2012).[15] P. Windpassinger and K. Sengstock, “Engineering noveloptical lattices,” Rep. Prog. Phys. , 086401 (2013).[16] M. Tomza, K. Jachymski, R. Gerritsma, A. Negretti,T. Calarco, Z. Idziaszek, and P. S. Julienne, “Cold hy-brid ion-atom systems,” Rev. Mod. Phys. , 035001(2019).[17] T. Lahaye, C. Menotti, L. Santos, M. Lewenstein, andT. Pfau, “The physics of dipolar bosonic quantum gases,” Rep. Prog. Phys. , 126401 (2009).[18] O. Dutta, M. Gajda, P. Hauke, M. Lewenstein, D. S. Lh-mann, B. A. Malomed, T. Sowinski, and J. Zakrzewski,“Non-standard Hubbard models in optical lattices: a re-view,” Rep. Prog. Phys. , 066001 (2015).[19] L. Tanzi, E. Lucioni, F. Fam, J. Catani, A. Fioretti,C. Gabbanini, R. N. Bisset, L. Santos, and G. Modugno,“Observation of a Dipolar Quantum Gas with MetastableSupersolid Properties,” Phys. Rev. Lett. , 130405(2019).[20] F. Bttcher, J. N. Schmidt, M. Wenzel, J. Hertkorn,M. Guo, T. Langen, and T. Pfau, “Transient SupersolidProperties in an Array of Dipolar Quantum Droplets,”Phys. Rev. X , 011051 (2019).[21] L. Chomaz, D. Petter, P. Ilzhfer, G. Natale, A. Traut-mann, C. Politi, G. Durastante, R. M. W. van Bijnen,A. Patscheider, M. Sohmen, M. J. Mark, and F. Fer-laino, “Long-Lived and Transient Supersolid Behaviors inDipolar Quantum Gases,” Phys. Rev. X , 021012 (2019).[22] G. Natale, R. M. W. van Bijnen, A. Patscheider, D. Pet-ter, M. J. Mark, L. Chomaz, and F. Ferlaino, “Excita-tion Spectrum of a Trapped Dipolar Supersolid and ItsExperimental Evidence,” Phys. Rev. Lett. , 050402(2019).[23] L. Tanzi, S. M. Roccuzzo, E. Lucioni, F. Fam, A. Fioretti,C. Gabbanini, G. Modugno, A. Recati, and S. Stringari,“Supersolid symmetry breaking from compressional os-cillations in a dipolar quantum gas,” Nature , 382(2019).[24] M. Guo, F. Bttcher, J. Hertkorn, J. N. Schmidt, M. Wen-zel, H. P. Bchler, T. Langen, and T. Pfau, “The low-energy Goldstone mode in a trapped dipolar supersolid,”Nature , 386 (2019).[25] S. Baier, M. J. Mark, D. Petter, K. Aikawa, L. Chomaz,Z. Cai, M. Baranov, P. Zoller, and F. Ferlaino, “Ex-tended Bose-Hubbard models with ultracold magneticatoms,” Science , 201 (2016).[26] G. G. Batrouni and R. T. Scalettar, “Phase Separationin Supersolids,” Phys. Rev. Lett. , 1599 (2000).[27] F. Hbert, G. G. Batrouni, R. T. Scalettar, G. Schmid,M. Troyer, and A. Dorneich, “Quantum phase tran-sitions in the two-dimensional hardcore boson model,”Phys. Rev. B , 014513 (2001).[28] 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).[29] L. Savary and L. Balents, “Quantum spin liquids: a re-view,” Rep. Prog. Phys. , 016502 (2016).[30] Y. Zhou, K. Kanoda, and T. K. Ng, “Quantum spinliquid states,” Rev. Mod. Phys. , 025003 (2017).[31] Y. H. Chan and L. M. Duan, “Evidence of a spin liquidwith hard-core bosons in a square lattice,” New J. Phys. , 113039 (2012).[32] P. W. Anderson, “The Resonating Valence Bond Statein La CuO and Superconductivity,” Science , 1196–1198 (1987).[33] P. Corboz, T. M. Rice, and M. Troyer, “CompetingStates in the t − J Model: Uniform d -Wave State versusStripe State,” Phys. Rev. Lett. , 046402 (2014).[34] W. L. Tu and T. K. Lee, “Genesis of charge orders inhigh temperature superconductors,” Sci. Rep. , 18675(2016). [35] B. X. Zheng, C. M. Chung, P. Corboz, G. Ehlers, M. P.Qin, R. M. Noark, H. Shi, S. R. White, S. Zhang, andG. K. L. Chan, “Stripe order in the underdoped region ofthe two-dimensional Hubbard model,” Science , 1155(2017).[36] W. L. Tu and T. K. Lee, “Evolution of Pairing Or-ders between Pseudogap and Superconducting Phases ofCuprate Superconductors,” Sci. Rep. , 1719 (2019).[37] K. K. Ng and Y. C. Chen, “Supersolid phases in thebosonic extended Hubbard model,” Phys. Rev. B ,052506 (2008).[38] L. Dang, M. Boninsegni, and L. Pollet, “Vacancy super-solid of hard-core bosons on the square lattice,” Phys.Rev. B , 132512 (2008).[39] T. Suzuki and N. Kawashima, “Supersolid of hardcoreBosons on the face-centered cubic lattice,” Phys. Rev. B , 180502(R) (2007).[40] S. J. Dong, W. Liu, X. F. Zhou, G. C. Guo, Z. W.Zhou, Y. J. Han, and L. He, “Peculiar supersolid phasesinduced by frustrated tunneling in the extended Bose-Hubbard model,” Phys. Rev. B , 045119 (2017).[41] Y. C. Chen and M. F. Yang, “Two supersolid phasesin hard-core extended Bose-Hubbard model,” J. Phys.Commun. , 035009 (2017).[42] J. Jordan, R. Ors, G. Vidal, F. Verstraete, and J. I.Cirac, “Classical Simulation of Infinite-Size QuantumLattice Systems in Two Spatial Dimensions,” Phys. Rev.Lett. , 250602 (2008).[43] T. Nishino and K. Okunishi, “Corner Transfer MatrixRenormalization Group Method,” J. Phys. Soc. Jpn ,891 (1996).[44] R. Ors and G. Vidal, “Simulation of two-dimensionalquantum systems on an infinite lattice revisited: Cor-ner transfer matrix for tensor contraction,” Phys. Rev. B , 094403 (2009).[45] P. Corboz, “Variational optimization with infinite pro-jected entangled-pair states,” Phys. Rev. B , 035133(2016).[46] H. J. Liao, J. G. Liu, L. Wang, and T. Xiang, “Differen-tiable Programming Tensor Networks,” Phys. Rev. X ,031041 (2019).[47] H. C. Jiang, Z. Y. Weng, and T. Xiang, “Accurate De-termination of Tensor Network State of Quantum Lat-tice Models in Two Dimensions,” Phys. Rev. Lett. ,090603 (2008).[48] H. N. Phien, J. A. Bengua, H. D. Taun, P. Corboz,and R. Ors, “Infinite projected entangled pair states al-gorithm improved: Fast full update and gauge fixing,”Phys. Rev. B , 035142 (2015).[49] S. Wessel, I. Niesen, J. Stapmanns, B. Normand, F. Mila,P. Corboz, and A. Honecker, “Thermodynamic prop-erties of the Shastry-Sutherland model from quantumMonte Carlo simulations,” Phys. Rev. B , 174432(2018).[50] T. Matsubara and H. Matsuda, “A Lattice Model of Liq-uid Helium, I,” Prog. Theor. Phys. , 569 (1956).[51] P. Corboz, R. Ors, B. Bauer, and G. Vidal, “Simulationof strongly correlated fermions in two spatial dimensionswith fermionic projected entangled-pair states,” Phys.Rev. B , 165104 (2010).[52] M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S.Fisher, “Boson localization and the superfluid-insulatortransition,” Phys. Rev. B , 546 (1989). [53] G. Musial and J. Rogiers, “On the possibility of nonuni-versal behavior in the 3D Ashkin-Teller model,” Phys.Stat. Sol. (b) , 335 (2006).[54] Janez Jakliˇc and Peter Prelovˇsek, “Finite-temperatureproperties of doped antiferromagnets,” Advances inPhysics , 1–92 (2000).[55] O. I. Motrunich and M. P. A. Fisher, “ d -wave correlatedcritical Bose liquids in two dimensions,” Phys. Rev. B , 235116 (2007).[56] D. N. Sheng, O. I. Motrunich, S. Trebst, E Gull, andM. P. A. Fisher, “Strong-coupling phases of frustrated bosons on a two-leg ladder with ring exchange,” Phys.Rev. B , 054250 (2008).[57] K. K. Ng, “Thermal phase transitions of supersolids inthe extended Bose-Hubbard model,” Phys. Rev. B ,184505 (2010).[58] A. Eckardt, “Colloquium: Atomic quantum gases in pe-riodically driven optical lattices,” Rev. Mod. Phys.89