Quantum criticality at finite temperature for two-dimensional JQ_3 models on the square and the honeycomb lattices
QQuantum criticality at finite temperature for two-dimensional
J Q models on thesquare and the honeycomb lattices J.-H. Peng, D.-R. Tan, L.-W. Huang, and F.-J. Jiang ∗ Department of Physics, National Taiwan Normal University, 88, Sec.4, Ting-Chou Rd., Taipei 116, Taiwan
We study the quantum criticality at finite temperature for three two-dimensional (2D) JQ modelsusing the first principle nonperturbative quantum Monte Carlo calculations (QMC). In particular,the associated universal quantities are obtained and their inverse temperature dependence are inves-tigated. The considered models are known to have quantum phase transitions from the N´eel orderto the valence bond solid. In addition, these transitions are shown to be of second order for twoof the studied models, with the remaining one being of first order. Interestingly, we find that theoutcomes obtained in our investigation are consistent with the mentioned scenarios regarding thenature of the phase transitions of the three investigated models. Moreover, when the temperaturedependence of the studied universal quantities is considered, a substantial difference between thetwo models possessing second order phase transitions and the remaining model is observed. Re-markably, by using the associated data from both the models that may have continuous transitions,good data collapses are obtained for a number of the considered universal quantities. The findingspresented here not only provide numerical evidence to support the results established in the litera-ture regarding the nature of the phase transitions of these JQ models, but also can be employed ascertain promising criterions to distinguish second order phase transitions from first order ones forthe exotic criticalities of the JQ -type models. Finally, based on a comparison between the resultscalculated here and the corresponding theoretical predictions, we conclude that a more detailedanalytic calculation is required in order to fully catch the numerical outcomes determined in ourinvestigation. PACS numbers:
I. INTRODUCTION
Studying the characteristics of transitions between var-ious phases of matters has long been an active researchtopic in the condensed matter physics [1–4]. In partic-ular, with the advances of experimental techniques andtheoretical approaches, some exotic phases and (or) ab-normal transitions between various states are observedand (or) proposed. One such a noticeable example is thetransition between the N´eel phase and the valence bondsolid [5].The so-called SU (2) JQ -type models of (quantum)spins in two-dimension (2D) [5–7], which will be intro-duced later, are systems hosting the exotic phase transi-tion between the N´eel and the valence bond solid states.The models have two types of couplings, namely the J -term and the Q -term. In addition, when the magnitudesof the couplings increase, the J -term and the Q -term fa-vor the formation of singlets and pairs of singlets (valencebond solid, VBS), respectively. On the one hand, whenthe magnitude of Q is much smaller than that of J , the SU (2) symmetry is broken spontaneously and the sys-tem is in the antiferromagnetic long-range order. On theother hand, the translation symmetry is broken in thecase of Q (cid:29) J . If one starts with a value of Q/J ( Q > J is fixed to be 1) so that the long-range order existsand the translation symmetry of the system is unbro- ∗ [email protected] ken, then by increasing the magnitude of Q one expectsthat the restoration of the broken SU (2) symmetry andthe breaking of the translation symmetry should occurat some couplings Q/J . In particular, these two distinctphysical incidents may take place at different values of
Q/J . If the two mentioned transitions occur at the same(
Q/J ) c , which is the most interesting scenario, then with-out fine-tuning this phase transition should be first orderaccording to the famous Landau-Ginsburg paradigm [8].Although the intriguing exotic phase transition intro-duced in the previous paragraph should be first order ingeneral, it is argue theoretically that this phase transi-tion, as well as those which belong to some other modelspossessing similar characteristics of the JQ -type models,can be second order [9, 10]. Such an exotic critical be-havior is referred to as the deconfined quantum criticality(DQC) in the literature.Many numerical investigations have been carried outin order to understand DQC, and various conclusions areobtained [11–17]. Although consistent critical exponentsare obtained by simulating square lattices as large as 448 [18, 19], one still cannot rule out the possibility that theconsidered phase transition of the SU (2) JQ model (onthe square lattice) is a weakly first order one. Indeed,it is shown in Ref. [20] that while the phase transitionof the ferromagnetic 5-state Potts model on the squarelattice is a weakly first order phase transition, yet a highquality data collapse with certain values of critical expo-nents, which should be a feature of a second order phasetransition, is reached. It is also unexpected that whenmore and more data of large lattices are included in the a r X i v : . [ c ond - m a t . s t r- e l ] F e b associated analyses, the determined values of ( Q/J ) c andthe correlation length exponent ν change accordingly ina sizeable manner [5, 12, 15, 18, 19]. For instance, whenthe largest linear system sizes are L = 32, 64 (128), 256,and 448, the calculated values of ν are given by 0.78,0.68, 0.446, and 0.448 (or 0.455), respectively. Hence, itwill extremely interesting to explore further the issue ofthe nature of these exotic phase transitions. In conclu-sion, at the moment it is a consensus that if these exoticphase transitions are second order, then they must re-ceive abnormally large corrections (to the scaling at thecorresponding criticalities).Studies of the SU (2) JQ -type models have been tar-geted at uncovering the nature of DQC. Moreover, due tothe famous Mermin-Wagner theorem [21–24], investigat-ing the finite-temperature properties of these models arefocusing on the VBS perspective. Yet by exploring thequantum critical regime (QCR) and the associated uni-versal quantities, which are features of the antiferromag-netic phase of the JQ -type models [25–36], one may gaincertain aspects of the nature of DQC [37–39]. Indeed,studying the QCR of various JQ -type models which haveeither a first order or a second order phase transitionsmay shed some light on the DQC.Due to the motivations described in the previous para-graphs, in this investigation we study three 2D SU (2) JQ models (Which will be introduced in detail later).In particular, two of these three systems are on the squarelattice and the remaining one is on the honeycomb lattice[6, 7, 17, 39]. Moreover, among these three considered JQ models, one is shown to have a first order phasetransition and the other two are argued to possess sec-ond order phase transitions [6, 7, 17, 39]. As a result,studying these three systems may reveal discrepancies ofthe features of QCR that result from the nature of thephase transitions.In our calculations, the associated universal quanti-ties of QCR, namely S ( π, π ) / ( χ s T ), χ u c /T and ρ s L/c are obtained. Here S ( π, π ), χ s , T , χ u , c , ρ s , and L are the staggered structure factor, the staggered sus-ceptibility, the temperature, the uniform susceptibility,the spinwave velocity, the spin stiffness, and the lin-ear system size, respectively. Remarkably, by simulat-ing large systems (448 spins for the square lattice andmore than 300 spins for the honeycomb lattice), we findthat there is a substantial difference in the temperaturedependence of these considered universal quantities be-tween models with first and second order phase transi-tions. Specifically, the features of QCR only show up forthe two JQ models which are shown to possess secondorder phase transitions in the literature. In addition, the S ( π, π ) / ( χ s T ) data of both the models having secondorder phase transitions reach the same plateau when re-garded as a function of the inverse temperature β . Asimilar scenario occurs for χ u c /T as well.Remarkably, by considering the S ( π, π ) / ( χ s T ) data ofboth the models having the characteristics of QCR as afunction of βc , a high quality data collapse shows up. J -terms Q -terms(a)(b)(c) FIG. 1: Three two-dimensional (2D) JQ models consideredin this study. For model (b), two more Q -terms, which canbe obtained by a 90 degree rotation of those depicted above,are not shown here. The quantity χ u c /T also demonstrates the same behav-ior when treated as a function of βc . The outcome associ-ated with ρ s L/c for the model on the honeycomb latticerequires a careful interpretation and can be satisfactoryaccounted for by taking into account a logarithmic cor-rection to the corresponding theoretical predictions.Although one cannot completely rule out the possibil-ity of weakly first order phase transitions, the obtainedresults shown here provide evidence to support the hy-pothesis that the targeted phase transitions, from theN´eel to the valence bond solid states, are likely continu-ous for two of the three studied JQ models.Finally, we also make a comparison between the nu-merical results obtained here and the related theoreticalpredictions in Refs. [27] and [37], and have concluded thatthe leading order large N -expansion calculations carriedout in Refs. [27] and [37] are not sufficient to match theoutcomes determined in our study.The findings presented in our investigation are not onlyinteresting in themselves, but also can be adopted as use-ful criterions to distinguish first order phase transitionsfrom the second order ones associated with DQC.The rest of this paper is organized as follows. In Sect.II the studied JQ models and the relevant observables ofQCR are introduced. We then demonstrate our numer-ical results, particularly the mentioned dramatical dif-ference of the temperature dependence of the calculateduniversal quantities are demonstrated in Sec. III. Finally,in Sec. IV we present discussions and conclusions. II. MICROSCOPIC MODELS ANDOBSERVABLES
The Hamiltonians of the three investigated 2D SU (2) JQ models on the square and the honeycomb latticeshave the same expression which is given by H = J (cid:88) (cid:104) ij (cid:105) (cid:126)S i · (cid:126)S j − Q (cid:88) (cid:104) ijklmn (cid:105) P ij P kl P mn ,P ij = 14 − (cid:126)S i · (cid:126)S j , (1)where in Eq. (1) J (which is set to be 1 here) and Q arethe couplings for the two-spin and the six-spin interac-tions, respectively, (cid:126)S i is the spin-1/2 operator at site i ,and P ij is the singlet pair projection operator betweennearest neighbor sites i and j . Figure 1 contains the car-toon representations of the studied models. In the follow-ing, the top, the middle, and the bottom models of fig. 1will be named the ladder, the staggered, and the hon-eycomb JQ models respectively, if no confusion arises.We would also like to emphasize the fact that althoughthe boundary condition (BC) for the honeycomb latticeconsidered here is different from that of Ref. [39], bulkproperties such as the critical point should be indepen-dent of the implemented BC. Hence, whenever ( Q/J ) c for the JQ model on the honeycomb lattice is needed,the one determined in Ref. [39] will be used.To study the quantum critical regime associated withthe investigated JQ models, particularly to calculate thecorresponding universal quantities, several relevant ob-servables, such as the staggered structure factor S ( π, π ),the staggered susceptibility χ s , the uniform susceptibility χ u , the spinwave velocity c , the spin stiffness ρ s , and thetemporal and spatial winding numbers squared ( (cid:104) W t (cid:105) and (cid:104) W i (cid:105) for i = 1 ,
2) are measured in our calculations.The formal expressions of these observables (on a L by L lattice) are as follows S ( π, π ) = 3 L L (cid:104) ( m zs ) (cid:105) , (2) m zs = (cid:88) i ( − i + i S zi , (3) χ s = 3 L L (cid:90) β (cid:104) m zs ( τ ) m zs (0) (cid:105) dτ, (4) χ u = βL L (cid:104) ( m z ) (cid:105) , (5) m z = (cid:88) i S zi , (6) ρ s = 34 β (cid:0) (cid:104) W (cid:105) + (cid:104) W (cid:105) (cid:1) , (7)where β and L i for i = 1 , i -spatial direction)used in the simulations, respectively, and S zi is the third-component of the spin-half operator (cid:126)S i at site i . III. THE NUMERICAL RESULTS
To calculate the desired physical observables, we havecarried out large-scale quantum Monte Carlo calculations(QMC) using the stochastic series expansion (SSE) algo-rithm with very efficient operator-loop update [40, 41].In addition, the simulations are performed at (close to)the associated critical points established in the literature.Specifically, we conduct the simulations at (
Q/J ) c =1 . JQ models, respectively. Finally,outcomes of up to 448 (over 300 ) spins on the square(honeycomb) lattice are obtained.The calculated numerical results regarding the univer-sal quantities introduced previously will be described indetail in the following. Particularly, in the analyses weassume that the associated dynamic critical exponents z are all given by z = 1 [9, 10]. A. S ( π, π ) / ( χ s T ) The observable S ( π, π ) / ( χ s T ) as functions of β for L = 24, 32, 40,..., 192, 256, and 448 (from top to bottom)for the ladder and the staggered (only up to L =256) JQ models are demonstrated in fig. 2. Similarly, fig. 3 con-tains the results of S ( π, π ) / ( χ s T ) versus β for the hon-eycomb JQ model.As can been seen from the top panel of fig. 2 whichis associated with the ladder JQ model, although mild(finite) temperature and system size dependence maystill be there for the largest lattice, it is clear that the S ( π, π ) / ( χ s T ) of L = 448 reaches a (more or less)plateau, i.e., a constant value close to 1.25, for β > β dependence of S ( π, π ) / ( χ s T ) for the staggered JQ model is demonstrated in the bottom panel of fig. 2.The results in the figure indicate S ( π, π ) / ( χ s T ) reachesthe value of 1 rapidly, both at high and moderate lowtemperatures. Moreover, no sign of the appearance of aplateau implies the absence of QCR. While further ver-ification is on request, most likely this is a feature of afirst order phase transition.Similar to the outcomes of the ladder JQ model (onthe square lattice), the quantity S ( π, π ) / ( χ s T ) of thehoneycomb JQ model reaches a plateau for β >
15, seefig. 3.It is interesting to notice that the results shown infigs. 2 and 3 imply that when compared with the laddermodel, larger values of β are required for the quantity S ( π, π ) / ( χ s T ) of the honeycomb JQ model to saturateto its plateau.Remarkably, the plateau for both the ladder and thehoneycomb JQ models take the same value. In partic- S ( , ) / ( s T ) L = 448 S ( , ) / ( s T ) L = 256 FIG. 2: S ( π, π ) / ( χ s T ) as functions of β for various L for theladder (top) and the staggered (bottom) JQ models. S ( , ) / ( s T ) L = 116, L = 68 L = 136, L = 80 L = 272, L = 160 L = 408, L = 240 FIG. 3: S ( π, π ) / ( χ s T ) as functions of β for various L for thehoneycomb JQ models. ular, if the S ( π, π ) / ( χ s T ) data on large lattices of bothmodels are plotted as functions of βc (The calculationsof c will be detailed shortly), a good data collapse out-come is obtained, see fig. 4. This indicates it is highlyprobable that QCR is indeed found in these two modelsand such an observation in turn provides an evidence forthe associated phase transitions to be continuous.
20 40 60 80 100 120 c S ( , ) / ( s T ) L = 256, square L = 448, square L = 272, L = 160, honeycomb L = 408, L = 240, honeycomb FIG. 4: S ( π, π ) / ( χ s T ) as functions of βc for the ladder andthe honeycomb JQ models. The scenarios reached here are consistent with theknown outcomes in the literature regarding the exoticphase transitions of these three investigated JQ mod-els. B. χ u c /T
1. Determination of the spinwave velocity c Since the low-energy constant spinwave velocity c isrequired for calculating the numerical value of χ u c /T ,we have determined c firstly.The method used here for the calculations of c is thewinding numbers squared [35, 37, 42, 43]. Specifically,for each considered box size L , the β is adjusted in thesimulations so that three winding numbers squared (cid:104) W i (cid:105) for i = 1 , (cid:104) W t (cid:105) take the same value. When thiscondition is met, the c corresponding to the simulated L is L/β .It should be pointed out that on the honeycomb lat-tice, square-shape area can only be reached with certainaspect ratios L /L . As a result, the simulations associ-ated with the honeycomb lattice are performed on theseaspect ratios (of L /L ) so the (almost) square-shapearea is obtained.Finally, the calculations are done at Q/J = 1.4925,1.1923, 1.19 for the ladder, the staggered, and the hon-eycomb JQ models, respectively.The spatial and the temporal winding numberssquared as functions of β for the ladder and the stag-gered JQ models are shown in the top ( L = 12) and thebottom ( L = 20) panels of fig. 5. In addition, the esti-mated c related to various L for the ladder JQ modelis demonstrated in fig. 6. A fit of applying the ansatz c + a/L + b/L ( c + a/L + b/L + d/L ) to the data offig. 6 leads to c = 5 . c = 5 . JQ model. This is because as we will show shortly,with the fact that c is a constant, the related χ u c /T does W x W t W x W t FIG. 5: (cid:104) W x (cid:105) and (cid:104) W t (cid:105) as functions of β for the ladder (top, L =12) and the staggered (bottom, L =20) JQ models. Here (cid:104) W x (cid:105) = 0 . (cid:0) (cid:104) W (cid:105) + (cid:104) W (cid:105) (cid:1) . L c FIG. 6: c as a function of L for the ladder JQ models. not possess the expected feature of QCR. As a result, cal-culating the bulk c for this model is no long needed.Regarding the c of the honeycomb JQ model, the out-come on a lattice of L = 136 and L = 80 is quoted hereand is given by c ∼ .
92. We have also carried out simu-lations on a 34 ×
20 honeycomb lattice and the resulting c is estimated to be c = 2 . . / .
92 is over 0.97 and the ratio of the linear box sizes of these twohoneycomb lattices is √ × / √ ×
20 = 4, one ex-pects negligible finite-lattice effects for the quoted valueof c ∼ .
2. Determination of χ u c /T After determining the values of the spinwave velocity c for all the three investigated JQ model, we turn tostudying the universal quantity χ u c /T . χ u c /T as a function of β for the ladder JQ modelis shown in fig. 7. The figure suggests convincingly thatthe observable χ u c /T of this model reaches a plateaufor β ≥
2, which is a feature of QCR. Similar to that ofthe quantity S ( π, π ) / ( χ s T ), this characteristic providesyet another evidence that the corresponding phase tran-sition is likely a second order one. The horizontal solid(2.7185) and dashed (1.7125) lines in fig. 7 represent thetheoretical calculations obtained in Refs. [27] and [37],respectively. It is clear that the next to leading orderlarge N -expansion computations carried out in these ref-erences are not sufficient to match the numerical value of χ u c /T determined in this study.In fig. 8 χ u /T is considered as a function of β for thestaggered JQ model. For a first order phase transi-tion, one expects the magnitude of the winding numbersquared grows with β . Hence the rising of χ u /T in thelow temperature region is a signal of first order phasetransition for the staggered JQ model. This is consis-tent with the conclusion established in the literature.Similar to the outcomes associated with S ( π, π )/( χ s T ),the plateaus of χ u c /T for both the ladder and the hon-eycomb JQ models have the values of about the samemagnitude, and a good data collapse result is obtained ifthe data on large lattices of both models are consideredas functions of βc , again see fig. 9. This observation canalso be regarded as an evidence that both the transitionsof the ladder and the honeycomb JQ model are likelycontinuous.Interestingly, if the data of χ u c /T of both the ladderand the honeycomb JQ models are plotted as functionsof β , a slightly downgrade quality scaling than that offig. 9 appears as well.Finally, it should be pointed out that the trend of mov-ing downward for the quantity χ u c /T at large valuesof β shown in figs. 7 and 9 is due the observable χ u and is an artifact. In particular, to obtain the correctbulk χ u , taking the appropriate order of limits, namely,lim β →∞ lim L →∞ is essential. C. The spin stiffness ρ s For the ladder JQ model, the quantity ρ s L/c as func-tions of β for several L are demonstrated in fig. 10. Theresults in fig. 10 indicate that in the limit LT → ∞ , ρ s L/c approaches a value between 0.32 and 0.33. This u c / T L = 448 FIG. 7: χ u c /T as functions of β for various L for the ladder JQ models. The horizontal solid and dashed lines representthe outcomes of 2.7185 and 1.7125 determined in Refs. [27]and [37], respectively. u / T L = 256 FIG. 8: χ u /T as a function of β for various L for the staggered JQ models. c u c / T L = 256, square L = 448, square L = 272, L = 160, honeycomb L = 408, L = 240, honeycomb FIG. 9: χ u c /T as functions of βc for various L for the ladderand the honeycomb JQ models. s L / c L = 16 L = 32 L = 48 L = 64 L = 96 FIG. 10: ρ s L/c as functions of β for various L for the ladder JQ models. s L L = 16 L = 48 L = 96 L = 192 FIG. 11: ρ s L/c as functions of β for various L for the stag-gered JQ models. number differ significantly from 0 . ρ s = β (cid:0) (cid:104) W (cid:105) + (cid:104) W (cid:105) (cid:1) is used, then our calculationslead to ρ s L/c ∼ .
215 in the LT → ∞ limit, which isclose to the predicted result of 0 . JQ model, it is anticipated that ρ s L/c will grow with L (in the LT → ∞ limit), and thisis indeed what’s observed in our data, see fig. 11.In fig. 12, ρ s L/c as functions of β is shown for severalhoneycomb lattices. Although no convergence (as LT →∞ ) is observed, it is interesting to notice that a good datacollapse emerges when a logarithmic correction is takeninto account. Specifically, if ( ρ s L/ log ( L/L )) /c, L =0 .
37 is considered as functions of βc , then a high qualityscaling appears, see fig. 13. This scenario is consistentwith that obtained in Ref. [39]. Finally, since a free pa-rameter is allowed in the logarithmic correction, namelyone can use a log ( L/L ) with a being some constant in-stead of log ( L/L ), we do not contrast the results relatedto the spin stiffness ρ s on the honeycomb lattice with thaton the square lattice. s L / c L = 34, L = 20 L = 68, L = 40 L = 136, L = 80 FIG. 12: ρ s L/c as functions of β for various L for the honey-comb JQ models. s L / c / ( l o g ( L / L )) L = 34, L = 20 L = 68, L = 40 L = 136, L = 80 FIG. 13: ρ s L/c/ (log (
L/L )) as functions of β for various L for the honeycomb JQ models. L is 0.37. IV. DISCUSSIONS AND CONCLUSIONS
Using the first principle quantum Monte Carlo simu-lations, we investigate the QCR of three 2D SU (2) JQ models on both the square and the honeycomb lattices.Particularly, three universal quantities associated withQCR, namely S ( π, π ) / ( χ s T ), χ u c /T and ρ s L/c are ob-tained and their β ( βc ) dependence is studied.Among the three studied models, the ladder and thehoneycomb JQ models are shown to possess second or-der phase transitions from the N´eel phase to the VBSphase [6, 17, 39], and for the staggered JQ model theassociated transition is first order [7]. Our investigationrelated to the QCR leads to results consistent with these scenarios.Specifically, while for the ladder and the honeycomb JQ models, the quantity S ( π, π ) / ( χ s T ) ( χ u c /T ) ofboth models reach the same plateau value and a goodscaling emerges when considered as a function of βc ,no such behavior is found for the staggered JQ mod-els. Moreover, ρ s L/c obtains a constant outcome in the LT → ∞ limit for the ladder model. Similar situationapplies to the honeycomb model when a logarithmic cor-rection to the observable ρ s L/c is taken into account.For the staggered model both the quantities χ u c /T and ρ s L/c grow with L or β . This is in agreement with thefact that the targeted phase transition of the staggeredmodel is first order.It is surprising to observe that to obtain a good scaling,the universal quantity associated with ρ s on the honey-comb lattice requires a logarithmic correction. This is notthe case for the same observable on the square lattice. Itwill be interesting to understand this from a theoreticalperspective.Apart from the three universal quantities consideredhere, one frequently studied observable associated withQCR is the Wilson ratio W = χ u T /C V . Here C V is thespecific heat and can be calculated directly or through theinternal energy density E . However, due to the subtletyof both approaches for the determination of C V shown inRefs. [35, 44], here we do not attempt to carry out a highprecision calculation of the quantity W . Nevertheless,the investigation conducted here provide evidence to sup-port the scenario that for the ladder and the honeycomb JQ models, the associated quantum phase transitionsfrom the N´eel phase to the VBS phase are likely to becontinuous.The results presented here are not only interesting inthemselves, but also can be used as criterions to dis-tinguish second order phase transitions from first orderones for exotic phase transitions such as the DQC studiedhere.This study is partially supported by the Ministry ofScience and Technology of Taiwan. [1] Nigel Goldenfeld, Lectures On Phase Transitions AndThe Renormalization Group (Frontiers in Physics) (Addison-Wesley, 1992).[2] J. Cardy,
Scaling and Renormalization in StatisticalPhysics , Cambridge University Press, Cambridge, UK, 1996.[3] Lincoln D. Carr,
Understanding Quantum Phase Transi-tions (Condensed Matter Physics) (CRC Press, 2010).[4] S. Sachdev,
Quantum Phase Transitions (CambridgeUniversity Press, Cambridge, 2nd edition, 2011). [5] A. W. Sandvik, Phys. Rev. Lett. , 227202 (2007).[6] J. Lou, A. W. Sandvik, and N. Kawashima, Phys. Rev.B. 80, 180414(R) (2009).[7] Arnab Sen and Anders W. Sandvik Phys. Rev. B ,174428 (2010).[8] V.L. Ginzburg and L.D. Landau, Zh. Eksp. Teor. Fiz. , 1064 (1950). English translation in: L. D. Landau,Collected papers (Oxford: Pergamon Press, 1965) p. 546.[9] T. Senthil, L. Balents, S. Sachdev, A. Vishmanath andM. P. A. Fisher, Science , 144407 (2004).[11] Anatoly Kuklov, Nikolay Prokof’ev, Boris Svistunov,Matthias Troyer, Annals of Physics Volume 321, Issue7, July 2006, Pages 1602-1621.[12] R. G. Melko and R. K. Kaul, Phys. Rev. Lett. ,017203 (2008).[13] F.-J. Jiang, M. Nyfeler, S. Chandrasekharan, and U.-J.Wiese, J. Stat. Mech. (2008) P02009.[14] A. B. Kuklov, M. Matsumoto, N. V. Prokof’ev, B. V.Svistunov, M. Troyer, Phys. Rev. Lett. , 050405(2008).[15] A. W. Sandvik, Phys. Rev. Lett. , 177201 (2010).[16] Kun Chen, Yuan Huang, Youjin Deng, A. B. Kuklov, N.V. Prokof’ev, and B. V. Svistunov Phys. Rev. Lett. ,185701 (2013).[17] Kenji Harada, Takafumi Suzuki, Tsuyoshi Okubo,Haruhiko Matsuo, Jie Lou, Hiroshi Watanabe, SyngeTodo, and Naoki Kawashima Phys. Rev. B 88, 220408(R)(2013).[18] Hui Shao, Wenan Guo, Anders W. Sandvik, Science 352,213 (2016).[19] Anders W. Sandvik, Bowen Zhao, Chin. Phys. Lett. ,057502 (2020)[20] Shumpei Iino, Satoshi Morita, Anders W. Sandvik, NaokiKawashima, J. Phys. Soc. Jpn. 88, 034006 (2019).[21] N. D. Mermin and H. Wagner, Phys. Rev. Lett. , 1133(1966).[22] P. C. Hohenberg, Phys. Rev. , 383 (1967).[23] Sidney Coleman, Communications in MathematicalPhysics volume 31, pages 259–264 (1973).[24] Axel Gelfert and Wolfgang Nolting, J. Phys.: Condens. Matter 13 R505 (2001).[25] A. V. Chubukov and S. Sachdev, Phys. Rev. Lett. ,169 (1993).[26] A. V. Chubukov and S. Sachdev, Phys. Rev. Lett. ,2680 (1993).[27] A. V. Chubukov, S. Sachdev, and J. Ye, Phys. Rev. B , 11919 (1994).[28] A. W. Sandvik, A. V. Chubukov, and S. Sachdev, Phys.Rev. B , 16483 (1995)[29] M. Troyer, H. Kantani, and K. Ueda, Phys. Rev. Lett. , 3822 (1996).[30] Matthias Troyer, Masatoshi Imada, and Kazuo Ueda, J.Phys. Soc. Jpn. 66, 2957 (1997).[31] Jae-Kwon Kim and Matthias Troyer, Phys. Rev. Lett. , 2705 (1998).[32] Y. J. Kim, R. J. Birgeneau, M. A. Kastner, Y. S. Lee,Y. Endoh, G. Shirane, and K. Yamada, Phys. Rev. B 60,3294 (1999).[33] Y. J. Kim and R. J. Birgeneau, Phys. Rev. B , 6378(2000).[34] A. W. Sandvik, V. N. Kotov, and O. P. Sushkov, Phys.Rev. Lett. , 207203 (2011).[35] A. Sen, H. Suwa, and A. W. Sandvik, Phys. Rev. B ,195145 (2015).[36] D.-R. Tan and F.-J. Jiang, Phys. Rev. B , 245111(2018).[37] Ribhu K. Kaul, Roger G. Melko, Phys. Rev. B , 014417(2008).[38] Ribhu K. Kaul and Subir Sachdev, Phys. Rev. B ,155105 (2008).[39] Sumiran Pujari, Kedar Damle, and Fabien Alet, Phys.Rev. Lett. , 087203 (2013).[40] A. W. Sandvik, Phys. Rev. B , R14157 (1999).[41] A. W. Sandvik, AIP Conf. Proc. 1297, 135 (AIP, NewYork, 2010).[42] F.-J. Jiang, Phys. Rev. B , 024419 (2011).[43] F.-J. Jiang and U.-J. Wiese, Phys. Rev. B , 155120(2011).[44] Jhao-Hong Peng, L.-W. Huang, D.-R. Tan, and F.-J.Jiang, Phys. Rev. B101