Table Of Content

Z spin liquid phase on the kagome lattice: a new saddle point 2 Tao Li Department of Physics, Renmin University of China, Beijing 100872, P.R.China (Dated: January 12, 2016) We haveperformed large scale variational search for the best RVBansatz for thespin-1 kagome 2 antiferromagnetwithbothnearest-neighboring(NN)andnext-nearest-neighboring(NNN)exchanges assuming only translational symmetry. We find the best RVB state is always fully symmetric and hasanmeanfieldansatzthatisgaugeequivalenttoapreviousproposedZ2 spinliquidansatz. The Z2 state is found to be slightly more stable than the extensively studied U(1) gapless Dirac spin 6 liquid state in both the J2 6=0 and the J2 =0 case and to possess a small spinon gap for J2 <0.2. 1 ThebreakingoftheU(1)gaugesymmetryintheZ2stateisfoundtoincreasewithJ2andtobequite 0 substantialaroundJ2 =0.15. However,wefindtheZ2 stateisalwaysveryclosetothegaplessU(1) 2 Dirac spin liquid state, although they have very different RVB parameters. We argue the kagome antiferromagnet should be better understood as a near critical system, rather than a system deep n inside a gapped spin liquid phase with well established Z2 topological order. a J PACSnumbers: 9 ] The search for spin liquid state in frustrated quan- Another scenario on the spin liquid ground state of l e tum antiferromagnets has attracted a lot of attention in the kagome antiferromagnets is provided by variational - the strongly correlated electron system community. The studies based on the RVB theory. It was found pre- r t spin-1 kagome antiferromagnet is a particularly promis- viously that the best RVB state for the ground state s 2 . ingsystemforthispurposebecauseofthestronggeomet- of the kagome antiferromagnet with NN exchange cou- t a rical frustration and the low coordinate number of the plingisaU(1)spinliquidstatewithaDirac-typespinon m lattice. Early studies on small clusters find the ground dispersion25,26. Extensive studies on this state find it - state of the kagome antiferromagnet with NN exchange is actually a very robust saddle point provided that the d has no signature of any symmetry breaking and is very symmetry of the system is unbroken27–29. Attempts to n likely a quantum spin liquid1–8. However, the exact na- break the U(1) gauge symmetry(and to introduce spin o ture of such a spin liquid phase is still under intense gap)alwaysresultinincreaseofenergyinthe thermody- c [ debate9,10. Especially,itis foundthatthe gapfortriplet namic limit. Morerecently,study onthe J1−J2 Heisen- excitationabovethegroundstateisverysmall,although bergmodelonthekagomelatticehasbeenperformedand 1 the spin correlation length is only of the order of the it is found that the U(1) Dirac spin liquid state remains v lattice constant. More curiously, it is find that the spin themoststablephaseevenwhenJ isquitelarge30. This 5 2 6 singlet channel of the system may host an exponentially is strange since strong evidence for Z2 topological order 1 large number of low energy excitations below the excep- has already been reported in DMRG simulation in this 2 tional small spin gap11–14. case. The U(1) Dirac spin liquid state is thus argued to 0 formastablephaseinthephasediagramaroundthe NN 1. The quest for the true nature of the ground state Heisenberg model point, rather than a single quantum 0 of the kagome antiferromagnet becomes even more ur- critical point. 6 gentafter its materialrealizationinthe ZnCu (OH) Cl 3 6 2 1 systems15–17. Recently,extensiveDMRGsimulationshas In this work, we reinvestigate the possible spin liq- : been performed on the kagome antiferromagnets both uidgroundstateofthe kagomeantiferromagnetwiththe v i withonlyNNexchangeandthatwithmoreextendedex- RVBtheory. Toreducebiasinourstudy,weonlyassume X change couplings18–24. The DMRG results seem to indi- translational symmetry for the RVB state at the begin- r catethatthekagomeantiferromagnetwithNNexchange ning. We have performed large scale variational search a has a small but finite gap in both spin triplet and spin for the best RVB ground state for the J −J model on 1 2 singlet channel, which is in support of a gapped Z spin the kagome lattice with up to the second neighbor RVB 2 liquid scenario. However, it is found that the ln2 topo- order parameters. We find the best RVB state is a fully logical entanglement entropy expected for a gapped Z symmetric spin liquid with a mean field ansatz that is 2 spinliquidcanonlybefaithfullydemonstratedwhenone gaugeequivalenttoaZ ansatzproposedpreviouslyfrom 2 introduce a finite NNN exchange coupling21. The abun- the analysis of projective symmetry group(PSG) on the dance of low energy singlet excitations in the system11 kagome lattice27. Different froma previous study on the also raises the suspicion that the DMRG results may same state, we find a finite spin gap and Z gauge struc- 2 not represent the intrinsic property of the system in the ture can be indeed be stabilized in both the J 6= 0 and 2 thermodynamic limit, since open boundary condition is the J =0 case, although the energy advantage over the 2 adopted in all DMRG simulations and symmetry break- U(1) gapless state is very small. We find the extent of ingperturbationsinthesingletchannelcanbesuperrel- U(1) gauge symmetry breaking increases with J . How- 2 evant. ever, the spin gap is found to follow the opposite trend 2 and vanishes around J = 0.2. The size of the spin gap In our study, we will adopt a L×L×3 cluster as il- 2 is found to be very small. lustrated in Fig.1. The doubled unit cell is illustrated The model we study in this paper is given by by the region within the pink parallelogram. According to the PSG of the state, the RVB order parameters on H =J1 Si·Sj +J2 Si·Sj (1) the dashed bonds should be multiplied by an additional <Xi,j> <<Xi,j>> minussignascomparedtotheRVBorderparameterson the solid bonds translated from them. The on-site RVB The motivation to introduce J is to perturb away the 2 order parameters should be manifestly translational in- systemfromthe U(1)Diracspinliquidphaseandtosta- variant. As a result, we have 3 independent on-site RVB bilizethepossibleZ spinliquidphase21–24. Here<i,j > 2 order parameters, 6 independent RVB order parameters denotes NN pair of sites, << i,j >> denotes NNN pair on the NN bonds and 6 independent RVB order param- of sites. To describe the possible Z spin liquid phase on 2 eters on the NNN bonds. Using the remaining gauge thekagomelattice,werewritethespinoperatorsinterms degree of freedom allowed by the translational symme- of the slave particle operators S = 1 f†σ f and 2 α,β α α,β β try, we can transform the RVB order parameters on the introduce the following mean field aPnsatz for the slave graysites andthe RVB orderparameterson the two NN Fermion f α bonds from the gray sites into the form of ατ . Finally, 3 we cantransformthe pairing orderparameteron a third H = ψ†U ψ (2) MF i i,j j NN bonds into real. Taking all these considerations into Xi,j account,we are left with 48 independent real variational parameters. inwhichψ =(f ,f† )istheNambuspinorconstructed i i,↑ i,↓ We have performed large scale variational optimiza- χ ∆∗ fromthe slaveFermion. Ui,j =(cid:18)∆ii,,jj −χi,∗ij,j (cid:19)is a2×2 otinonthoenktahgeomabeolvaettiRcVe.BWanesahtazvefoardtohpeteJd1t−heJS2tomcohdase-l matrixwithitsmatrixelementχ and∆ representing i,j i,j tic Reconfiguration method in our variational search32. the hopping and pairing type RVB order parameters for The optimized RVB ansatz is illustrated in Fig.1. More the ansatz. The RVB state is constructed by Gutzwiller specifically, the on-site RVB order parameter are given projectionofthe meanfield groundstate andisgivenby by |RVBi=P |{χ ,∆ }i, (3) G i,j i,j µτ gray sites 3 U = (4) here |{χ ,∆ }i denotes the meanfield groundstate of i,i i,j i,j (cid:26)µ~n ·~τ navy sites H and is a generalized BCS-type state. P denotes φ1 MF G the Gutzwiller projectionto removethe doublyoccupied in which ~τ = (τ ,τ ,τ ) are the Pauli matrixes, ~n = configuration in the mean field ground state. The re- 1 2 3 φ (sin(φ),0,cos(φ)) is a vectorof unit length in the τ −τ sultant RVB state can be simulated efficiently with the 1 3 plane. µ is a real number. The RVB order parameters variational Monte Carlo method. on the nearest-neighboring(NN) bonds are given by In our RVB state we will keep RVB order parameters U up to the second neighboring bonds. To represent i,j a fully symmetric spin liquid state, these RVB order pa- τ3 black bonds U =−s (5) rametersshouldbe invariantunderthe symmetryopera- i,j i,j(cid:26)~n ·~τ red bonds φ2 tions of the Hamiltonian up to some SU(2) gauge trans- formations of the form Ui,j −→ GiUi,jG†j. Here Gi is a Heresi,j =±1onthesolid(dashed)NNbonds. TheRVB SU(2) matrix acting on the Nambu spinor ψi. The sym- order parameter on the next-nearest-neighboring(NNN) metricRVBstatesoconstructedcanbeclassifiedintodis- bonds are given by tinct PSGs accordingto the gaugetransformations{G } i needed to restore the form of the RVB order parameters η~n ·~τ blue bonds afterthesymmetryoperations. Herewewillonlyassume Ui,j =−νi,j φ3 (6) (cid:26)η~n ·~τ green bonds the translational symmetry at the beginning. According φ4 to the scheme of the PSG classification, there are only two ways to realize the translational symmetry for a Z2 Here η is a real number, νi,j = ±1 on the solid(dashed) spin liquid31. In the first class, the mean field ansatz NNN bonds. should be manifestly translational invariant. In the sec- Theaboveansatzismanifestlytimereversalsymmetric ondclass,theunitcellofthemeanfieldansatzshouldbe but may in general break the point group symmetry of doubled and a Z gauge transformation is needed to re- the kagome lattice. However, we find at the optimum 2 storethetranslatedansatz. Inthisstudywewillconcen- of the variational energy, the following relations always trated on the second class since it is energetically much hold: φ =2φ ,φ =φ +φ . Insuchasituation,wecan 1 2 4 3 2 more favorable than the first class. It can also be shown actuallyprovethattheansatzdescribesafullysymmetric thatonlythesecondkindofRVBstatecanbeconnected spin liquid state. In fact, we can show that the above continuously to the U(1) Dirac spin liquid phase. RVBansatzisgaugeequivalenttotheZ ansatzproposed 2 3 -0.4 (a) 1.2 (b) 1.0 -0.6 0.8 2 0.6 -0.8 0.4 0.2 -1.0 0.0 0.00 0.04 0.08 0.12 0.16 0.20 0.00 0.04 0.08 0.12 0.16 0.20 J2 J2 3.0 0.30 2.5 (c) 0.25 (d) 2.0 0.20 3 1.5 0.15 1.0 0.10 0.5 0.05 0.0 0.00 0.00 0.04 0.08 0.12 0.16 0.20 0.00 0.04 0.08 0.12 0.16 0.20 FIG.1: ThemeanfieldansatzoftheZ2spinliquidstatestud- iedinthispaper. Theregionwithinthepinkparallelogram is FIG. 2: The optimized variational parameters as a function theunitcell ofthemeanfield Hamiltonian. Thedots(bonds) of J2. Thesystem size is L=12 in thecalculation in the same color have the same value for the RVB order . parameter Ui,j, except for an additional minus sign on the dashed bonds. parameters. However,our numerical optimization shows in Ref.[1] by the following gauge transformation thattheNNNRVBorderparameteriscrucialtostabilize a Z gauge structure. More specifically, if we set η = 0, 2 φ τ then in the optimized ansatz we always have φ =0. As exp(−i 2 2) i∈gray sites 2 G = 2 (7) we will see below, the optimized value of both φ2 and η i  exp(iφ2τ2) i∈navy sites increasewithJ2. ThisisconsistentwithpreviousDMRG studieswhichfindthattheconvergenceofthetopological 2  entanglement entropy becomes better for larger J2. The gauge transformed ansatz takes the form of Now we present our results of the variational calculation on the J −J model with the above RVB ansatz. 1 2 µ~n ·~τ on−site φ2 Our calculation is done on a L×L×3 cluster with L Ui,j = −si,jτ3 NN bonds (8) ranging from 12 to 24. We use periodic- anti-periodic −ν η~n ·~τ NNN bonds boundary condition in all calculations. Previous studies i,j φ3 on the same ansatz find that the U(1) Dirac spin liquid  When both φ2 and η are set to be zero, the above is the best variational state even for quite large J2 on ansatz is reduced to the U(1) Dirac spin liquid ansatz. sufficiently large system. What we find here is different. To break the U(1) gauge symmetry, we can introduce a We find the U(1) Dirac spin liquid is not the most sta- nonzeropairingtermeither onasite orona NNN bond. ble state for both zero and non-zero J2 on the largest Morespecifically,ifthereexiststwoclosedloopsstarting system we have attempted, which is much larger than from the same point, say, site i, on which the successive that is used in previous studies. In the fully symmetric productofthe RVBorderparametersP =U U ...U state, there are four variational parameters, namely µ, i i,j j,k l,i do not commute, the U(1) gauge symmetry is broken. φ2, φ3 and η to be determined by optimization of the To determine the gauge structure of the above ansatz, variationalenergy. In Fig.2 we plot the optimized values we can define the following two closed loops. The first of the variational parameters as functions of J2(J1 =1). oneisalength-zeroloopconsistingofthesiteionly. The Fromtheresults,weseetheextentofU(1)gaugesymme- productofRVB orderparametersonthisloopisjustthe trybreaking(whichisdictatedby the angleφ2)increases on-site RVB order parameter Pi1 = Ui,i = µ~nφ2 ·~τ. The monotonicallywithJ2. AtJ2 =0.15,thebreakingofthe second closed loop is the elementary triangle from site U(1) gauge symmetry is very significant. i. The product of RVB order parameters on this loop To show that the Z saddle point we find is indeed 2 is given by P2 = U U U = −τ . Here i,j,k are more stable than the U(1) saddle point found in previ- i i,j j,k k,i 3 the three vertices of the elementary triangle. Thus, the ous studies, we plot the interpolation of the variational U(1) gauge symmetry is broken when φ 6= 0 or π? . energybetweenthesetwosaddlepoints. TheU(1)saddle 2 The deviation of φ from 0 or π provides us a measure point is reached when we set φ = 0 and φ = 0. The 2 2 3 of the extent of the U(1) gauge symmetry breaking. For valueofη isingeneralnonzerofornonzeroJ andshould 2 the above ansatz, the U(1) gauge symmetry is broken beoptimized. Wedefineaninterpolationbetweentheop- when φ 6= 0,π even in the absence of NNN RVB order timized U(1) ansatz and the optimized Z saddle point 2 2 4 -0.43463 0.07 -0.43464 J2=0.15, L=12 0.06 -0.43465 0.05 p )-0.43466 ga0.04 E(-0.43467 on pin0.03 -0.43468 s 0.02 -0.43469 0.01 -0.43470 0.0 0.2 0.4 0.6 0.8 1.0 1.2 0.00 0.00 0.05 0.10 0.15 0.20 J2/J1 FIG. 3: The interpolation of the variational energy between the U(1) saddle point and the Z2 saddle point. The calcula- FIG. 5: The mean field spinon gap in the thermodynamic tion is done for J2 =0.15 on a L=12 system. limitcalculatedwiththeoptimizedRVBparametersobtained on the L=12 cluster. systems. Fig.4 shows the size dependence of φ up to 2 0.20 L = 24, which is the largest size that we can deal with. We find although there is a tendency in φ to decrease 2 0.16 with increasing L, the slope is very small and it is hard conclude if it vanish on any finite system that we can 0.12 treat. We note the largest L we have treated is already 2 significantly larger than the circumference of the cylin- 0.08 drical clusters used in DMRG calculations. Thus it is possiblethattheweaksignatureofZ spinliquidclaimed 2 byDMRGsimulationsatJ =0isstillafinitesizeeffect. 0.04 2 The true ground state at J = 0 in the thermodynamic 2 limitmaystillbeveryclosetoagaplessU(1)spinliquid, 0.00 10 12 14 16 18 20 22 24 or a quantum critical point33. L Nowlet’sinvestigatetheevolutionofthespingapwith J . Usually, the spin gap opens when the mean field 2 ansatzbreakstheU(1)gaugesymmetry. Wethusshould FIG. 4: The optimized value of φ2 at J2 = 0 for different expect the spin gap to increase with the extent of U(1) system sizes. gauge symmetry breaking. However, we find this is not true for our state. In principle, it is impossible to ex- tract the spin gap directly from the variational calcula- by mixing their variational parameters in the following tion on the ground state. The spin correlation length in way: x = αx(0) +(1−α)x(1). Here α is an interpola- i i i the groundstate is also nota goodindicatorfor the spin tionparameter,x(0) andx(1) denote thesetofoptimized excitation gap for the kagome antiferromagnet. Here we i i variational parameters(µ,φ2,φ3,η) for the U(1) and Z2 will use the mean field result for the spinon gap as an saddle point. The variationalenergy as a function of the indicatorofthespinexcitationgap. Theevolutionofthe interpolationparameterαatJ2 =0.15isshowninFig.3. mean field spinon gap with J2 is shown in Fig.5 (note The Z2 saddle point is obviously more robust than the themagnitudeoftheNNRVBorderparameterhasbeen U(1) saddle point. However, the energy difference be- set to be one). To reduce finite size effect, we have cal- tween the two saddle points is extremely small(or the culated the mean field dispersion in the thermodynamic order of 10−4J1 per site). This is very unusual since the limit with the optimized RVB parameters obtained on differenceinthegaugestructureofthetwosaddlepoints the L = 12 cluster. We find the spinon gap actually de- is so significant. creases with increasing J and approaches zero around 2 At the same time, we find the Z2 phase angle φ2 re- J2 = 0.2, a signature which may imply the proximity of mains nonzero even for J2 =0(note the optimized value the system to a magnetic ordered state at large J224. for η is significantly larger than previous report, which Finally,letusreturntothe problemoftheexceptional is 0.0186(2)28). To check if this is a finite size effect, insensitivityofthevariationalenergyintheRVBparam- wehaveperformedthevariationaloptimizationonlarger eters as we find in Fig.3. Similar behavior has also been 5 only when E > 1. This implies that the main difference between the two states is in their high energy behavior, which for some reason depend only very weakly on the structure of the ground state. This is possible U(1) for a multi-band system, in which totally different high nit) Z2 energy excitation spectrum can be realized on the same u ary ground state. It is interesting to see if similar behavior bitr alsooccurs for RVB states onother on-primitive Bravais ar lattices. E) ( D( In summary, we have performed large scale search for Z spin liquid state for the spin-1 kagome antiferromag- 2 2 netwithNNandNNN exchangesassumingonlytransla- tional invariance. We find the best RVB state is always 0 1 2 3 4 fully symmetric and is gauge equivalent to a Z ansatz 2 E proposed previously from symmetry analysis. The Z 2 state is found to be slightly more stable than the exten- FIG. 6: The density of state of the U(1) and Z2 ansatz at sivelystudiedU(1)gaplessDiracspinliquidstateinboth J2 =0.15. the J2 6= 0 and the J2 = 0 case and to possess a small spinongapforJ <0.2. WefindwhiletheextentofU(1) 2 gauge symmetry breaking in the Z state increases with 2 reported in previous variational studies of the kagome J , the spinon gap follows the opposite trend. These re- 2 antiferromagnet. A possible origin for such a strange sults are qualitatively consistent with the findings of the behavior is that the saddle points with distinct RVB DMRG simulations on the J −J model on the kagome 1 2 parameters may correspond to similar spin liquid state. lattice. However, we note the Z state we find is always 2 This implies that some RVB parameters may be quasi- very close to the gapless U(1) Dirac spin liquid state, al- redundant. To check if this is the case, we have calcu- though they have very different RVB parameters. Thus lated the overlapbetween the U(1) spin liquid state and we think the kagome antiferromagnet should be better the Z2 spin liquid state at J2 = 0.15. To our surprise, understood as a near critical system, rather than a sys- we find the overlap is larger than 0.99 on the L = 12 tem deep inside a gapped spin liquid phase with well cluster(the overlap at J2 is even larger). To find out the established Z2 topological order. origin for such a quasi-redundancy, we plot in Fig.6 the density of state of spinon excitation for the two ansatzs This work is supported by NSFC Grant No. 11034012 at J = 0.15. Except for the small gap feature around and the Research Funds of Renmin University of China. 2 E = 0.03 for the Z ansatz, the density of state of both The author acknowledge the discussion with Tomonori 2 ansatzs are almost identical below E = 1. The den- Shirakawa,SeijiYunoki,Yuan-MingLu,NormandBruce sity of state for both states become drastically different and Fa Wang in different stages of this work. 1 V. Elser, Phys.Rev.Lett. 62, 2405 (1989). 11 C. Waldtmann, H.-U. Everts, B. Bernu, C. Lhuillier, P. 2 J.T.ChalkerandJ.F.Eastmond,Phys.Rev.B46,14201 Sindzingre,P.Lecheminant,andL.Pierre,Eur.Phys.J.B (1992). 2, 501 (1998); P. Sindzingre, G. Misguich, C. Lhuillier, B. 3 P. W. Leungand V.Elser, Phys. Rev.B 47, 5459 (1993). Bernu,L.Pierre,Ch.Waldtmann,andH.-U.Everts,Phys. 4 N.ElstnerandA.P.Young,Phys.Rev.B50,6871(1994). Rev. Lett. 84, 2953 (2000); G. Misguich and B. Bernu, 5 P. Lecheminant, B. Bernu, C. Lhuillier, L. Pierre, and P. Phys. Rev. B, 71, 014417(2005); A. M. Läuchli and C. Sindzingre, Phys.Rev.B 56, 2521 (1997). Lhuillier, arXiv:0901.1065. 6 P. Sindzingre and C. Lhuillier, Europhys. Lett. 88, 27009 12 F.Mila,Phys.Rev.Lett.81,2356(1998); M.Mambriniand (2009). F. Mila, Eur. Phys. J. B 17, 651 (2000). 7 H. Nakano and T. Sakai, J. Phys. Soc. Jpn. 80, 053704 13 R. Budnik and A.Auerbach, Phys. Rev.Lett. 93, 187205 (2011). (2004). 8 A. M. Läuchli, J. Sudan, and E. S. Sørensen, Phys. Rev. 14 D. Poilblanc, M. Mambrini, and D. Schwandt,Phys. Rev. B 83, 212401 (2011). B 81, 180402 (2010). 9 R.R.P.Singh and D.A.Huse,Phys.Rev.Lett.68, 1766 15 J. S. Helton, K. Matan, M. P.Shores, E. A.Nytko,B. M. (1992); R. R. P. Singh and D. A. Huse, Phys. Rev. B 76, Bartlett, Y. Yoshida, Y. Takano, A. Suslov, Y. Qiu, J.-H. 180407(R) (2007); R. R. P. Singh and D. A. Huse, Phys. Chung,D.G.Nocera, and Y.S.Lee,Phys.Rev.Lett.98, Rev.B 77,144415 (2008). 107204 (2007). 10 G. Evenbly and G. Vidal, Phys. Rev. Lett. 104, 187203 16 P.Mendels,F.Bert,M.A.deVries,A.Olariu,A.Harrison, (2010). F.Duc,J.C.Trombe,J.S.Lord,A.Amato,andC.Baines, 6 Phys. Rev.Lett. 98, 077204 (2007). 25 M. B. Hastings, Phys. Rev.B 63, 014413 (2000). 17 T. H. Han, J. S. Helton, S. Chu, D. G. Nocera, J. A. 26 Y.Ran,M.Hermele,P.A.Lee,andX.G.Wen,Phys.Rev. Rodriguez-Rivera, C. Broholm, Y. S. Lee Nature 492, Lett. 98, 117205 (2007). 406410 (2012). 27 Y. M.Lu, Y. Ran, and P. A. Lee, Phys.Rev.B 83, 224413 18 H. C. Jiang, Z. Y. Weng, and D. N. Sheng, Phys. Rev. (2011). Lett. 101, 117203 (2008). 28 Y. Iqbal, F. Becca, and D. Poilblanc, Phys. Rev. B 83, 19 S. Yan, D. A. Huse, and S. R. White, Science 332, 1173 100404(R) (2011). (2011). 29 Y. Iqbal, F. Becca, and D. Poilblanc, Phys. Rev. B 84, 20 S.Depenbrock,I.P.McCulloch, andU.Schollwöck,Phys. 020407(R) (2011); New J. Phys. 14, 115031 (2012); Y. Rev.Lett. 109, 067201 (2012). Iqbal, F. Becca, S. Sorella, and D. Poilblanc, Phys. Rev. 21 H. C. Jiang, Z. Wang, and L. Balents, Nat. Phys. 8, 902 B 87, 060405(R) (2013); Y. Iqbal, D. Poilblanc, and F. (2012). Becca, Phys. Rev.B 89, 020407(R) (2014). 22 S. S. Gong, W. Zhu, and D. N. Sheng, Sci. Rep. 4, 6317 30 Y. Iqbal, D. Poilblanc, and F. Becca, Phys. Rev. B 91, (2014). 020402(R) (2015). 23 F.Kolley,S.Depenbrock,I.P.McCulloch, U.Schollwöck, 31 X. G. Wen,Phys. Rev.B 65 165113 (2002). V. Alba, Phys.Rev.B 91, 104418 (2015). 32 S.YunokiandS.Sorella, Phys.Rev.B74,014408 (2006). 24 S. S. Gong, W. Zhu, L. Balents, and D. N. Sheng, Phys. 33 Tao Li, arXiv:1106.6134. Rev.B 91, 075112 (2015).