Plaquette Ordered Phase and Quantum Phase Diagram in the Spin-1/2 J1-J2 Square Heisenberg Model
Shou-Shu Gong, Wei Zhu, D. N. Sheng, Olexei I. Motrunich, Matthew P. A. Fisher
PPlaquette Ordered Phase and Quantum Phase Diagram in the Spin- J - J Square HeisenbergModel
Shou-Shu Gong , Wei Zhu , D. N. Sheng , Olexei I. Motrunich , Matthew P. A. Fisher Department of Physics and Astronomy, California State University, Northridge, California 91330, USA Department of Physics, California Institute of Technology, Pasadena, California 91125, USA Department of Physics, University of California, Santa Barbara, California 93106-9530, USA
We study the spin- Heisenberg model on the square lattice with first- and second-neighbor antiferromag-netic interactions J and J , which possesses a nonmagnetic region that has been debated for many years andmight realize the interesting Z spin liquid. We use the density matrix renormalization group approach with ex-plicit implementation of SU (2) spin rotation symmetry and study the model accurately on open cylinders withdifferent boundary conditions. With increasing J , we find a N´eel phase and a plaquette valence-bond (PVB)phase with a finite spin gap. From the finite-size scaling of the magnetic order parameter, we estimate that theN´eel order vanishes at J /J (cid:39) . . For . < J /J < . , we find dimer correlations and PVB textureswhose decay lengths grow strongly with increasing system width, consistent with a long-range PVB order inthe two-dimensional limit. The dimer-dimer correlations reveal the s -wave character of the PVB order. For . < J /J < . , spin order, dimer order, and spin gap are small on finite-size systems, which is consistentwith a near-critical behavior. The critical exponents obtained from the finite-size spin and dimer correlationscould be compatible with the deconfined criticality in this small region. We compare and contrast our resultswith earlier numerical studies. PACS numbers: 73.43.Nq, 75.10.Jm, 75.10.Kt
Introduction.—
Quantum spin liquid (SL) is an exotic stateof matter where a spin system does not form magneti-cally ordered state or break lattice symmetries even at zerotemperature[1]. Understanding spin liquids is important infrustrated magnetic systems and may also hold clues to un-derstanding the non-Fermi liquid of doped Mott materials andhigh- T c superconductivity[2]. While the exciting propertiesof SL such as deconfined quasiparticles and fractional statis-tics have been revealed in many artificially constructed sys-tems [3–12], the possibility of finding spin liquids in realis-tic Heisenberg models has attracted much attention over thepast 20 years due to its close relation to experimental materi-als. The prominent example is the kagome antiferromagnet,where recent density matrix renormalization group (DMRG)studies point to a gapped Z SL[10, 13–16] characterized bya Z topological order and fractionalized spinon and visonexcitations[17–21].One of the candidate models for SL is the spin- J - J square Heisenberg model (SHM) with the Hamiltonian H = J (cid:88) (cid:104) i,j (cid:105) S i · S j + J (cid:88) (cid:104)(cid:104) i,j (cid:105)(cid:105) S i · S j , (1)where the sums (cid:104) i, j (cid:105) and (cid:104)(cid:104) i, j (cid:105)(cid:105) run over all the nearest-neighbor (NN) and the next-nearest-neighbor bonds, respec-tively. We set J = 1 . The frustrating J couplings suppressthe N´eel order and induce a nonmagnetic region around thestrongest frustration point J = 0 . [22–47]. Different candi-date states have been proposed based on approximate methodsor small-size exact diagonalization calculations, such as pla-quette valence-bond (PVB) state[26, 29, 32, 33, 35, 38, 46],the columnar valence-bond (CVB) state[24, 25, 28], or a gap-less SL[30, 31, 44, 45]. However, the true nature of the non-magnetic phase remains unresolved. J m s m s Δ Τ Δ Τ Neel
PVB gaplessregion J J J pin RC4-6 cylinder
FIG. 1: (color online) Phase diagram of spin- J - J SHM obtainedby our SU (2) DMRG studies. With growing J , the model has aN´eel phase for J < . and a PVB phase for . < J < . .Between these two phases, the finite-size magnetization and spin gapappear small in our calculations, consistent with a near-critical be-havior. The main panel shows N´eel order parameter m s and spin gap ∆ T in the thermodynamic limit. The inset is a sketch of a RC4-6cylinder; J pin shows the modified odd vertical bonds providing theboundary pinning for dimer orders. Recent DMRG study of the J - J SHM [40] proposed agapped Z SL for . ≤ J ≤ . by establishing the ab-sence of the magnetic and dimer orders, and by measuring apositive topological entanglement entropy term close to thevalue ln 2 expected for a Z SL[48, 49]. Very recent varia-tional Monte Carlo (VMC) work[45] proposed a gapless Z SL for . (cid:46) J (cid:46) . . On the other hand, recent DMRGstudies[50–52] of another bipartite frustrated system—the J - J spin- / honeycomb Heisenberg model—found a PVB a r X i v : . [ c ond - m a t . s t r- e l ] N ov phase in the nonmagnetic region, with a possible SL phasebetween the N´eel and PVB phases[52] or with a direct N´eel toPVB transition characterized by deconfined quantum critical-ity [50–54]. These studies[51, 52] also found that in the non-magnetic region the convergence of DMRG in wider systems,which is controlled by the number of states kept, is crucial fordetermining the true nature of the ground state.In this Letter, we reexamine the nonmagnetic region ofthe J - J SHM using DMRG with SU (2) spin rotationsymmetry[55]. We obtain accurate results on wide cylindersby keeping as many as U (1) -equivalent states. We finda N´eel phase below J (cid:39) . and a nonmagnetic region for . < J < . by finite-size scaling of the magnetic or-der parameter. In the nonmagnetic region, we establish a PVBorder for J > . —in contrast to the previous proposal of agapped Z SL[40]—by observing that the PVB decay lengthgrows strongly with increasing system width. We identify thePVB order as the s -wave plaquette[33] by studying dimer-dimer correlations. For . < J < . , we find that themagnetic order, valence-bond crystal (VBC) orders, and spinexcitation gap are small on finite-size systems, suggesting anear-critical behavior. The magnetic and dimer critical ex-ponents at J = 0 . are roughly similar to the values foundfor the deconfined criticality in the J - Q models on the squareand honeycomb lattices[56–63], which is consistent with thedeconfined criticality scenario conjectured also for the J - J model in Ref. [64].We establish the phases based on high accuracy DMRG re-sults on cylinders[65]. The first cylinder is the rectangularcylinder (RC) with closed boundary in the y direction andopen boundaries in the x direction. We denote it as RC L y - L x , where L y and L x are the number of sites in the y and x directions; the width of the cylinder is W y = L y (see the insetof Fig. 1). To study the dimers oriented in the y direction, wecan induce such an order near the open boundaries by mod-ifying every other NN vertical bond on the boundary to be J pin (cid:54) = J as illustrated in Fig. 1. The second geometry is thetilted cylinder (TC), as shown in Fig. 4(a), when discussingVBC order. N´eel order.—
The N´eel order parameter m s is defined as m s = N (cid:80) i,j (cid:104) S i · S j (cid:105) e i(cid:126)q · ( (cid:126)r i − (cid:126)r j ) ( N is the total site number),with (cid:126)q = ( π, π ) . We calculate m s from the spin correlationsof the L × L sites in the middle of the RC L - L cylinder, whichefficiently reduces boundary effects[40, 66]. In Fig. 2(a), weshow m s for different systems with L = 4 − [67]. We showthe obtained two-dimensional (2D) limit m s, ∞ in the inset ofFig. 2(a). Such an analysis suggests that the N´eel order van-ishes for J > . .The estimated J of spin order vanishing is different fromthe point J = 0 . where the PVB order develops as found be-low. One possibility is an intermediate SL phase[44, 45]. An-other possibility is that the system is near critical for .
We introduce the “pinning” bonds J pin (cid:54) = J on boundaries to induce a vertical dimer pattern, and measurethe decay length of the dimer order parameter (DOP) texturefrom the edge to the middle of the cylinder[40, 64]. The ver-tical DOP (vDOP) is defined as the difference between thestrong and weak vertical bond energies. In Fig. 3(a), we showa log-linear plot of the vDOP for J = 0 . and J pin = 2 . on long cylinders. We find that, although the amplitude ofthe vDOP texture changes with J pin , the decay length ξ y isindependent of J pin (see Supplemental Material[68]). In theinset of Fig. 3(a) we compare our ξ y with those in Ref. [40]. d v DO P RC4-48RC6-48RC8-48RC10-48RC12-48 4 6 8 10 12 W y ξ y our resultRef. [40] J =0.5 (a) J pin =2.0 W y ξ y J =0.35J =0.4J =0.45J =0.48J =0.5J =0.52J =0.55 (b) W y ξ x (c) FIG. 3: (color online) (a) Log-linear plot of vDOP for J = 0 . and J pin = 2 . on the RC cylinder. The inset is the comparison of widthdependence of the vertical dimer decay length ξ y with Ref. [40]. (b),(c) ξ y and ξ x versus W y on RC cylinders with J pin = 2 . for a rangeof J shown with the same symbols in both panels. We observe consistency for W y ≤ , but disagreement for W y ≥ [69]. The disagreement might originate from lessgood convergence in Ref. [40]. Our results are fully convergedby keeping 16000 (24000) states for L y = 10 ( ) systems. InFig. 3(b), we show the width dependence of ξ y for various J . ξ y grows slowly and saturates on wide cylinders for J < . ,demonstrating the vanishing VBC order. For J > . , ξ y - . - . - . - . - . - . - . - . - . - . . . - . - . . . - . - . . . - . - . . . . - . - . . . - . - . . . - . - . . . - . . . - . - . . . - . - . . . - . - . . . - . - . - . . . - . - . . . - . - . . . - . - . . - . - . . . - . - . . . - . - . . . - . - . . . . - . - . . . - . - . . . - . - . . . - . . . - . - . . . - . - . . . - . - . . . - . - . - . . . - . - . . . - . - . . . - . - . . - . - . . . - . - . . . - . - . . . - . - . . . . - . - . . . - . - . . . - . - . . . - . . . - . - . . . - . - . . . - . - . . . - . - . - . . . - . - . . . - . - . . . - . - . . - . - . . . - . - . . . - . - . . . - . - . . . . - . - . . . - . - . . . - . - . . . - . . . - . - . . . - . - . . . - . - . . . - . - . - . . . - . - . . . - . - . . . - . - . . - . - . . . - . - . . . - . - . . . - . - . . . . - . - . . . - . - . . . - . - . . . - . . . - . - . . . - . - . . . - . - . . . - . - . - . . . - . - . . . - . - . . . - . - . . - . - . . . - . - . . . - . - . . . - . - . . . . - . - . . . - . - . . . - . - . . . - . . . - . - . . . - . - . . . - . - . . . - . - . - . . . - . - . . . - . - . . . - . - . . (a) TC8-25, J =0.55, left half lattice E s E s E s E w E w W y ξ P J =0.35J =0.4J =0.45J =0.48J =0.5J =0.52J =0.55 (b) FIG. 4: (color online) (a) The NN J bond energy for J = 0 . on the left half of the tilted TC8-25 cylinder, where we have sub-tracted a constant − . from all the bond energies. Note that theTC cylinders have the square lattice rotated by 45 deg compared toFig. 1. We trim every other site on both boundaries to make latticeselect unique PVB order. E s ( E w ) denotes the sum of four NN bondenergies of the red (blue) plaquette with negative (positive) numbers.(b) Dependence of the pDOP decay length ξ P on the cylinder width W y . grows faster than linear, suggesting nonzero vDOP in the 2Dlimit. In addition to the vertical dimer, the system also has thehorizontal bond dimer with an exponentially decaying hori-zontal DOP (hDOP). In Fig. 3(c), we show that the hDOPdecay length ξ x also grows strongly for J > . . The coex-isting nonzero horizontal and vertical dimer orders suggest aPVB state.We also study the dimer structure factors S VBC and S col defined in Ref. [33]; the former detects both the PVB andCVB orders while the latter is nonzero only for the CVB or-der. We take RC L - L cylinders without pinning and calculatethe structure factors using the dimer correlations of the middle L × L sites. The picture of the dimer correlations is consistentwith the s -wave plaquette state[33]. The finite-size extrapo-lations show that while S VBC /N possibly approaches finitevalues for J > . , S col /N clearly approaches zero, which Δ Τ J =0.4, RCJ =0.46, RCJ =0.5, RCJ =0.55, RC E S=0,M -E S=0, ∞ ε E S=1,M -E S=0, ∞ (c)(b)(a) J =0.5, RC10-20J =0.5, RC10-20 FIG. 5: (color online) (a), (b) Ground-state energies for RC10-20cylinder at J = 0 . in the S = 0 ( E S =0 ,M ) and S = 1 ( E S =1 ,M )sectors versus the DMRG truncation error ε . All the energies havesubtracted the ground-state energy E S =0 , ∞ = − . . M isthe number of kept U (1) -equivalent DMRG states and is indicatednext to the symbols. (c) Finite-size extrapolations of spin gap ∆ T on RC L - L cylinders ( L = 4 , , , ). For J < . , the data arefitted using the formula ∆ T ( L ) = ∆ T ( ∞ )+ α/L + β/L + γ/L ,while for J ≥ . , we fit the data using ∆ T ( L ) = ∆ T ( ∞ )+ a/L + b/L + c/L . We estimate ∆ T ( ∞ ) = 0 . ± . and . ± . for J = 0 . and . , respectively. definitely excludes the CVB order.To explicitly demonstrate PVB order, we study the TC ob-tained by cutting the cylinder along the J direction and trim-ming every other site on the boundary as shown in Fig. 4(a).We label it as TC L y - L x , where L y and L x denote the numberof square plaquettes along the y and x directions; the widthof the cylinder is W y = √ L y . The trimmed edges inducestrong PVB order on the boundaries. We denote the sum ofthe four NN bond energies of a “strong” red (“weak” blue)plaquette as E s ( E w ), and define the plaquette DOP (pDOP)as the difference between E s and E w , which is found to de-cay exponentially with a decay length ξ P . In Fig. 4(b), we findstrong growth of ξ P with W y for J > . , consistent with aPVB state. By studying the log-log plot of VBC order pa-rameter versus system size (see Supplemental Material[68]),we estimate the anomalous exponent of dimer correlations η VBC (cid:39) . at J = 0 . , which is not far from estimates in thedeconfined criticality scenario in the J - Q models on square( η VBC (cid:39) . )[56–62] and honeycomb ( η VBC (cid:39) . )[63]lattices. Spin gap and ground state energy.—
We calculate the spingap ∆ T on the RC L - L cylinders up to L = 10 followingthe method from Ref. [14]: We sweep the ground state first,and then target the S = 1 sector sweeping the middle L × L sites to avoid edge excitations. In Figs. 5(a) and 5(b), we showenergies versus the DMRG truncation error for the RC10-20cylinder at J = 0 . in the S = 0 and S = 1 sectors. In bothplots we have subtracted the ground-state energy − . obtained by the extrapolation in Fig. 5(a). We find that weneed about twice as many states to achieve the same energy error in the S = 1 sector as in the S = 0 sector. The difficultyto reach the convergent energy in the S = 1 sector may ex-plain the overestimate of the spin gap in the earlier work [40]:We find ∆ T (cid:39) . while Ref. [40] estimates ∆ T (cid:39) . .We obtain accurate spin gaps by keeping up to states at L y = 10 , which sets the limit of our simulations.Figure 5(c) shows the finite-size extrapolations of ∆ T . Inour fits, we find ∆ T extrapolating vanishing for J ≤ . ,consistent with the N´eel order for J ≤ . . For J = 0 . and . , ∆ T ( L → ∞ ) is fitted to . ± . and . ± . , respectively; this is compatible with a VBC phase. Summary and discussion.—
We have studied the groundstate of spin- J - J SHM by accurate SU (2) DMRG sim-ulations. We find a N´eel order persisting up to J = 0 . .Contrary to the previous proposals of gapped Z SL fromDMRG[40] or gapless Z SL from VMC calculations[45],we establish an s -wave PVB state for J > . by observ-ing rapidly growing characteristic lengths of both the verticaland horizontal dimer orders on different cylinders. Betweenthe N´eel and PVB phases, we find a near-critical region thatcould be compatible with the deconfined criticality scenario.However, since the system in this region has large correlationlength scales that can be comparable to or even larger than thesystem widths we can approach, we cannot exclude a possi-ble gapless SL region proposed in variational studies[44, 45].We hope that future studies on larger system size, either push-ing DMRG further or using new techniques such as tensornetwork will be able to resolve between these scenarios moreclearly.We would like to particularly thank H.-C. Jiang and L.Balents for extensive discussions. We also acknowledgestimulating discussions with K. S. D. Beach, Z.-C. Gu,W.-J. Hu, L. Wang, S. White, Z.-Y. Zhu, and A. Sand-vik. This research is supported by the National ScienceFoundation through grants DMR-1205734 (S.S.G.), DMR-0906816 (D.N.S.), DMR-1206096 (O.I.M.), DMR-1101912(M.P.A.F.), the U.S. Department of Energy, Office of Ba-sic Energy Sciences under grant No. DE-FG02-06ER46305(W.Z), and by the Caltech Institute of Quantum Informa-tion and Matter, an NSF Physics Frontiers Center with sup-port of the Gordon and Betty Moore Foundation (O.I.M. andM.P.A.F.). [1] L. Balents, Nature (London) , 199 (2010).[2] P. A. Lee, N. Nagaosa, and X. G. Wen, Rev. Mod. Phys. , 17(2006).[3] R. Moessner and S. L. Sondhi, Phys. Rev. Lett. , 1881 (2001).[4] C. Nayak and K. Shtengel, Phys. Rev. B , 064422 (2001).[5] T. Senthil and O. Motrunich, Phys. Rev. B , 205104 (2002).[6] L. Balents, M. P. A. Fisher, and S. M. Girvin, Phys. Rev. B ,224412 (2002).[7] D. N. Sheng and L. Balents, Phys. Rev. Lett. , 146805 (2005).[8] S. V. Isakov, Y. B. Kim, and A. Paramekanti, Phys. Rev. Lett. , 207204 (2006). [9] S. V. Isakov, M. B. Hastings, and R. G. Melko, Nat. Phys. ,772 (2011).[10] Y. C. He, D. N. Sheng, and Y. Chen, Phys. Rev. B , 075110(2014).[11] A. Kitaev, Ann. Phys. (Amsterdam) , 2 (2006).[12] H. Yao and S. A. Kivelson, Phys. Rev. Lett. , 247206 (2012).[13] H. C. Jiang, Z. Y. Weng, and D. N. Sheng, Phys. Rev. Lett. ,117203 (2008).[14] S. Yan, D. Huse, and S. R. White, Science , 1173 (2011).[15] S. Depenbrock, I. P. McCulloch, and U. Schollw¨ock, Phys. Rev.Lett. , 067201 (2012).[16] H. C. Jiang, Z. H. Wang, and L. Balents, Nat. Phys. , 902(2012).[17] N. Read and S. Sachdev, Phys. Rev. Lett. , 1773 (1991).[18] X. G. Wen, Phys. Rev. B , 2664 (1991).[19] X. G. Wen, Phys. Rev. B , 7387 (1989).[20] L. Balents, M. P. A. Fisher, and C. Nayak, Phys. Rev. B ,1654 (1999).[21] T. Senthil and M.P.A. Fisher, Phys. Rev. B , 7850 (2000);Phys. Rev. Lett. , 292 (2001).[22] P. Chandra and B. Doucot, Phys. Rev. B , 9335(R) (1988).[23] L. B. Ioffe and A. I. Larkin, Int. J. Mod. Phys. B , 203 (1988).[24] S. Sachdev and R. N. Bhatt, Phys. Rev. B , 9323 (1990).[25] Andrey V. Chubukov and Th. Jolicoeur, Phys. Rev. B , 12050(1991).[26] M. E. Zhitomirsky and K. Ueda, Phys. Rev. B , 9007 (1996).[27] A. E. Trumper, L. O. Manuel, C. J. Gazza, and H. A. Ceccatto,Phys. Rev. Lett. , 2216 (1997).[28] R. R. P. Singh, Z. Weihong, C. J. Hamer, and J. Oitmaa, Phys.Rev. B , 7278 (1999).[29] L. Capriotti and S. Sorella, Phys. Rev. Lett. , 3173 (2000).[30] L. Capriotti, F. Becca, A. Parola, and S. Sorella, Phys. Rev. Lett. , 097201 (2001).[31] G. M. Zhang, H. Hu, and L. Yu, Phys. Rev. Lett. , 067201(2003).[32] K. Takano, Y. Kito, Y. Ono, and K. Sano, Phys. Rev. Lett. ,197202 (2003).[33] M. Mambrini, A. L¨auchli, D. Poilblanc, and F. Mila, Phys. Rev.B , 144422 (2006).[34] R. Darradi, O. Derzhko, R. Zinke, J. Schulenburg, S. E. Kr¨uger,and J. Richter, Phys. Rev. B , 214415 (2008).[35] L. Isaev, G. Ortiz, and J. Dukelsky, Phy. Rev. B , 024409(2009).[36] J. Richter and J. Schulenburg, Eur. Phys. J. B , 117 (2010).[37] J. Reuther and P. W¨olfle, Phys. Rev. B , 144410 (2010).[38] J. F. Yu, Y. J. Kao, Phys. Rev. B , 094407 (2012).[39] L. Wang, Z. C. Gu, F. Verstraete, and X. G. Wen,arXiv:1112.3331.[40] H. C. Jiang, H. Yao, and L. Balents, Phys. Rev. B , 024424(2012).[41] K. S. D. Beach, Phys. Rev. B , 224431 (2009).[42] F. Mezzacapo, Phys. Rev. B , 045115 (2012).[43] T. Li, F. Becca, W. J. Hu, and S. Sorella, Phys. Rev. B ,075111 (2012).[44] L. Wang, D. Poilblanc, Z. C. Gu, X. G. Wen, and F. Verstraete,Phys. Rev. Lett. , 037202 (2013).[45] W. J. Hu, F. Becca, A. Parola, and S. Sorella, Phys. Rev. B ,060402 (2013).[46] R. L. Doretto, Phys. Rev. B , 104415 (2014).[47] Y. Qi and Z. C. Gu, Phys. Rev. B , 235122 (2014).[48] A. Kitaev, J. Preskill, Phys. Rev. Lett. , 110404 (2006).[49] M. Levin, X.-G. Wen, Phys. Rev. Lett. , 110405 (2006).[50] R. Ganesh, Jeroen van den Brink, and S. Nishimoto, Phys. Rev.Lett. , 127203 (2013). [51] Zhenyue Zhu, D. A. Huse, and S. R. White, Phys. Rev. Lett. , 127205 (2013).[52] S. S. Gong, D. N. Sheng, Olexei I. Motrunich, and Matthew P.A. Fisher, Phys. Rev. B , 165138 (2013).[53] T. Senthil, L. Balents, S. Sachdev, A. Vishwanath, and MatthewP. A. Fisher, Phys. Rev. B , 144407 (2004).[54] T. Senthil, A. Vishwanath, L. Balents, S. Sachdev, and MatthewP. A. Fisher, Science , 1490 (2004).[55] I. P. McCulloch and M. Gul´acsi, Europhys. Lett. , 852 (2002);I. P. McCulloch, J. Stat. Mech. (2007) P10014.[56] A. W. Sandvik, Phys. Rev. Lett. , 227202 (2007).[57] R. G. Melko and R. K. Kaul, Phys. Rev. Lett. , 017203(2008).[58] F. Jiang, M. Nyfeler, S. Chandrasekharan, and U. Wiese, J. Stat.Mech. (2008) 02009.[59] J. Lou, A. W. Sandvik, and N. Kawashima, Phys. Rev. B ,180414(R) (2009).[60] J. Lou and A. W. Sandvik, Phys. Rev. B , 212406 (2009).[61] A. W. Sandvik, Phys. Rev. Lett. , 177201 (2010).[62] Matthew S. Block, R. G. Melko, and R. K. Kaul, Phys. Rev.Lett. , 137202 (2013).[63] S. Pujari, K. Damle, and F. Alet, Phys. Rev. Lett. , 087203(2013).[64] A. W. Sandvik, Phys. Rev. B , 134407 (2012).[65] We get much better convergence on cylinder than on torus. Wecan achieve the truncation error × − for L y = 10 cylinderand × − for L y = 12 cylinder. For example, for J = 0 . and L y = 10 , we get the truncation error × − on cylinderby keeping states, while the error is much larger × − on × torus even when we keep states.[66] S. R. White and A. L. Chernyshev, Phys. Rev. Lett. , 127004(2007).[67] We get converged m s for L ≤ by keeping more than states, while the results for L = 14 are obtained through ex-trapolations with the DMRG truncation error.[68] See Supplementary Material for more details.[69] In U (1) DMRG calculations for L y ≥ , one needs to extrap-olate ξ y with the DMRG truncation error, which may give someuncertainty to ξ y . 1 y -0.54-0.53-0.52-0.51-0.5-0.49 e J =0.5, torusJ =0.5, RCJ =0.5, RC, J pin =2J =0.5, TCJ =0.55, torusJ =0.55, RCJ =0.55, RC, J pin =2J =0.55, TC FIG. 6: (color online) DMRG ground-state energy per site for J =0 . and . on torus, on RC cylinder without pinning or with ver-tical dimer boundary pinning J pin = 2 . , as well as on TC cylin-der. The energies on torus are obtained through extrapolation withDMRG truncation error (see Table I). On cylinder, we get bulk en-ergy by subtracting the energies of two long cylinders with differentsystem lengths. With growing system width, the energies on differ-ent samples approach each other, giving the estimates of ground-stateenergy in the 2D limit for J = 0 . and . as e ∞ (cid:39) − . and − . , respectively. Supplementary Information
DMRG ground-state energies for J = 0 . and . .— Weshow our DMRG ground-state energies for J = 0 . and . in Fig. 6. We study L × L torus systems with L = 4 , , , .We keep more than M = 32000 optimal states for DMRGsweeping, and estimate the energy through extrapolation offinite- M energies via DMRG truncation error (see data in Ta-ble I below). For cylinders, we obtain bulk energy by subtract-ing the energies of two long cylinders with different systemlengths to eliminate boundary effects.As shown in Fig. 6, the energies per site of all samples in-crease slowly with increasing system width W y and approachclose to each other for W y (cid:38) . The energies on torus arelower than those on cylinder, and the difference decreases withincreasing W y . The bulk energy on RC cylinder is essentiallyindependent of the boundary pinning J pin . As the ground-state energy appears close to convergence for W y ≥ , wetake a simple straight line fitting of the large-size results togive estimates of the energy in the 2D limit e ∞ as shownby the dashed lines in Fig. 6. We find e ∞ (cid:39) − . and − . for J = 0 . and . , respectively. Horizontal dimer order on RC cylinder without pinning.—
On RC cylinder without pinning, the open edges break thelattice translation symmetry only in the x direction. The hori-zontal nearest-neighbor (NN) bond energies have the “strong-weak” dimer pattern as shown in Fig. 7(a). We define thehDOP as the difference of the adjacent horizontal NN bondenergies, which decays exponentially with a decay length ξ x .In Fig. 7(b), we show the hDOP decay length ξ x versus sys-tem width. For J < . , ξ x grows more slowly than linearly - . -0.0351 - . -0.0352 - . -0.0351 - . -0.0351 - . -0.0351 - . -0.0351 - . -0.0350 - . -0.0350 - . -0.0350-0.0350 - . +0.0386 - . +0.0386 - . +0.0387 - . +0.0388 - . +0.0389 - . +0.0390 - . +0.0391 - . +0.0391 - . +0.0391+0.0390 - . -0.0207 - . -0.0207 - . -0.0207 - . -0.0206 - . -0.0205 - . -0.0205 - . -0.0205 - . -0.0205 - . -0.0205-0.0206 - . +0.0172 - . +0.0172 - . +0.0173 - . +0.0173 - . +0.0174 - . +0.0174 - . +0.0173 - . +0.0173 - . +0.0173+0.0172 - . -0.0115 - . -0.0114 - . -0.0114 - . -0.0114 - . -0.0113 - . -0.0113 - . -0.0113 - . -0.0114 - . -0.0114-0.0115 - . +0.0088 - . +0.0089 - . +0.0090 - . +0.0090 - . +0.0090 - . +0.0090 - . +0.0090 - . +0.0090 - . +0.0089+0.0089 - . -0.0061 - . -0.0061 - . -0.0061 - . -0.0060 - . -0.0060 - . -0.0060 - . -0.0060 - . -0.0060 - . -0.0061-0.0061 - . +0.0047 - . +0.0047 - . +0.0048 - . +0.0048 - . +0.0048 - . +0.0048 - . +0.0048 - . +0.0048 - . +0.0047+0.0047 - . -0.0032 - . -0.0032 - . -0.0032 - . -0.0031 - . -0.0031 - . -0.0031 - . -0.0031 - . -0.0031 - . -0.0032-0.0032 - . +0.0025 - . +0.0025 - . +0.0025 - . +0.0025 - . +0.0026 - . +0.0026 - . +0.0026 - . +0.0025 - . +0.0025+0.0024 - . -0.0018 - . -0.0017 - . -0.0017 - . -0.0017 - . -0.0017 - . -0.0016 - . -0.0017 - . -0.0017 - . -0.0017-0.0018 - . - . - . - . - . - . - . - . - . (a) RC10-60, x=(1,12), J =0.55 W y ξ x J =0.4J =0.44J =0.48J =0.5J =0.55 (b) FIG. 7: (color online) (a) Subtracted NN bond energies for J =0 . on RC10-60 cylinder without pinning; the subtracted number − . is the average horizontal bond energy in the bulk of thelattice. Here we show the left columns. The alternation of red(negative number) and blue (positive number) bonds indicates hori-zontal dimer texture. (b) hDOP decay length ξ x versus system widthon RC cylinder without pinning. The extracted decay lengths ξ x aresimilar to those in Fig. 3(c) of the main text obtained on RC cylinderwith the vertical dimer pinning. and approaches saturation on large size, which is consistentwith vanishing dimer order. However, for J > . , ξ x growsfast, suggesting nonzero bulk hDOP on wider cylinders. Thishorizontal dimer order supports our claim of the VBC statefor J > . . We also find that the ξ x obtained here are al-most the same as those in Fig. 3(c) of the main text where thecylinder systems have the vertical bond pinning. Pinning independence of the vDOP decay length ξ y .— Inthe main text, we introduced modified vertical bonds J pin onboundaries to break the lattice translational symmetry in the y direction, allowing us to study the vDOP and the width d v DO P RC10-24RC10-30RC10-48RC10-60 ξ y ~3.531 J =0.5, J pin =0.0 FIG. 8: (color online) Log-linear plot of the vDOP on RC10 cylin-ders with different system lengths at J = 0 . and J pin = 0 . . Theexponential fitting gives decay length ξ y = 3 . . d vDOP J =0.5, J pin =0.0 ξ y ~1.81J =0.5, J pin =2.0 ξ y ~1.81 (a) RC6-48 d vDOP J =0.5, J pin =0.0 ξ y ~3.53J =0.5, J pin =2.0 ξ y ~3.53 (b) RC10-48 d vDOP J =0.55, J pin =0.0 ξ y ~3.12J =0.55, J pin =2.0 ξ y ~3.06 (c) RC8-48 d vDOP J =0.55, J pin =0.0 ξ y ~4.88J =0.55, J pin =2.0 ξ y ~4.88 (d) RC10-60 FIG. 9: (color online) Comparisons of the vDOP textures on RCcylinder with different boundary pinnings. We have studied severaldifferent J pin and found that although the vDOP varies with J pin , thedecay length ξ y is almost independent of the pinning strength. Herewe show the results with J pin = 0 . and . . dependence of the vDOP decay length. A direct question iswhether the pinning strength affects these quantities. We havecompared the vDOP and its decay length for several differentpinning strengths, from weak pinning J pin = 1 . , . , . to strong pinning J pin = 2 . and J pin = 0 . . First of all,in Fig. 8 we show that our results are obtained on quite longcylinders, thus minimizing the influence of finite-size effectson the decay length. Next, in Fig. 9 we show some exam-ples of varying boundary pinning at J = 0 . , . ; we findthat although the amplitude of the vDOP texture varies with J pin , the decay length ξ y is almost independent of the pinningstrength, indicating that our results with pinning are robustproperties of the bulk (infinitely long cylinder) phase. Horizontal dimer order on RC cylinder with odd L y .— On - . -0.4314 - . -0.4314 - . -0.4314 - . -0.4314-0.4314 - . -0.2328 - . -0.2328 - . -0.2328 - . -0.2328-0.2328 - . -0.4312 - . -0.4312 - . -0.4312 - . -0.4312-0.4312 - . -0.2329 - . -0.2329 - . -0.2329 - . -0.2329-0.2329 - . -0.4311 - . -0.4311 - . -0.4311 - . -0.4311-0.4311 - . -0.2330 - . -0.2330 - . -0.2330 - . -0.2330-0.2330 - . -0.4311 - . -0.4311 - . -0.4311 - . -0.4311-0.4311 - . -0.2330 - . -0.2330 - . -0.2330 - . -0.2330-0.2330 - . -0.4311 - . -0.4311 - . -0.4311 - . -0.4311-0.4311 - . -0.2331 - . -0.2331 - . -0.2331 - . -0.2331-0.2331 - . -0.4310 - . -0.4310 - . -0.4310 - . -0.4310-0.4310 - . - . - . - . (a) RC5-56, x=(17,28), J =0.5 y h DO P J =0.0J =0.4J =0.44J =0.48J =0.5J =0.52J =0.55 (b) FIG. 10: (color online) (a) The NN bond energies for RC5-56 cylin-der at J = 0 . , showing bonds with x from to . The systemhas a spontaneous bulk horizontal dimer order, and the bulk hDOPis defined as the difference of the strong and weak bond energies inthe middle of cylinder. (b) Width dependence of the hDOP on theodd- L y RC cylinders, showing the data for W y = 3 , , , . finite-size odd- L y RC cylinder, the system spontaneously de-velops a nonzero horizontal dimer order in the bulk, whichhappens both when the 2D phase is VBC or Z SL[12, 40].For a Z SL in the 2D limit, the dimer order would decay ex-ponentially with growing L y . On the other hand, for a VBCstate, it should go to a finite value in the 2D limit[12, 40]. Westudy the horizontal dimer order on odd- L y RC cylinder with L y up to and L x up to to get the results representing L x → ∞ cylinders. We define the absolute difference of thestrong and weak horizontal bond energies in the bulk as hDOP,see Fig. 10(a). We show thus measured hDOP versus /W y in Fig. 10(b). For J < . the hDOP decays fast with thecylinder width and appears to extrapolate to zero, while for J > . the hDOP has a slow decay and seems to saturate toa finite value. The nonzero hDOP does not support a Z SL,but indicates a VBC state for J > . . Dimer structure factors on RC cylinder.—
The CVB orderbreaks rotational symmetry, while the PVB order preserves it.Following Ref. [33], we consider two structure factors S vbc and S col obtained from the dimer-dimer correlations. S vbc diverges in both the CVB and PVB states, while S col diverges +++++ ++++ +++++ _____ ____ + + + + ++ + + + ++ ++ + + + ++ + + + + ____ (a) (b) FIG. 11: (color online) Phase factors ε λ ( k, l ) for (a) “vbc” and (b)“col” dimer structure factors defined in Eq. (2). The solid bondshave phase factor ε λ ( k, l ) = 1 , while dashed bonds have ε λ ( k, l ) = − . The ( k, l ) bonds that are nearest neighbors to the reference bond ( i, j ) (central black solid bond) are omitted in the calculation of thestructure factors. - . +0.0092 - . +0.0095 - . +0.0099 - . +0.0103 - . +0.0122 - . +0.0104 - . +0.0101 - . +0.0097 - . +0.0094+0.0094 + . -0.0153 + . -0.0162 + . -0.0178 + . -0.0192 + . -0.0243 + . -0.0194 + . -0.0180 + . -0.0164 + . -0.0155-0.0152 - . +0.0233 - . +0.0263 + . +0.0298 - . +0.0407 - . +0.0777 + . +0.0408 - . +0.0299 - . +0.0263 - . +0.0233+0.0234 + . -0.0352 + . -0.0425 + . -0.0531 + . -0.0758 + . + . -0.0760 + . -0.0534 + . -0.0428 + . -0.0355-0.0348 - . +0.0620 - . +0.0886 + . +0.1184+0.4059 + . +0.4060 - . +0.1184 - . +0.0885 - . +0.0619+0.0621 - . -0.0355 - . -0.0428 + . -0.0534-0.0760 + . -0.0758 - . -0.0531 - . -0.0425 - . -0.0351-0.0346 + . +0.0233 + . +0.0263 + . +0.0298 + . +0.0406 + . +0.0773 + . +0.0405 + . +0.0297 + . +0.0262 + . +0.0233+0.0234 - . -0.0153 - . -0.0162 + . -0.0179 - . -0.0192 - . -0.0242 - . -0.0191 - . -0.0177 - . -0.0160 - . -0.0152-0.0150 + . +0.0092 + . +0.0096 + . +0.0099 + . +0.0103 + . +0.0120 + . +0.0102 + . +0.0097 + . +0.0093 + . +0.0090+0.0091 - . -0.0059 - . -0.0060 - . -0.0062 - . -0.0063 - . -0.0065 - . -0.0062 - . -0.0060 - . -0.0058 - . -0.0057-0.0057 - . + . - . + . - . - . + . - . + . - . RC10-20, J =0.55 FIG. 12: (color online) Dimer-dimer correlation function for J =0 . on RC10-20 cylinder. The black bond in the middle of the cylin-der denotes the reference bond ( i, j ) . The blue and red bonds indi-cate the positive and negative correlations, respectively. Here themiddle × lattice dimer correlations are shown, which are usedto calculate the dimer structure factors. only in the CVB state. The structure factors are defined as S λ = (cid:88) ( k,l ) ε λ ( k, l ) C ijkl , (2)where λ is either “vbc” or “col”. The phase factors ε λ ( k, l ) are shown in Fig. 11, which reproduces Fig. 7 in Ref. [33].Dimer-dimer correlation function C ijkl is defined as C ijkl = 4 [ (cid:104) ( S i · S j )( S k · S l ) (cid:105) − (cid:104) S i · S j (cid:105)(cid:104) S k · S l (cid:105) ] . (3)We calculate dimer-dimer correlation function on the RC L - ε S vb c / N b J =0.55, L=8 (a) ε S vb c / N b J =0.55, L=12 (b) FIG. 13: (color online) Dimer structure factor S vbc /N b versus theDMRG truncation error ε for J = 0 . on (a) RC8-16 and (b)RC12-24 cylinders. For the RC12-24 cylinder at J = 0 . , wereach the truncation error × − by keeping up to U (1) equivalent states, and we extrapolate the data to estimate the re-sult for ε = 0 using the extrapolation function S vbc ( ε ) /N b = S vbc (0) /N b + aε + bε . The red lines indicate the extrapolationerror bar. L cylinder with a reference bond ( i, j ) in the middle ofthe cylinder (we have considered both horizontal and verti-cal reference bonds but will show only the former). Figure 12shows the dimer-dimer correlations on the RC10-20 cylinderat J = 0 . with the reference bond ( i, j ) oriented horizon-tally in the middle of the cylinder. The red and blue bonds in-dicate negative and positive dimer correlations, respectively.We see alternating red and blue horizontal bonds of compa-rable strengths, while the vertical bonds show significantlyweaker correlations; this picture looks much more like the pat-tern of the pure s -wave plaquette state (PVB) in Table III ofRef. [33] rather than the pattern of the pure columnar state.Figures 13 and 14 shows our best estimates of the structurefactors S vbc /N b and S col /N b obtained with a horizontal ref-erence bond ( i, j ) and normalized by the number of bonds N b used to calculate the structure factors. We note that we needto pay much attention to the DMRG convergence to obtain re-liable results. Figure 13 shows the structure factor S vbc /N b versus the DMRG truncation error ε at J = 0 . . FromFig. 13(a), we can obtain accurate results (within . ) for S vb c / N b J =0.42J =0.46J =0.48J =0.5J =0.52J =0.55 (a) S c o l / N b J =0.42J =0.46J =0.48J =0.5J =0.52J =0.55 (b) L S vb c / N b J =0.42J =0.46J =0.48J =0.5J =0.52J =0.55 (c) FIG. 14: (color online) Size dependence of dimer structure factors (a) S vbc /N b and (b) S col /N b . S vbc /N b appears to extrapolate to finitevalues for J > . , while S col /N b decays quite fast and approacheszero with increasing system size thus excluding the CVB order. (c)Log-log plot of S vbc /N b versus width L . The error bars for L = 12 data denote the extrapolation uncertainties with DMRG truncationerror as shown in Fig. 13(b). the RC8-16 cylinder by just using the data with the smallest ε . We can similarly obtain accurate results for the RC10-20cylinder (not shown) with ε (cid:39) × − . However, for theRC12-24 cylinder, our truncation error is still ε (cid:39) × − even when we keep as many as U (1) equivalent states.In this case, we extrapolate the data with ε to estimate the result for ε = 0 as shown in Fig. 13(b). The extrapolationfunction is S vbc ( ε ) /N b = S vbc (0) /N b + aε + bε . Withoutsuch an extrapolation, we would underestimate the magnitudeof the VBC order parameter by about . The error bar inFig. 13(b) is from the extrapolation uncertainty.In Fig. 14(a), we see that S vbc /N b approaches zero uponincreasing system size for J < . and possibly extrapolatesto finite values for J > . if we fit the large-size data us-ing polynomials of /N . This suggests PVB or CVB ordersat J > . . In Fig. 14(b), we see that S col /N b decays quitefast with system size and always approaches zero in the ther-modynamic limit, which implies vanishing CVB order. Thus,the behavior of these two structure factors reveals the possiblePVB order at J > . and clearly excludes the CVB order.We observe similar results with a vertical reference bond ( i, j )(not shown).We also notice that when we plot S vbc /N b versus /L ( L = √ N ), the data could be extrapolated to zero or smallvalues also for J > . , which would be similar to the anal-ysis in Ref. [40]. However, for J = 0 . , the extrapolationfunction to zero is almost linear in /L (plot not shown), whilein a phase with no VBC order we would expect S vbc /N b tovanish as /N ∼ /L . Thus this data is not consistent withvanishing VBC order.In Fig. 14(c), we show log-log plot of S vbc /N b versuscylinder width L , which provides more insight about the tran-sition region. For J < . , the accelerated decay of S vbc /N b is consistent with vanishing dimer order. The finite-size dataat J = 0 . and . appear close to critical, with power lawbehavior S vbc /N b ∼ L − . . On the other hand, at J = 0 . the S vbc /N b decays significantly more slowly, which is con-sistent with a VBC order. Overall, our structure factor re-sults are consistent with the phase diagram in the main textobtained from the studies of the VBC texture decay lengths. Dimer order of the next-nearest-neighbor bonds.—
We alsostudied the dimer order of the next-nearest-neighbor (NNN) J bonds by investigating the NNN bond energy textures andhow their decay length depends on the cylinder width.On the RC cylinder without additional boundary pinning,the diagonal NNN bond energy is the same inside each columnbut depends on the column distance from the boundary. Wedefine the NNN dimer order parameter as the difference ofthe NNN bond energies in adjacent columns. By studying theNNN dimer order parameter on long cylinders, we find thatit decays exponentially from the boundary to the bulk witha decay length ξ NNN . As shown in Fig. 15(a), ξ NNN growsfaster than linearly for J > . , while it saturates for J < . ; these results are consistent with the behavior of the decaylength of the NN bond texture.On the TC cylinder, the NNN bonds are either horizontalor vertical. We find that in both cases the textures have thesame decay length, which exhibits the same behavior withincreasing system width as the NN bond decay length, seeFig. 15(b). Thus, similar to the NN bonds, the NNN bondsalso appear to have a PVB order in the 2D limit. Such re-sults for the NNN bonds are expected and do not represent a0 W y ξ NNN J =0.46J =0.5J =0.55 RC cylinder (a) W y ξ NNN J =0.46J =0.5J =0.55 TC cylinder (b)
FIG. 15: (color online) (a) Dependence of the decay length of theNNN dimer order ξ NNN with system width W y on the RC cylinderwithout boundary pinning. (b) Dependence of ξ NNN with systemwidth W y on the TC cylinder. On this cylinder, we find that both thehorizontal and vertical dimer orders have the same ξ NNN . On bothcylinders, ξ NNN grows faster than linearly for J > . and saturatesfor J < . , consistent with the behavior of the decay length of theNN dimer orders. new VBC state but generally confirm the VBC order estab-lished from the NN bond energy studies. Indeed, in the PVBphase, we expect the NNN bonds on the strong plaquettes tohave somewhat different bond energy than the NNN bonds onthe weak plaquettes, and our results are consistent with suchexpectations. Comparisons of torus energies from DMRG and VMC.—
We have compared our DMRG ground-state energies on L × L tori to VMC results with additional Lanczos improvementsteps from Ref. [45]. Since the torus system is extremely dif-ficult to fully converge for × or larger sizes, we keep upto states and extrapolate the energy with the DMRGtruncation error[65]. The extrapolated DMRG and Lanczos-VMC results are quite close to each other in the possible SL ε -0.4995-0.499-0.4985-0.498-0.4975-0.497-0.4965-0.496 e J =0.5, L=8 (a) ε -0.499-0.498-0.497-0.496-0.495-0.494-0.493-0.492 e J =0.5, L=10 (b) FIG. 16: (color online) Ground-state energy per site e versus DMRGtruncation error ε on the L × L torus systems for (a) J = 0 . , L = 8 ,and (b) J = 0 . , L = 10 . The numbers in the figures denotethe kept SU (2) states M for obtaining the energy. We extrapolatethe data with a straight line fitting and denote the ε → intercept(corresponding to M → ∞ ) as DMRG ( ∞ ). region . < J < . , indicating that the gapless Z SL ofRef. [45] has very competitive energies in this region.Table I shows energy comparisons of DMRG and VMC for J = 0 . , . , . , and . ; it includes DMRG results ob-tained by keeping , , and SU (2) states [equiv-alent to about , , and U (1) optimal states],as well as VMC results with Lanczos improvement steps fromRef. [45]. DMRG ( ∞ ) denotes the DMRG energy extrapo-lated with truncation error; as illustrated in Fig. 16, we extrap-olate the data points using a straight line fitting. VMC ( p = ∞ )denotes the VMC energy extrapolations with the variance inRef. [45]. The overall agreement shows, on one hand, that theDMRG is performing reasonably well even in the most chal-lenging torus geometry. Here we emphasize that all results inthe main text are obtained using cylinder geometry where theDMRG measurements are much better converged[65] and rep-resent essentially exact unbiased measurements. On the otherhand, the excellent performance of the Lanczos-VMC methodis also notable. It would be interesting to see this method triedin the cylinder geometries and results subjected to the finite-size scaling analysis as in the present work.1 J = 0 . DMRG ( ) DMRG ( ) DMRG ( ) DMRG ( ∞ ) VMC ( p = 0 ) VMC ( p = 1 ) VMC ( p = 2 ) VMC ( p = ∞ ) L = 6 − . − . − . − . − . − . − . − . L = 8 − . − . − . − . − . − . − . − . L = 10 − . − . − . − . − . − . − . − . J = 0 . DMRG ( ) DMRG ( ) DMRG ( ) DMRG ( ∞ ) VMC ( p = 0 ) VMC ( p = 1 ) VMC ( p = 2 ) VMC ( p = ∞ ) L = 6 − . − . − . − . − . − . − . − . L = 8 − . − . − . − . − . − . − . − . L = 10 − . − . − . − . − . − . − . − . J = 0 . DMRG ( ) DMRG ( ) DMRG ( ) DMRG ( ∞ ) VMC ( p = 0 ) VMC ( p = 1 ) VMC ( p = 2 ) VMC ( p = ∞ ) L = 6 − . − . − . − . − . − . − . − . L = 8 − . − . − . − . − . − . − . − . L = 10 − . − . − . − . − . − . − . − . J = 0 . DMRG ( ) DMRG ( ) DMRG ( ) DMRG ( ∞ ) VMC ( p = 0 ) VMC ( p = 1 ) VMC ( p = 2 ) VMC ( p = ∞ ) L = 6 − . − . − . − . − . − . − . − . L = 8 − . − . − . − . − . − . − . − . L = 10 − . − . − . − . − . − . − . − . TABLE I: DMRG and VMC ground-state energies on L × L tori with J = 0 . , . , . and . . DMRG energies are obtained by keeping , , and SU (2) states. DMRG ( ∞ ) is obtained from the straight line energy extrapolation with DMRG truncation error asillustrated in Fig. 16. For L = 10 , we do not show the error bar of extrapolation because the DMRG truncation error is large ( (cid:15) (cid:39) × − )and thus our estimation of the error bar is not accurate. The extrapolated results for L = 10 are obtained from the linear fitting of the databy keeping , and SU (2) states. The VMC energies are from Ref. [45]; p denotes the Lanczos step; and VMC ( p = ∞ ) isobtained from extrapolation with the variance. W y -1012345 e n t r opy J =0.46, RCJ =0.55, RCJ =0.46, TCJ =0.55, TC FIG. 17: (color online) Entanglement entropy as a function of systemwidth on RC and TC cylinders. For each width, we obtain the entropyby extrapolating measurements on long cylinders to L x → ∞ limit. Entanglement entropy.—
For gapped quantum stateswith topological order, topological entanglement entropy(TEE) γ is proposed to characterize non-local feature ofentanglement[48, 49]. The Renyi entropies of a subsys-tem A with reduced density matrix ρ A are defined as S n = (1 − n ) − ln(Tr ρ nA ) ; n → limit gives the VonNeuman entropy. For a topologically ordered state, Renyientropies have the form S n = αL − γ , where L is theboundary of the subsystem, and all other terms vanish in the large L limit; α is a non-universal constant, while a positive γ is a correction to the area law of entanglement and reaches auniversal value determined by total quantum dimension D ofquasiparticle excitations[48, 49]. Previous DMRG study[40]found γ ≈ ln 2 in the intermediate region of J consistentwith a Z SL for this model. We compute the entanglemententropy (EE) on long cylinders by partitioning the system inthe middle along the vertical direction. For each fixed L y ,we fit the entropy to L x → ∞ limit to find the entropy of apossible minimum entropy state[16].In Fig. 17, we show our DMRG results for the EE at J = 0 . and . on both TC and RC cylinders. We obtainaccurate EE when W y < . For W y = 12 , we extrapolatethe EE with the DMRG truncation error, which has signifi-cant uncertainty from the extrapolation. On RC cylinder, weperform linear fit of the EE versus W y using the three largestsizes. We find the TEE at J = 0 . is close to ln 2 , whileat J = 0 . is close to − . . However, the system appearsto have large finite-size effects, which can be seen by com-paring the results on the RC and TC cylinders. On the TCcylinder, the linear fits of the EE vs W y give the TEE close tozero, which is different from the RC cylinder. Similar effecthas also been observed in the J - J model on the honeycomblattice[51, 52]. Because of such strong finite-size effects, theTEE obtained by fitting EE on our small sizes may not be ableto distinguish different quantum phases in the J - J2