## Abstract

Quantum spin liquids (QSLs) are exotic phases of matter that host fractionalized excitations. It is difficult for local probes to characterize QSL, whereas quantum entanglement can serve as a powerful diagnostic tool due to its nonlocality. The kagome antiferromagnetic Heisenberg model is one of the most studied and experimentally relevant models for QSL, but its solution remains under debate. Here, we perform a numerical Aharonov-Bohm experiment on this model and uncover universal features of the entanglement entropy. By means of the density matrix renormalization group, we reveal the entanglement signatures of emergent Dirac spinons, which are the fractionalized excitations of the QSL. This scheme provides qualitative insights into the nature of kagome QSL and can be used to study other quantum states of matter. As a concrete example, we also benchmark our methods on an interacting quantum critical point between a Dirac semimetal and a charge-ordered phase.

## INTRODUCTION

Quantum spin liquids (QSLs) are highly entangled states of matter with exotic excitations behaving as fractions of fundamental particles (*1*, *2*). In the vigorous search for candidate materials, herbertsmithite ranks as one of the most promising ones (*3*). Although it displays several signatures of a QSL, consensus between experimental and theoretical studies is hindered not only by disorder but also by the lack of understanding for the minimal model (*3*). A starting point for a theoretical description of this correlated material is the antiferromagnetic Heisenberg model on the kagome lattice, which is built out of corner-sharing triangles that frustrate the anti-alignment of spins favored on each bond. Frustration renders this model difficult to solve (*4*–*17*), and debates between different theoretical scenarios persist. For instance, in a numerical study using the density matrix renormalization group (DMRG) (*18*), a gapped ground state without magnetic order was found (*11*). More recent numerical studies (*15*–*17*) suggest a gapless QSL. In this direction, evidence for a Dirac QSL (*5*, *7*) was obtained using DMRG simulations (*17*).

The fractional excitations of a QSL cannot be characterized by local order parameters. Instead, entanglement may directly reveal fractionalization by virtue of its nonlocal nature. For example, the so-called topological entanglement entropy (EE) can be used to detect gapped QSL [e.g., see the review (*19*)]. Less comprehensive but nevertheless interesting results for the entanglement properties of gapless systems were obtained (*20*–*26*). For example, attempts have been made to understand the entanglement response to a flux insertion (*21*, *27*, *28*). However, our understanding of interacting systems remains limited, especially for realistic QSL models.

Here, we investigate the quantum entanglement of the QSL in the kagome antiferromagnetic Heisenberg model in response to a magnetic flux. We perform large-scale DMRG simulations (see Methods) on infinitely long cylinders through which the flux is threaded (see Fig. 1A). The main finding is that the EE is highly sensitive to the flux and is consistent with emergent Dirac cones for fractionalized spinon excitations. This constitutes new evidence that the gapless Dirac QSL is the ground state of the kagome antiferromagnetic Heisenberg model. Moreover, to illustrate these entanglement signatures in a simpler setting, we will begin by studying a strongly interacting quantum phase transition between a Dirac semimetal and a charge-ordered state. These results not only help with the interpretation of the data for the debated kagome QSL but also shed new light on quantum critical states of matter.

## RESULTS

### Entanglement scaling of Dirac fermions

We consider a general quantum system on an infinitely long cylinder, and we calculate the von Neumann EE *S* of the ground state by partitioning the system into two halves, as shown in Fig. 1A. *S* quantifies the amount of quantum entanglement between the halves and takes the form (*19*): , where *a* ≪ *L*_{y} is a microscopic scale such as a lattice spacing. The first term arises for the ground states of most Hamiltonians and is called the “boundary law” because it scales with the length of the partition, the circumference *L*_{y} of the cylinder. The boundary law term is of little interest in itself because it is not universal, as it depends on the microscopic scale *a*. In contrast, the subleading term γ is a low-energy property and does not depend on *a*. Additional information can be extracted by inserting a flux Φ in the cylinder and studying the response of γ to Φ. The flux dependence can then be used as a fingerprint to identify the quantum state (*21*, *27*).

For a two-component free (gapless) Dirac fermion on the cylinder, the EE takes the following form in the continuum (see the supplementary materials)(1)where *B* = 1/6 (*27*, *28*). One way to understand the above scaling behavior is to realize that the transverse momenta are quantized on a cylinder and the flux Φ will move those quantized momenta toward or away from the Dirac point. Intuitively, the subleading term quantifies how far the quantized momenta are from the Dirac point. When Φ → 0, 2π, one momentum exactly hits the Dirac point, leading to a diverging γ. When Φ = π, the momenta are farthest away from the Dirac point so *S* becomes minimal.

To study interacting systems on a lattice, based on Eq. 1, we propose the following ansatz for the EE scaling function(2)

There are three new ingredients. First, we can have *N* > 1 different Dirac fermions. Second, the momenta of the corresponding Dirac points in the Brillouin zone can be different, which is encoded in the shift . is proportional to the flux at which the *n*th Dirac fermion’s gap vanishes on the cylinder; 1/*s* is the proportionality constant. Third, the Dirac fermions can carry a fractional charge *s* (e.g., *s* = ±1/2 in the kagome QSL); hence, they will feel a flux *s*Φ instead of Φ. The flux response of *S* thus gives a clear way to identify fractionalization, which is notoriously difficult using conventional approaches.

### Quantum criticality

We start our discussion by explaining the salient entanglement features of a quantum critical transition between a Dirac semimetal and an interaction-driven insulator with charge order. By virtue of its universality, such a transition is relevant in contexts such as charge density wave transitions in graphene (*29*). We consider the π-flux square lattice model with a short-ranged repulsion *V*(3)where is the creation operator for a spinless fermion on site *i* and *n*_{i} is the particle number operator. The phase factor generates a π flux on each square plaquette as shown in Fig. 1B. We consider the filling factor one-half, where Fermi level is set to be zero. In the noninteracting limit, *V* = 0, the band structure hosts two Dirac cones located at (*k*_{x}, *k*_{y}) = (±π/2, π/2), as shown in Fig. 1C. The repulsive interaction between nearest neighbors drives a quantum phase transition from the Dirac semimetal to a charge density wave phase through the strongly interacting Gross-Neveu-Yukawa quantum critical point (*29*), where the Dirac quasiparticles are destroyed by quantum fluctuations. Numerical studies studying local observables were performed (*30*, *31*), but the entanglement properties near the critical point have not been investigated.

Figure 1 shows the EE *S* of model (Eq. 3) on infinite cylinders threaded by a flux Φ. First, we find that *S*(Φ) agrees with the scaling function Eq. 2. Specifically, we have two (nonfractionalized) Dirac fermions: *N* = 2 and *s* = 1. For one type of cylinder (YC8-0; see the caption of Fig. 1), the momenta are quantized such that they hit the Dirac points when Φ = 0. Therefore, in the scaling function Eq. 2. For the other type of cylinder (YC8-2), the quantized momenta hit the Dirac points when Φ = π; hence, we have .

The scaling behavior is robust in the entire Dirac semimetal phase *V* < *V*_{c} ≈ 1.3*t* (Fig. 1, F to I), despite the decrease of the prefactor *B* as the quantum critical point is approached. In the charge-ordered phase *V* > *V*_{c}, the entropy does not follow the scaling behavior Eq. 2 anymore (see the supplementary materials). Last, we emphasize that the scaling behavior is robust against changes to the circumference and cylinder type. As shown in Fig. 1J, for various cylinder sizes and types, the scaling parameter *B* follows the same decreasing trend as the critical point is approached. As a consistency check, we observe that *B* approaches its noninteracting value *B* = 1/6 when *V* → 0. To confirm that the above properties of the EE are universal, we analyzed another model on the honeycomb lattice and reached identical conclusions (see the supplementary materials).

The deviation of *B* from the expectation of free theory is a consequence of the increasing quantum fluctuations as the critical point is approached. To explain this fact, we now invoke a field theory description. To make this theory tractable, we extend the number of Dirac fermions to *N* ≫ 1. In this limit, it was recently shown (*28*) that the subleading correction γ vanishes at leading order. This represents a marked reduction compared with a weakly interacting Dirac semimetal, where γ is directly proportional to the number of Dirac fermions, *N*. Extrapolating to finite *N*, we conjecture that as the quantum critical point of Eq. 3 is approached, *B* is suppressed. Our data in Fig. 1 corroborate this conclusion.

### Kagome spin liquid

We now tackle our main objective, the spin-1/2 antiferromagnetic Heisenberg model on the kagome lattice (Fig. 2)(4)where *J*_{1}, *J*_{2} > 0 are nearest and next-to-nearest neighbor antiferromagnetic couplings, respectively. Although we focus on the *J*_{2} = 0 case, we shall also consider the effects of a small *J*_{2}/*J*_{1}, which makes the numerical results more stable. Figure 3 shows the flux dependence of the EE at *J*_{2}/*J*_{1} = 0, 0.05, 0.1, as well as for two cylinder types, YC8-0 and YC8-2. We first note that in all cases, *S* strongly depends on Φ, which is a hallmark of low-energy excitations. In contrast, a state with a large gap would be essentially insensitive to Φ. The data in Fig. 3 can be accurately fitted with the scaling function Eq. 2. The parameters (*N*, *s*, and ) are chosen to match the π-flux Dirac QSL (*5*, *7*), in which four two-component Dirac spinons (*N* = 4 accounts for the spin and valley degrees of freedom) carry fractionalized spin *s* = 1/2. The shifts depend on the cylinder type (Fig. 2) and are given in Table 1.

In Fig. 3, we observe that the scaling function Eq. 2 accurately fits the data for all the geometries and couplings considered. When *J*_{2} = 0, the fitting parameter *B* takes the value 0.219 on YC8-0 and 0.240 on YC8-2. This is larger than the free fermion value *B* = 1/6. This deviation from the noninteracting value can be understood by the fact that the low-energy description of the Dirac QSL is in terms of Dirac spinons strongly coupled to an emergent photon (*7*). In contrast, in the large *N* limit, when the gauge fluctuations become suppressed, the leading order EE is that of *N* free Dirac fermions (see the supplementary materials). Our data suggest that the value of *B* becomes renormalized at finite *N* but the flux dependence remains largely unchanged.

At *J*_{2}/*J*_{1} = 0.05 and 0.1, we see that the fits are more accurate. *B* increases slightly to 0.250 for YC8-0. On YC8-2, *B* first decreases to 0.166 and then to 0.125. This tendency for *B* to be smaller on type 2 cylinders was seen above for the interacting Dirac semimetal on the square lattice. Here, the finite size effects for YC8-2 may be further amplified compared to YC8-0. One reason may be that the allowed momentum lines are always closer to the gapless Dirac points on the YC8-2 cylinder. As *J*_{2}/*J*_{1} is increased, a conventional ordered state becomes favored and the resulting phase transition is expected to leave an imprint on the entanglement. A more detailed analysis in this direction will likely lead to new insights into the properties of the kagome QSL and its phase transitions.

## DISCUSSION

In summary, by monitoring the EE response to a flux threaded in a cylindrical geometry, we were able to gain new insights about two physical systems: (i) a quantum critical phase transition of itinerant electrons and (ii) the frustrated kagome Heisenberg model. In the first case, the EE tracks the evolution of the Dirac fermions as the quantum critical point is approached. For the kagome model, the flux dependence of the EE unambiguously points to four emergent Dirac cones of fractionalized excitations (spinons). The robust features we have identified on various cylinder types and values of the Heisenberg couplings strongly suggest that the kagome Heisenberg model is a gapless Dirac QSL. These new insights will help with the modeling of candidate materials such as herbertsmithite. Our two concrete examples give us confidence that entanglement signatures will become a valuable tool in the investigation of a broad class of quantum states of matter.

In closing, we would like to make remarks about numerical finite size effects. While we found strong evidence for gapless features in the kagome Heisenberg model for the system sizes that we can access, caution is needed when making a statement about the thermodynamic limit. These finite-size calculations are not optimal to distinguish a gapless ground state from a ground state with a large spatial correlation length (compared to the circumference of the cylinder). Thus, we cannot rule out the possibility of a gapped ground state with a small but finite gap for the kagome Heisenberg model in the thermodynamic limit. To resolve this issue, numerical simulations beyond what is currently feasible are required. However, if a small gap is found, the gapless Dirac spin liquid description will likely be the most practical in many situations, such as when temperature exceeds the gap scale.

## METHODS

The ground states of Eqs. 3 and 4 were determined using the DMRG (*18*), which is a powerful algorithm to determine, in an unbiased fashion, the low-lying states of quantum systems. In our simulations, we worked on infinitely long cylinders with a finite circumference (*32*). We could reach circumferences of four unit cells on the kagome lattice, which is close to the current computational limit. In our simulations, matrix product states of bond dimension 1600 were sufficient to describe the EE of the π-flux model on the square lattice, whereas for the kagome model, a bond dimension of 6000 was used, as the subsystem entanglement is significantly larger. The numerical flux insertion experiment was performed by adiabatically changing (twisting) boundary conditions in the Hamiltonian. In the simulations, we imposed twisted boundary conditions along the circumference of the cylinder by replacing the terms () for all bonds crossing the *y* boundary with () (*33*).

Once the ground state |Ψ(Φ)〉 was computed, we partitioned the cylinder into two halves, *A* and *B*, and calculated the von Neumann EE *S*(Φ) = − ∑_{i} λ_{i}(Φ)ln λ_{i}(Φ), where λ_{i} are the eigenvalues of reduced density matrix of the *A* half, ρ_{A}(Φ>) = *Tr*_{B}|Ψ(Φ)〉〈Ψ(Φ)|. The EE measures the amount of quantum entanglement between a region *A* and its complement. To obtain the EE at different Φ, we used an adiabatic scheme: The ground state |Ψ(Φ)〉 is taken as the initial state for the calculation at Φ + ΔΦ.

The fit of *S*(Φ) was based on the least-squares method. The data points near the entropy minimum were used in the fitting as these are the most reliable. For the kagome lattice, the entropy in the range |Φ − Φ^{min}| < 0.24π was used for the fitting process. For the square lattice, the entropy in the range |Φ − Φ^{min}| < 0.4π was used in the fits, where Φ^{min} is the flux value where the entropy is minimal. We had verified that all the fits are stable and independent of the data range we select, except for the YC8-2 *J*_{2} = 0 case (Fig. 3D). For YC8-2 *J*_{2} = 0, the scaling parameter *B* could vary from 0.20 to 0.25, as we changed the data regime from |Φ − Φ^{min}| < 0.2π to |Φ − Φ^{min}| < 0.4π. This could be attributed to the stronger finite-size effects at this coupling on the YC8-2 cylinder.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/4/11/eaat5535/DC1

Section S1. EE for free Dirac fermions on the cylinder

Section S2. EE of interacting Dirac fermions at large *N*

Section S3. Quantum critical point of Dirac fermions on the honeycomb lattice

Section S4. EE in the gapped phase

Fig. S1. Bipartition of an infinite cylinder, which is threaded by a flux Φ.

Fig. S2. Entanglement of quantum critical Dirac fermions on the honeycomb lattice.

Fig. S3. Entanglement in the charge-ordered phase.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We are grateful for discussions with H. J. Changlani, E. Fradkin, J. Maciejko, S. Sachdev, C. Wang, and S. Whitsitt. Y.-C.H. thanks M. Zaletel, M. Oshikawa, and F. Pollmann for previous collaboration on a related project.

**Funding:**W.Z. was supported by the U.S. Department of Energy through the Los Alamos National Laboratory LDRD Program. X.C. was supported by a postdoctoral fellowship from the Gordon and Betty Moore Foundation, under the EPiQS initiative, grant GBMF4304, at the Kavli Institute for Theoretical Physics. Y.-C.H. is supported by the Gordon and Betty Moore Foundation under the EPiQS initiative, GBMF4306, at Harvard University. Y.-C.H. was supported, in part, by Perimeter Institute for Theoretical Physics. Research at Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development Canada and by the Province of Ontario through the Ministry of Research, Innovation and Science. W.W.-K. was funded by a Discovery Grant from NSERC and by a Canada Research Chair. The work was initiated at a Moore funding postdoc symposium in Aspen. Part of the work was performed at the Aspen Center for Physics, which is supported by NSF grant PHY-1066293.

**Author contributions:**X.C. and Y.-C.H. initiated the project. W.Z. and Y.-C.H. performed the DMRG simulations. All authors contributed equally to the analysis of the data and writing of the manuscript.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2018 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).