## Abstract

We use the density matrix renormalization group to perform accurate calculations of the ground state of the nearest-neighbor quantum spin *S* = 1/2 Heisenberg antiferromagnet on the kagome lattice. We study this model on numerous long cylinders with circumferences up to 12 lattice spacings. Through a combination of very-low-energy and small finite-size effects, our results provide strong evidence that, for the infinite two-dimensional system, the ground state of this model is a fully gapped spin liquid.

We consider the quantum spin *S* = 1/2 kagome Heisenberg antiferromagnet (KHA) with only nearest-neighbor isotropic exchange interactions (Hamiltonian *1*–*3*). A spin liquid is a magnetic system that has “melted” in its ground state because of quantum fluctuations, so it has no spontaneously broken symmetries (*4*). A key problem in searching for spin liquids in two-dimensional (2D) models is that there are no exact or nearly exact analytical or computational methods to solve infinite 2D quantum lattice systems. For 1D systems, the density matrix renormalization group (DMRG) (*5*, *6*), the method we use here, serves in this capacity. In addition to its interest as an important topic in quantum magnetism, the search for spin liquids thus serves as a test-bed for the development of accurate and widely applicable computational methods for 2D many-body quantum systems.

The KHA has been studied with approximate approaches for decades (*1*), with proposals for spin-liquid and valence-bond–crystal (VBC) ground states (*2*, *3*). In the past few years, numerical evidence had suggested that the KHA ground state might be a VBC with a 36-site unit cell of bonds with spatially nonuniform spin-spin correlations in the ground state that break the translational symmetry of the kagome lattice. First proposed in (*2*), and explored in more detail in (*7*), this honeycomb VBC (HVBC) was also studied with a perturbative series expansion (*8*, *9*). The series for the ground-state energy rapidly converges to a low energy (*8*, *9*). Subsequently, Evenbly and Vidal (*10*) used the multiscale entanglement renormalization ansatz (MERA), obtaining the HVBC and an energy close to the series result. Another approximate numerical approach based on an effective quantum dimer model also yielded the HVBC as the ground state, but noted that a *Z*_{2} spin liquid is close by in a generalized parameter space (*11*). A recent DMRG study (*12*), in contrast, found a spin-liquid ground state, but on the largest lattices the energy obtained was substantially above that of the HVBC, suggesting that the method had not found the true ground state. The simple KHA is clearly near a number of different possible ground-state phases upon adding small perturbations to its Hamiltonian [such as a second-neighbor exchange interaction (*13*)], and it is this proximity of many competing states that has made this such a challenging system. A variational state that may serve as a multicritical state between all these possible ground-state phases is the critical gapless spin liquid of (*14*).

We performed an extensive DMRG study of the KHA with important differences in technique from the previous study. Most importantly, we studied long cylinders with open ends, thus avoiding fully periodic (toroidal) boundary conditions, which are known to greatly magnify the truncation errors in the DMRG (*6*). We also show why the series expansions appear to converge well, yet reach the incorrect HVBC state: The path connecting the series expansion starting point and the KHA is interrupted by a first-order phase transition (*15*). For details on the competition between the spin liquid and HVBC states, see figs. S1 to S3.

DMRG efficiently finds the ground state for long 1D systems, but for 2D systems one must study cylinders or strips of limited width. We studied mostly cylinders with a maximum circumference of 12 lattice spacings, obtaining ground-state energies with uncertainties of less than 0.1%, with much higher accuracy for the narrower cylinders. We labeled the cylinders by their orientation, circumference, and any shift in wrapping the cylinder periodically (Fig. 1A) (*15*). For example, “YC9-2” denotes a cylinder (C) with some of the bonds oriented in the *y* direction (Y), with circumference of nine lattice spacings and a shift of two columns when connected periodically. We estimate the energy per site for each type of infinitely long cylinder by subtracting energies of cylinders of different length to remove end effects.

Our DMRG results for ground-state energies, along with results from other approaches, are shown in Fig. 1B. The DMRG energies are consistent with the Lanczos results (*16*, *17*) and well below the energies of MERA (*10*) and the series expansion for the HVBC (*8*, *9*). The previous DMRG result (*12*) is close to the true ground state (*17*) for a torus with a circumference of six (circumference-6) lattice spacings. The entanglement across a cut that separates a circumference-six torus into two parts is roughly the same as that across a cut that separates a circumference-12 cylinder—our current limit. Thus, it is not surprising that the DMRG energies for tori (*12*) are overestimates for circumferences larger than six lattice spacings. We also obtain a new rigorous upper bound on the ground-state energy of the infinite 2D system using the DMRG on strips with open boundary conditions (*15*); this is well below the previous bound from MERA (*10*).

The 36-site unit cell of the HVBC does fit on the XC8 cylinder (1/*c* ≅ 0.14, where *c* is the circumference), allowing a direct comparison between the HVBC series and our DMRG: the DMRG energy is lower by 0.004(1). Comparing this cylinder to the corresponding torus, the energy shows strong finite-size effects in the HVBC series (*8*, *9*), but the same comparison between Lanczos (*16*) and our DMRG shows that the true finite-size effect is much smaller (*15*). These finite-size effects remain small for even smaller circumferences (*17*). This is consistent with the small correlation length that we find: less than 1.5 lattice spacings (*15*). We conclude that the ground-state energies of our widest cylinders have minimal finite-size effects and thus provide reliable estimates of the 2D energy.

The ground state that we find has only short-range correlations, with a nonzero energy gap for any excitations, including spin-singlet excitations. We have tested the response of this ground state to many sorts of perturbations that would select out ordered states, if they exist, without detecting any signs of any ordering. Thus, we conclude that the ground state is a gapped spin liquid. Such a state can be represented as some sort of short-range resonating–valence-bond (RVB) state (*18*–*24*). The shortest resonant loops of singlet dimers in a nearest-neighbor RVB state on a kagome lattice each surround only one hexagon of the lattice (*7*). If the RVB state had all dimer covers equally weighted, all 32 of these elementary resonant loops would be equally present. Instead, the ground state appears to substantially overemphasize certain eight-site loops of a diamond shape (Fig. 2B). To test the response of the ground state to enhancing each of these 32 elementary resonant loops, we slightly increased the exchange couplings along the bonds of such a loop at the center of a YC8 cylinder and measured how much this enhanced the spin-spin correlations along the loop and elsewhere. It is the eight-site diamond loop that elicits the strongest response (Fig. 2B). The six-site “perfect hexagon” loop (Fig. 2A) (*2*, *7*, *8*), shows a much smaller response, suggesting that this resonant loop is actually underweighted in the ground state.

One cannot tile a kagome lattice with just these favored resonant diamonds. However, a particular “diamond-pattern” VBC (Fig. 2C), appears to be closely related to the spin liquid, and it is useful to think of the spin liquid as a melted state of this crystal. We have measured the correlation length for VBC energy correlations in this diamond pattern along cylinder YC8 and find that it is less than 1.5 lattice spacings (*15*). Unlike the HVBC, which shows a first-order phase transition to the spin liquid (*15*), this diamond VBC evolves smoothly into the spin liquid without any phase transition as one changes the strengths of the exchanges that are altered to favor it (we call this “pinning it,” Fig. 2C). For the even cylinders on which this diamond VBC does fit, this process allows a careful production of the spin liquid by approaching it in a smooth fashion from the diamond VBC (*15*).

The infinitely long cylinders may be viewed as one-dimensional systems with a unit cell containing *N*_{c} spins. Even-*N*_{c} cylinders (e.g., *N*_{c} = 12 for YC8) are compatible with the diamond VBC, and the ground state of the infinite cylinder appears to be nondegenerate. Odd-*N*_{c} cylinders are not compatible with the diamond VBC, and the Lieb-Schultz-Mattis theorem implies that the ground state must be degenerate (*25*). Ground states on these odd cylinders weakly break translational invariance, spontaneously doubling the unit cell, which produces a pair of degenerate ground states, still with a gap to higher excited singlet states. The symmetry breaking is in a “striped” pattern (Fig. 3). For YC6 and YC10, the stripes run around the circumference, whereas for the other odd spiral cylinders the stripes are spirals.

The ends of our cylinders may have low-lying states below the bulk singlet and triplet gaps. The following DMRG procedure avoids these edge states: First target only one state, and sweep enough to obtain a high-accuracy ground state. Then restrict the range of bonds that are updated in the DMRG sweeps to the central half of the sample and target the two lowest-energy states, again sweeping to high accuracy, but keeping the end regions of the samples locally in the ground state. This technique is particularly important for obtaining the singlet gap. For the triplet gap, one can also keep the excitation away from the ends with local magnetic fields. With this approach, we can target both states together—one with total *S*_{z} = 0 and the other with *S*_{z} = 1—or run them separately. These different approaches allowed for fairly independent checks on the results. Figure 4 shows the measured bulk gaps. Gaps are more demanding than ground-state energies, so we do not have gap estimates for our widest cylinder (*15*). See table S1 for energies and gaps of all the cylinders studied.

The singlet gap is 0.050, within the errors, for the even XC8 and YC8 cylinders, and it remains near this value for the wider even XC12-2 cylinder. The odd cylinders come in at least two families: YC6 and YC10 are not spirals, whereas YC5-2, YC7-2, and YC9-2 are a series of spirals with increasing circumference. In each of these (small) families of odd cylinders, the singlet gap increases as the circumference increases, supporting our conclusion that the singlet gap is ~0.05 in the 2D limit. This is quite different from the exact diagonalization results, where there are many lower-lying singlets, and the lowest singlet gap is only ~0.01 on the standard 36-site torus (*26*). The singlet gap is a strong function of a second-neighbor exchange coupling (*J*_{2}) (*13*), with an apparent phase transition at a very small ferromagnetic value of *J*_{2}. The location of this nearby transition is quite sensitive to circumferences, and this produces the large differences in the singlet gap between the torus and our cylinders.

The triplet (spin) gap on the 36-site torus is 0.164 from exact diagonalization (*26*). Although XC8 and YC8 have gaps that are quite close to this, our results on other cylinders suggest that the 2D triplet gap is smaller (Fig. 4). The triplet excitations are composed of two spinons, but we cannot resolve whether or not the two spinons bind, although in some cylinders any binding must be very weak. This composite nature of the excitation seems to make the finite-size effects and variation between the cylinders more pronounced. We do not yet understand the details of these effects. As for the singlet gap, the spiral odd cylinders have the smallest triplet gaps, and the even cylinders the largest. The triplet gap remains above the singlet gap in all the systems we have studied; thus, we believe it remains nonzero in the 2D limit.

A nearest-neighbor RVB wave function is a linear combination of nearest-neighbor singlet dimers covering the kagome lattice (*24*). For a kagome lattice wrapped on a cylinder, such dimer covers are in two topologically distinct sectors that differ by a *Z*_{2} winding number, and dimer resonances on finite loops do not change this winding number (*19*, *20*, *22*–*24*). We can force our even cylinders to have one or the other of the two winding numbers by choosing how many spins to leave at each end. For a finite circumference *c*, these two sectors have different ground-state energies (*E*) per site; for YC8, their difference is δ*E* ≈ 0.00069(3). For odd cylinders, the two sectors are related by translation along the length of the cylinder and, thus, are degenerate. If an even cylinder is topologically ordered, one expects δ*E* ∼ exp(−*c*/ξ) (where ξ is a correlation length ) at large *c*, but we do not yet have reliable δ*E* results for larger cylinders to test this hypothesis. This is partly because the singlet gap above the ground state in the higher-energy sector is substantially smaller than the singlet gap above the overall ground state, slowing DMRG’s convergence.

A domain wall along the cylinder where the *Z*_{2} winding number changes is a spinon. For the odd cylinders, the degeneracy of the two sectors means the spinons are unconfined. However, two spinons might bind with a finite binding energy—it appears that this might be happening on cylinder YC10 and, thus, it might also occur in the 2D limit. For the even cylinders, on the other hand, the spinons are confined by an effective potential that grows linearly with the distance, because the domain between them is in the higher-energy sector. As a result, the excitation across the spin gap for the even cylinders is a bound spinon pair. It remains to be determined if they remain bound in the 2D limit.

Much remains to be understood concerning the low-energy behavior of the KHA, particularly the detailed structure, exchange statistics, and dispersion relations of the various excitations. It will be instructive to also explore the phase diagram in the vicinity of this simple nearest-neighbor–only Heisenberg model by changing the Hamiltonian in various ways (*13*) to find what other phases are nearby and perhaps to move “deeper” into this spin liquid where it might be easier to study.

## Supporting Online Material

www.sciencemag.org/cgi/content/full/science.1201080/DC1

Materials and Methods

SOM Text

Figs. S1 to S3

Table S1

References (*5*, *6*, *8*–*10*, *12*, *16*, *17*, *27*–*29*)

## References and Notes

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
Supporting material is available on
*Science*Online. - ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
**Acknowledgments:**We thank A. Lauchli for Lanczos results; R. Singh for series results; and C. Lhuillier, M. P. A. Fisher, T. Senthil, E. Sorensen, S. Sondhi, G. Vidal, X.-G. Wen, O. Tchernyshyov, and M. Stoudenmire for discussions. This work was supported by NSF grants DMR-0907500 and DMR-0819860.