# Plaquette Ordered Phase and Quantum Phase Diagram in the Spin-1/2 J1-J2 Square Heisenberg Model

## Abstract

We study the spin- Heisenberg model on the square lattice with first- and second-neighbor antiferromagnetic interactions and , which possesses a nonmagnetic region that has been debated for many years and might realize the interesting spin liquid. We use the density matrix renormalization group approach with explicit implementation of spin rotation symmetry and study the model accurately on open cylinders with different boundary conditions. With increasing , we find a Néel 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 the Néel order vanishes at . For , we find dimer correlations and PVB textures whose decay lengths grow strongly with increasing system width, consistent with a long-range PVB order in the two-dimensional limit. The dimer-dimer correlations reveal the -wave character of the PVB order. For , spin order, dimer order, and spin gap are small on finite-size systems, which is consistent with a near-critical behavior. The critical exponents obtained from the finite-size spin and dimer correlations could be compatible with the deconfined criticality in this small region. We compare and contrast our results with earlier numerical studies.

###### pacs:

73.43.Nq, 75.10.Jm, 75.10.KtIntroduction.—Quantum spin liquid (SL) is an exotic state of matter where a spin system does not form magnetically ordered state or break lattice symmetries even at zero temperature(1). Understanding spin liquids is important in frustrated magnetic systems and may also hold clues to understanding the non-Fermi liquid of doped Mott materials and high- superconductivity(2). While the exciting properties of SL such as deconfined quasiparticles and fractional statistics have been revealed in many artificially constructed systems (3); (4); (5); (6); (7); (8); (9); (10); (11); (12), the possibility of finding spin liquids in realistic Heisenberg models has attracted much attention over the past 20 years due to its close relation to experimental materials. The prominent example is the kagome antiferromagnet, where recent density matrix renormalization group (DMRG) studies point to a gapped SL(13); (14); (15); (16); (10) characterized by a topological order and fractionalized spinon and vison excitations(17); (18); (19); (20); (21).

One of the candidate models for SL is the spin- - square Heisenberg model (SHM) with the Hamiltonian

(1) |

where the sums and run over all the nearest-neighbor (NN) and the next-nearest-neighbor bonds, respectively. We set . The frustrating couplings suppress the Néel order and induce a nonmagnetic region around the strongest frustration point (22); (23); (24); (25); (26); (27); (28); (29); (30); (31); (32); (33); (34); (35); (36); (37); (38); (39); (40); (41); (42); (43); (44); (45); (46); (47). Different candidate states have been proposed based on approximate methods or small-size exact diagonalization calculations, such as plaquette valence-bond (PVB) state(26); (29); (32); (33); (35); (38); (46), the columnar valence-bond (CVB) state(24); (25); (28), or a gapless SL(30); (31); (44); (45). However, the true nature of the nonmagnetic phase remains unresolved.

Recent DMRG study of the - SHM (40) proposed a gapped SL for by establishing the absence of the magnetic and dimer orders, and by measuring a positive topological entanglement entropy term close to the value expected for a SL(48); (49). Very recent variational Monte Carlo (VMC) work(45) proposed a gapless SL for . On the other hand, recent DMRG studies(50); (51); (52) of another bipartite frustrated system—the - spin- honeycomb Heisenberg model—found a PVB phase in the nonmagnetic region, with a possible SL phase between the Néel and PVB phases(52) or with a direct Néel to PVB transition characterized by deconfined quantum criticality (50); (51); (52); (53); (54). These studies(51); (52) also found that in the nonmagnetic region the convergence of DMRG in wider systems, which is controlled by the number of states kept, is crucial for determining the true nature of the ground state.

In this Letter, we reexamine the nonmagnetic region of the - SHM using DMRG with spin rotation symmetry(55). We obtain accurate results on wide cylinders by keeping as many as -equivalent states. We find a Néel phase below and a nonmagnetic region for by finite-size scaling of the magnetic order parameter. In the nonmagnetic region, we establish a PVB order for —in contrast to the previous proposal of a gapped SL(40)—by observing that the PVB decay length grows strongly with increasing system width. We identify the PVB order as the -wave plaquette(33) by studying dimer-dimer correlations. For , we find that the magnetic order, valence-bond crystal (VBC) orders, and spin excitation gap are small on finite-size systems, suggesting a near-critical behavior. The magnetic and dimer critical exponents at are roughly similar to the values found for the deconfined criticality in the - models on the square and honeycomb lattices(56); (57); (58); (59); (60); (61); (62); (63), which is consistent with the deconfined criticality scenario conjectured also for the - model in Ref. (64).

We establish the phases based on high accuracy DMRG results on cylinders(65). The first cylinder is the rectangular cylinder (RC) with closed boundary in the direction and open boundaries in the direction. We denote it as RC-, where and are the number of sites in the and directions; the width of the cylinder is (see the inset of Fig. 1). To study the dimers oriented in the direction, we can induce such an order near the open boundaries by modifying every other NN vertical bond on the boundary to be as illustrated in Fig. 1. The second geometry is the tilted cylinder (TC), as shown in Fig. 4(a), when discussing VBC order.

Néel order.—The Néel order parameter is defined as ( is the total site number), with . We calculate from the spin correlations of the sites in the middle of the RC- cylinder, which efficiently reduces boundary effects(40); (66). In Fig. 2(a), we show for different systems with (67). We show the obtained two-dimensional (2D) limit in the inset of Fig. 2(a). Such an analysis suggests that the Néel order vanishes for .

The estimated of spin order vanishing is different from the point where the PVB order develops as found below. One possibility is an intermediate SL phase(44); (45). Another possibility is that the system is near critical for . In this case, to get some idea about the criticality, Fig. 2(b) shows the log-log plot of . approaches finite value in the Néel phase as seen for and . On the other hand, we expect at a critical point and in the nonmagnetic phase. The accelerated decay of at is consistent with vanishing Néel order: from the two largest sizes we estimate , which is quite close to . In the near-critical region, we fit the data to and the data () to . This range of is compatible with the findings in the - models on the square ()(56); (57); (58); (59); (60); (61); (62) and honeycomb ()(63) lattices, which show continuous Néel-to-VBC transition argued to be in the deconfined criticality class, so our model is compatible with this scenario as well.

VBC orders.—We introduce the “pinning” bonds on boundaries to induce a vertical dimer pattern, and measure the decay length of the dimer order parameter (DOP) texture from the edge to the middle of the cylinder(40); (64). The vertical DOP (vDOP) is defined as the difference between the strong and weak vertical bond energies. In Fig. 3(a), we show a log-linear plot of the vDOP for and on long cylinders. We find that, although the amplitude of the vDOP texture changes with , the decay length is independent of (see Supplemental Material(68)). In the inset of Fig. 3(a) we compare our with those in Ref. (40). We observe consistency for , but disagreement for (69). The disagreement might originate from less good convergence in Ref. (40). Our results are fully converged by keeping 16000 (24000) states for () systems. In Fig. 3(b), we show the width dependence of for various . grows slowly and saturates on wide cylinders for , demonstrating the vanishing VBC order. For , grows faster than linear, suggesting nonzero vDOP in the 2D limit. In addition to the vertical dimer, the system also has the horizontal bond dimer with an exponentially decaying horizontal DOP (hDOP). In Fig. 3(c), we show that the hDOP decay length also grows strongly for . The coexisting nonzero horizontal and vertical dimer orders suggest a PVB state.

We also study the dimer structure factors and defined in Ref. (33); the former detects both the PVB and CVB orders while the latter is nonzero only for the CVB order. We take RC- cylinders without pinning and calculate the structure factors using the dimer correlations of the middle sites. The picture of the dimer correlations is consistent with the -wave plaquette state(33). The finite-size extrapolations show that while possibly approaches finite values for , clearly approaches zero, which definitely excludes the CVB order.

To explicitly demonstrate PVB order, we study the TC obtained by cutting the cylinder along the direction and trimming every other site on the boundary as shown in Fig. 4(a). We label it as TC-, where and denote the number of square plaquettes along the and directions; the width of the cylinder is . The trimmed edges induce strong PVB order on the boundaries. We denote the sum of the four NN bond energies of a “strong” red (“weak” blue) plaquette as (), and define the plaquette DOP (pDOP) as the difference between and , which is found to decay exponentially with a decay length . In Fig. 4(b), we find strong growth of with for , consistent with a PVB state. By studying the log-log plot of VBC order parameter versus system size (see Supplemental Material(68)), we estimate the anomalous exponent of dimer correlations at , which is not far from estimates in the deconfined criticality scenario in the - models on square ()(56); (57); (58); (59); (60); (61); (62) and honeycomb ()(63) lattices.

Spin gap and ground state energy.—We calculate the spin gap on the RC- cylinders up to following the method from Ref. (14): We sweep the ground state first, and then target the sector sweeping the middle sites to avoid edge excitations. In Figs. 5(a) and 5(b), we show energies versus the DMRG truncation error for the RC10-20 cylinder at in the and sectors. In both plots we have subtracted the ground-state energy obtained by the extrapolation in Fig. 5(a). We find that we need about twice as many states to achieve the same energy error in the sector as in the sector. The difficulty to reach the convergent energy in the sector may explain the overestimate of the spin gap in the earlier work (40): We find while Ref. (40) estimates . We obtain accurate spin gaps by keeping up to states at , which sets the limit of our simulations.

Figure 5(c) shows the finite-size extrapolations of . In our fits, we find extrapolating vanishing for , consistent with the Néel order for . For and , is fitted to and , respectively; this is compatible with a VBC phase.

Summary and discussion.—We have studied the ground state of spin- - SHM by accurate DMRG simulations. We find a Néel order persisting up to . Contrary to the previous proposals of gapped SL from DMRG(40) or gapless SL from VMC calculations(45), we establish an -wave PVB state for by observing rapidly growing characteristic lengths of both the vertical and horizontal dimer orders on different cylinders. Between the Néel and PVB phases, we find a near-critical region that could be compatible with the deconfined criticality scenario. However, since the system in this region has large correlation length scales that can be comparable to or even larger than the system widths we can approach, we cannot exclude a possible gapless SL region proposed in variational studies(44); (45). We hope that future studies on larger system size, either pushing DMRG further or using new techniques such as tensor network will be able to resolve between these scenarios more clearly.

We would like to particularly thank H.-C. Jiang and L. Balents for extensive discussions. We also acknowledge stimulating discussions with K. S. D. Beach, Z.-C. Gu, W.-J. Hu, L. Wang, S. White, Z.-Y. Zhu, and A. Sandvik. This research is supported by the National Science Foundation 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 Basic Energy Sciences under grant No. DE-FG02-06ER46305 (W.Z), and by the Caltech Institute of Quantum Information and Matter, an NSF Physics Frontiers Center with support of the Gordon and Betty Moore Foundation (O.I.M. and M.P.A.F.).

Supplementary Information

DMRG ground-state energies for and .—We show our DMRG ground-state energies for and in Fig. 6. We study torus systems with . We keep more than optimal states for DMRG sweeping, and estimate the energy through extrapolation of finite- energies via DMRG truncation error (see data in Table 1 below). For cylinders, we obtain bulk energy by subtracting the energies of two long cylinders with different system lengths to eliminate boundary effects.

As shown in Fig. 6, the energies per site of all samples increase slowly with increasing system width and approach close to each other for . The energies on torus are lower than those on cylinder, and the difference decreases with increasing . The bulk energy on RC cylinder is essentially independent of the boundary pinning . As the ground-state energy appears close to convergence for , we take a simple straight line fitting of the large-size results to give estimates of the energy in the 2D limit as shown by the dashed lines in Fig. 6. We find and for and , respectively.

Horizontal dimer order on RC cylinder without pinning.—On RC cylinder without pinning, the open edges break the lattice translation symmetry only in the direction. The horizontal nearest-neighbor (NN) bond energies have the “strong-weak” dimer pattern as shown in Fig. 7(a). We define the hDOP as the difference of the adjacent horizontal NN bond energies, which decays exponentially with a decay length . In Fig. 7(b), we show the hDOP decay length versus system width. For , grows more slowly than linearly and approaches saturation on large size, which is consistent with vanishing dimer order. However, for , grows fast, suggesting nonzero bulk hDOP on wider cylinders. This horizontal dimer order supports our claim of the VBC state for . We also find that the obtained here are almost the same as those in Fig. 3(c) of the main text where the cylinder systems have the vertical bond pinning.

Pinning independence of the vDOP decay length .— In the main text, we introduced modified vertical bonds on boundaries to break the lattice translational symmetry in the direction, allowing us to study the vDOP and the width dependence of the vDOP decay length. A direct question is whether the pinning strength affects these quantities. We have compared the vDOP and its decay length for several different pinning strengths, from weak pinning to strong pinning and . First of all, in Fig. 8 we show that our results are obtained on quite long cylinders, thus minimizing the influence of finite-size effects on the decay length. Next, in Fig. 9 we show some examples of varying boundary pinning at ; we find that although the amplitude of the vDOP texture varies with , the decay length is almost independent of the pinning strength, indicating that our results with pinning are robust properties of the bulk (infinitely long cylinder) phase.

Horizontal dimer order on RC cylinder with odd .—On finite-size odd- RC cylinder, the system spontaneously develops a nonzero horizontal dimer order in the bulk, which happens both when the 2D phase is VBC or SL(40); (12). For a SL in the 2D limit, the dimer order would decay exponentially with growing . On the other hand, for a VBC state, it should go to a finite value in the 2D limit(40); (12). We study the horizontal dimer order on odd- RC cylinder with up to and up to to get the results representing cylinders. We define the absolute difference of the strong and weak horizontal bond energies in the bulk as hDOP, see Fig. 10(a). We show thus measured hDOP versus in Fig. 10(b). For the hDOP decays fast with the cylinder width and appears to extrapolate to zero, while for the hDOP has a slow decay and seems to saturate to a finite value. The nonzero hDOP does not support a SL, but indicates a VBC state for .

Dimer structure factors on RC cylinder.— The CVB order breaks rotational symmetry, while the PVB order preserves it. Following Ref. (33), we consider two structure factors and obtained from the dimer-dimer correlations. diverges in both the CVB and PVB states, while diverges only in the CVB state. The structure factors are defined as

(2) |

where is either “vbc” or “col”. The phase factors are shown in Fig. 11, which reproduces Fig. 7 in Ref. (33). Dimer-dimer correlation function is defined as

(3) |

We calculate dimer-dimer correlation function on the RC- cylinder with a reference bond in the middle of the cylinder (we have considered both horizontal and vertical reference bonds but will show only the former). Figure 12 shows the dimer-dimer correlations on the RC10-20 cylinder at with the reference bond () oriented horizontally in the middle of the cylinder. The red and blue bonds indicate negative and positive dimer correlations, respectively. We see alternating red and blue horizontal bonds of comparable strengths, while the vertical bonds show significantly weaker correlations; this picture looks much more like the pattern of the pure -wave plaquette state (PVB) in Table III of Ref. (33) rather than the pattern of the pure columnar state.

Figures 13 and 14 shows our best estimates of the structure factors and obtained with a horizontal reference bond () and normalized by the number of bonds used to calculate the structure factors. We note that we need to pay much attention to the DMRG convergence to obtain reliable results. Figure 13 shows the structure factor versus the DMRG truncation error at . From Fig. 13(a), we can obtain accurate results (within ) for the RC8-16 cylinder by just using the data with the smallest . We can similarly obtain accurate results for the RC10-20 cylinder (not shown) with . However, for the RC12-24 cylinder, our truncation error is still even when we keep as many as equivalent states. In this case, we extrapolate the data with to estimate the result for as shown in Fig. 13(b). The extrapolation function is . Without such an extrapolation, we would underestimate the magnitude of the VBC order parameter by about . The error bar in Fig. 13(b) is from the extrapolation uncertainty.

In Fig. 14(a), we see that approaches zero upon increasing system size for and possibly extrapolates to finite values for if we fit the large-size data using polynomials of . This suggests PVB or CVB orders at . In Fig. 14(b), we see that decays quite fast with system size and always approaches zero in the thermodynamic limit, which implies vanishing CVB order. Thus, the behavior of these two structure factors reveals the possible PVB order at and clearly excludes the CVB order. We observe similar results with a vertical reference bond () (not shown).

We also notice that when we plot versus (), the data could be extrapolated to zero or small values also for , which would be similar to the analysis in Ref. (40). However, for , the extrapolation function to zero is almost linear in (plot not shown), while in a phase with no VBC order we would expect to vanish as . Thus this data is not consistent with vanishing VBC order.

In Fig. 14(c), we show log-log plot of versus cylinder width , which provides more insight about the transition region. For , the accelerated decay of is consistent with vanishing dimer order. The finite-size data at and appear close to critical, with power law behavior . On the other hand, at the decays significantly more slowly, which is consistent with a VBC order. Overall, our structure factor results are consistent with the phase diagram in the main text obtained from the studies of the VBC texture decay lengths.

Dimer order of the next-nearest-neighbor bonds.—We also studied the dimer order of the next-nearest-neighbor (NNN) bonds by investigating the NNN bond energy textures and how 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 column but depends on the column distance from the boundary. We define the NNN dimer order parameter as the difference of the NNN bond energies in adjacent columns. By studying the NNN dimer order parameter on long cylinders, we find that it decays exponentially from the boundary to the bulk with a decay length . As shown in Fig. 15(a), grows faster than linearly for , while it saturates for ; these results are consistent with the behavior of the decay length of the NN bond texture.

On the TC cylinder, the NNN bonds are either horizontal or vertical. We find that in both cases the textures have the same decay length, which exhibits the same behavior with increasing system width as the NN bond decay length, see Fig. 15(b). Thus, similar to the NN bonds, the NNN bonds also appear to have a PVB order in the 2D limit. Such results for the NNN bonds are expected and do not represent a new VBC state but generally confirm the VBC order established from the NN bond energy studies. Indeed, in the PVB phase, we expect the NNN bonds on the strong plaquettes to have somewhat different bond energy than the NNN bonds on the weak plaquettes, and our results are consistent with such expectations.

Comparisons of torus energies from DMRG and VMC.—We have compared our DMRG ground-state energies on tori to VMC results with additional Lanczos improvement steps from Ref. (45). Since the torus system is extremely difficult to fully converge for or larger sizes, we keep up to states and extrapolate the energy with the DMRG truncation error(65). The extrapolated DMRG and Lanczos-VMC results are quite close to each other in the possible SL region , indicating that the gapless SL of Ref. (45) has very competitive energies in this region.

Table 1 shows energy comparisons of DMRG and VMC for , and ; it includes DMRG results obtained by keeping , , and states [equivalent to about , , and optimal states], as well as VMC results with Lanczos improvement steps from Ref. (45). DMRG () denotes the DMRG energy extrapolated with truncation error; as illustrated in Fig. 16, we extrapolate the data points using a straight line fitting. VMC (=) denotes the VMC energy extrapolations with the variance in Ref. (45). The overall agreement shows, on one hand, that the DMRG is performing reasonably well even in the most challenging torus geometry. Here we emphasize that all results in the main text are obtained using cylinder geometry where the DMRG measurements are much better converged(65) and represent essentially exact unbiased measurements. On the other hand, the excellent performance of the Lanczos-VMC method is also notable. It would be interesting to see this method tried in the cylinder geometries and results subjected to the finite-size scaling analysis as in the present work.

DMRG () | DMRG () | DMRG () | DMRG () | VMC () | VMC () | VMC () | VMC () | |

DMRG () | DMRG () | DMRG () | DMRG () | VMC () | VMC () | VMC () | VMC () | |

DMRG () | DMRG () | DMRG () | DMRG () | VMC () | VMC () | VMC () | VMC () | |

DMRG () | DMRG () | DMRG () | DMRG () | VMC () | VMC () | VMC () | VMC () | |

Entanglement entropy.—For gapped quantum states with topological order, topological entanglement entropy (TEE) is proposed to characterize non-local feature of entanglement(48); (49). The Renyi entropies of a subsystem with reduced density matrix are defined as ; limit gives the Von Neuman entropy. For a topologically ordered state, Renyi entropies have the form , where is the boundary of the subsystem, and all other terms vanish in the large limit; is a non-universal constant, while a positive is a correction to the area law of entanglement and reaches a universal value determined by total quantum dimension of quasiparticle excitations(48); (49). Previous DMRG study(40) found in the intermediate region of consistent with a SL for this model. We compute the entanglement entropy (EE) on long cylinders by partitioning the system in the middle along the vertical direction. For each fixed , we fit the entropy to limit to find the entropy of a possible minimum entropy state(16).

In Fig. 17, we show our DMRG results for the EE at and on both TC and RC cylinders. We obtain accurate EE when . For , we extrapolate the EE with the DMRG truncation error, which has significant uncertainty from the extrapolation. On RC cylinder, we perform linear fit of the EE versus using the three largest sizes. We find the TEE at is close to , while at is close to . However, the system appears to have large finite-size effects, which can be seen by comparing the results on the RC and TC cylinders. On the TC cylinder, the linear fits of the EE vs give the TEE close to zero, which is different from the RC cylinder. Similar effect has also been observed in the - model on the honeycomb lattice(51); (52). Because of such strong finite-size effects, the TEE obtained by fitting EE on our small sizes may not be able to distinguish different quantum phases in the - square lattice model.

### References

- L. Balents, Nature (London) 464, 199 (2010).
- P. A. Lee, N. Nagaosa, and X. G. Wen, Rev. Mod. Phys. 78, 17 (2006).
- R. Moessner and S. L. Sondhi, Phys. Rev. Lett. 86, 1881 (2001).
- C. Nayak and K. Shtengel, Phys. Rev. B 64, 064422 (2001).
- T. Senthil and O. Motrunich, Phys. Rev. B 66, 205104 (2002).
- L. Balents, M. P. A. Fisher, and S. M. Girvin, Phys. Rev. B 65, 224412 (2002).
- D. N. Sheng and L. Balents, Phys. Rev. Lett. 94, 146805 (2005).
- S. V. Isakov, Y. B. Kim, and A. Paramekanti, Phys. Rev. Lett. 97, 207204 (2006).
- S. V. Isakov, M. B. Hastings, and R. G. Melko, Nat. Phys. 7, 772 (2011).
- Y. C. He, D. N. Sheng, and Y. Chen, Phys. Rev. B 89, 075110 (2014).
- A. Kitaev, Ann. Phys. (Amsterdam) 321, 2 (2006).
- H. Yao and S. A. Kivelson, Phys. Rev. Lett. 108, 247206 (2012).
- H. C. Jiang, Z. Y. Weng, and D. N. Sheng, Phys. Rev. Lett. 101, 117203 (2008).
- S. Yan, D. Huse, and S. R. White, Science 332, 1173 (2011).
- S. Depenbrock, I. P. McCulloch, and U. Schollwöck, Phys. Rev. Lett. 109, 067201 (2012).
- H. C. Jiang, Z. H. Wang, and L. Balents, Nat. Phys. 8, 902 (2012).
- N. Read and S. Sachdev, Phys. Rev. Lett. 66, 1773 (1991).
- X. G. Wen, Phys. Rev. B 44, 2664 (1991).
- X. G. Wen, Phys. Rev. B 40, 7387 (1989).
- L. Balents, M. P. A. Fisher, and C. Nayak, Phys. Rev. B 60, 1654 (1999).
- T. Senthil and M.âP.âA. Fisher, Phys. Rev. B 62, 7850 (2000); Phys. Rev. Lett. 86, 292 (2001).
- P. Chandra and B. Doucot, Phys. Rev. B 38, 9335(R) (1988).
- L. B. Ioffe and A. I. Larkin, Int. J. Mod. Phys. B 02, 203 (1988).
- S. Sachdev and R. N. Bhatt, Phys. Rev. B 41, 9323 (1990).
- Andrey V. Chubukov and Th. Jolicoeur, Phys. Rev. B 44, 12050 (1991).
- M. E. Zhitomirsky and K. Ueda, Phys. Rev. B 54, 9007 (1996).
- A. E. Trumper, L. O. Manuel, C. J. Gazza, and H. A. Ceccatto, Phys. Rev. Lett. 78, 2216 (1997).
- R. R. P. Singh, Z. Weihong, C. J. Hamer, and J. Oitmaa, Phys. Rev. B 60, 7278 (1999).
- L. Capriotti and S. Sorella, Phys. Rev. Lett. 84, 3173 (2000).
- L. Capriotti, F. Becca, A. Parola, and S. Sorella, Phys. Rev. Lett. 87, 097201 (2001).
- G. M. Zhang, H. Hu, and L. Yu, Phys. Rev. Lett. 91, 067201 (2003).
- K. Takano, Y. Kito, Y. Ono, and K. Sano, Phys. Rev. Lett. 91, 197202 (2003).
- M. Mambrini, A. Läuchli, D. Poilblanc, and F. Mila, Phys. Rev. B 74, 144422 (2006).
- R. Darradi, O. Derzhko, R. Zinke, J. Schulenburg, S. E. Krüger, and J. Richter, Phys. Rev. B 78, 214415 (2008).
- L. Isaev, G. Ortiz, and J. Dukelsky, Phy. Rev. B 79, 024409 (2009).
- J. Richter and J. Schulenburg, Eur. Phys. J. B 73, 117 (2010).
- J. Reuther and P. Wölfle, Phys. Rev. B 81, 144410 (2010).
- J. F. Yu, Y. J. Kao, Phys. Rev. B 85, 094407 (2012).
- L. Wang, Z. C. Gu, F. Verstraete, and X. G. Wen, arXiv:1112.3331.
- H. C. Jiang, H. Yao, and L. Balents, Phys. Rev. B 86, 024424 (2012).
- K. S. D. Beach, Phys. Rev. B 79, 224431 (2009).
- F. Mezzacapo, Phys. Rev. B 86, 045115 (2012).
- T. Li, F. Becca, W. J. Hu, and S. Sorella, Phys. Rev. B 86, 075111 (2012).
- L. Wang, D. Poilblanc, Z. C. Gu, X. G. Wen, and F. Verstraete, Phys. Rev. Lett. 111, 037202 (2013).
- W. J. Hu, F. Becca, A. Parola, and S. Sorella, Phys. Rev. B 88, 060402 (2013).
- R. L. Doretto, Phys. Rev. B 89, 104415 (2014).
- Y. Qi and Z. C. Gu, Phys. Rev. B 89, 235122 (2014).
- A. Kitaev, J. Preskill, Phys. Rev. Lett. 96, 110404 (2006).
- M. Levin, X.-G. Wen, Phys. Rev. Lett. 96, 110405 (2006).
- R. Ganesh, Jeroen van den Brink, and S. Nishimoto, Phys. Rev. Lett. 110, 127203 (2013).
- Zhenyue Zhu, D. A. Huse, and S. R. White, Phys. Rev. Lett. 110, 127205 (2013).
- S. S. Gong, D. N. Sheng, Olexei I. Motrunich, and Matthew P. A. Fisher, Phys. Rev. B 88, 165138 (2013).
- T. Senthil, L. Balents, S. Sachdev, A. Vishwanath, and Matthew P. A. Fisher, Phys. Rev. B 70, 144407 (2004).
- T. Senthil, A. Vishwanath, L. Balents, S. Sachdev, and Matthew P. A. Fisher, Science 303, 1490 (2004).
- I. P. McCulloch and M. Gulácsi, Europhys. Lett. 57, 852 (2002); I. P. McCulloch, J. Stat. Mech. (2007) P10014.
- A. W. Sandvik, Phys. Rev. Lett. 98, 227202 (2007).
- R. G. Melko and R. K. Kaul, Phys. Rev. Lett. 100, 017203 (2008).
- F. Jiang, M. Nyfeler, S. Chandrasekharan, and U. Wiese, J. Stat. Mech. (2008) 02009.
- J. Lou, A. W. Sandvik, and N. Kawashima, Phys. Rev. B 80, 180414(R) (2009).
- J. Lou and A. W. Sandvik, Phys. Rev. B 80, 212406 (2009).
- A. W. Sandvik, Phys. Rev. Lett. 104, 177201 (2010).
- Matthew S. Block, R. G. Melko, and R. K. Kaul, Phys. Rev. Lett. 111, 137202 (2013).
- S. Pujari, K. Damle, and F. Alet, Phys. Rev. Lett. 111, 087203 (2013).
- A. W. Sandvik, Phys. Rev. B 85, 134407 (2012).
- We get much better convergence on cylinder than on torus. We can achieve the truncation error for cylinder and for cylinder. For example, for and , we get the truncation error on cylinder by keeping states, while the error is much larger on torus even when we keep states.
- S. R. White and A. L. Chernyshev, Phys. Rev. Lett. 99, 127004 (2007).
- We get converged for by keeping more than states, while the results for are obtained through extrapolations with the DMRG truncation error.
- See Supplementary Material for more details.
- In DMRG calculations for , one needs to extrapolate with the DMRG truncation error, which may give some uncertainty to . 1