Entanglement entropy as a portal to the physics of quantum spin liquids

Quantum spin liquids (QSLs) are phases of strongly interacting spins that do not order even at absolute zero temperature and, therefore, are not characterized by a Landau order parameter [1]. They are associated with remarkable phenomena such as fractional quantum numbers [2–5], transmutation of statistics (e.g. fermions appearing in a purely bosonic model) [6, 7] and the enabling of otherwise impossible quantum phase transitions [8], to name just a few. Spin liquids may be gapless or gapped. While current experimental candidates for spin liquids appear to have gapless excitations [9], gapped spin liquids are indicated in numerical studies on the Kagome [10, 11], square [12, 13] and honeycomb lattice [14]. A distinctive property of gapped spin liquids is that they are characterized by topological order—i.e. ground-state degeneracy that depends on the topology of the underlying space [15, 16]. The gapless spin liquids are relatively difficult to characterize and their low-energy description often involves strongly interacting matter-gauge theories [17].

In the theoretical search for QSLs, two qualitatively different kinds of questions can be asked:

• Given a particular Hamiltonian H, does it admit a QSL as its ground-state? If yes, can one further characterize the QSL?
• Given a ground-state wavefunction(s), can one tell if it corresponds to a QSL? If yes, can one further characterize the QSL?

These two questions are related: if a numerical scheme can provide one with access to the approximate ground-state wavefunction of a particular Hamiltonian H, then answering the second question leads to an answer to the first one as well. Recent years have seen remarkable progress in tackling these questions. Of particular usefulness has been the application of ideas from quantum information science to quantum many-body physics. Motivated by these developments, in this review we will focus on the second question above from an information theoretic point of view, and discuss some of the advances on the first question in section 7.

The usefulness of entanglement entropy (EE) for understanding the physics of QSLs becomes apparent when one introduces the notion of 'short-range' and 'long-range' entanglement. The heuristic picture of a short-range entangled phase is that its caricature ground-state wavefunction can be written as a direct product state in the real space, hence the nomenclature. The phases of matter that do not satisfy this property are thereby called 'long-range entangled'. From this definition, gapless phases, such as a Fermi liquid or a superfluid, that exhibit algebraic decay of local operators are always long-range entangled in the real space. Naively, one might think that the converse would hold; that is, existence of short-range correlations for local operators will imply short-range correlation. It is precisely here that the connection between gapped QSLs and ground-state entanglement properties comes into play. To wit, gapped QSLs are phases of matter that show long-range entanglement in their ground states despite the fact that the correlation functions of all local operators decay exponentially with distance.

To make the connection between gapped QSLs and long-range entanglement precise, let us define some basic terminology that will also be useful throughout this review. Given a normalized wavefunction, |ϕ〉, and a partition of the system into subsystems A and B, one can trace out the subsystem B to give a density matrix on A:ρ_A = Tr_B|ϕ〉〈ϕ|. The Renyi entropies are defined by

$$\begin{equation} S_n = \frac1{1-n}\log({\rm Tr}\rho_A^n). \end{equation} \tag{ 1 }$$

It is common to pay special attention to the von Neumann entropy, S₁ = −Tr[ρ_A log ρ_A] (obtained by taking the limit n → 1 in equation (1)). However, the Renyi entropies are also equally good measures of entanglement, and are often easier to compute.

In seminal papers by Kitaev and Preskill [18] and Levin and Wen [19], it was shown that the EE of a topologically ordered phase, such as a gapped QSL, as a function of the linear size of disc-shaped region A of linear size l, receives a universal subleading contribution:

$$\begin{equation} S(L) = \frac{l}{a} - b_0 \gamma + O(1/l), \end{equation} \tag{ 2 }$$

where b₀ is the number of connected components of the boundary of region A and the constant γ is a universal property of the phase of matter that is non-zero if and only if the phase is topologically ordered. One may wonder why does a non-zero γ imply long-range entanglement? One way to understand this [19] is to consider the combination S₀ = (S₁ − S₂) − (S₃ − S₄) for the geometry shown in figure 1. If the phase was in fact short-range entangled, then S₀ will vanish identically. From equation (2), S₀ = −2γ and, therefore, a non-zero γ implies long-range entanglement. A different way to derive the same conclusion is to define a quantity S_local such that it is obtained by locally patching the contributions to the EE from the boundary ∂A of region A [20]. If the phase is short-range entangled, then S_local is the entire S. As shown in [20], on general grounds, S_local cannot contain a constant term, and therefore, γ ≠ 0 implies long-range entanglement.

Zoom In Zoom Out Reset image size

Gapless spin liquids pose a different set of questions. The gapless excitations make it difficult to assign them a simple topological number such as γ above. One exception is gapless spin-liquids where the emergent gauge degrees of freedom are gapped and are coupled to gapless matter. In this case, there are known contributions to the EE that can characterize the nature of QSL [21, 22]. In this review, we will focus on cases where the emergent gauge fields are instead gapless. As an example, consider a two-dimensional (2D) QSL where the gapless matter fields are neutral fermions which form a Fermi surface, and are coupled to a gapless gauge field. As we discuss in section 6, empirically one finds that their EE scales as S ∼ l log l, akin to gapless fermions in two dimensions [23, 24], thus suggesting the presence of emergent fermions in a purely bosonic state.

One of the most interesting properties of EE is that $S(A) = S(\overline {A})$, where $\overline {A}$ is the complement of region A. This suggests that EE is a property of the boundary ∂A of region A since $\partial A = \partial \overline {A}$, which further suggests that for a short-range correlated system, EE scales linearly with the boundary ∂A: S(l) ∼ l^d−1, where l is the linear size of region A [25, 26]. Such a scaling of EE goes by the nomenclature of 'area law' or 'boundary law'. Indeed, for one-dimensional (1D) systems with a gap, this statement can be proved rigorously [27] and has also been observed numerically for various gapped 2D systems [28]. Furthermore, even when the correlations of local operators are long-range, for example, in a free fermionic system with a Fermi surface, the boundary law is violated only by a multiplicative logarithm, that is, S(l) ∼ l^d−1log(l) [23, 24], while for a conformal field theory (CFT), where the low-lying modes exist at only discrete point(s) in the momentum space, the boundary law is not violated at all in d > 1 [29].

To give a flavor of the information theoretic approach to many-body physics, let us consider one of the earliest such applications that concerns the relation between the central charge c of a 1D CFT, and the coefficient of the logarithmic term in the bipartite EE S for the corresponding CFT [30]. The result is

$$\begin{equation} S(l) = \frac{c}{6} \log\left(\frac{l}{a}\right), \end{equation} \tag{ 3 }$$

where l is the size of the subsystem A and a is a microscopic cutoff. The main observation to gather from this equation is that given the ground-state wavefunction of a CFT, its central charge can be extracted by calculating the EE. Note that rescaling l → l/α where α is a non-universal number does not change the coefficient of the logarithm. This makes physical sense, since c is a property of the low-energy theory and should be insensitive to the details of how the theory is regularized at the level of the lattice.

An extension of this result to a 2D CFT is [29]:

$$\begin{equation} S(l) = \frac{l}{a}- \Gamma , \end{equation} \tag{ 4 }$$

for a region A with linear size l and smooth boundaries. Here Γ is a constant that depends only on the shape of region A. Rescaling l → l/α does not change Γ akin to the case of 1D CFT, albeit in this case Γ will generically be shape dependent. As we will see, this result will be useful in characterizing gapless QSLs whose low-energy theory is a CFT (section 6).

The common feature of the above two examples is that EE as a function of system size l contains terms that are universal properties of the low-energy description of the system. Returning to our topic of interest, namely QSLs, one of the most interesting results [18, 19] is that for gapped systems in two dimensions, the subleading term in EE is universal and contains information regarding the presence or absence of topological order, and also partially characterizes the topological order:

$$\begin{equation} S(l) = \frac{l}{a} - \gamma + O(1/l). \end{equation} \tag{ 5 }$$

The subleading constant γ is called 'topological entanglement entropy' (TEE) (table 1) and, as already mentioned in the introduction, is a universal property of the phase of the matter, even though one happens to have access to only one of the ground states in the whole phase diagram. γ is non-zero if and only if the system has topological order. Furthermore, γ = log(D) where $D = \sqrt {\sum _i d^2_i}$ is the so-called total quantum dimension associated with the phase of matter whose quasi-particles carry a label d_i, the individual quantum dimension for the ith quasi-particle. As a concrete example, for Kitaev's Toric code model [32], there are four quasi-particles, the identity particle 1, the electric charge e, the magnetic charge m and the dyon em. The quantum dimension is one for each of these quasi-particles and hence the total quantum dimension $D = \sqrt {4} = 2$. We note that the universal contribution to EE was first observed in the context of discrete gauge theories in [31].

Table 1. TEE γ_calculated of spin-liquid states obtained using the Monte Carlo technique discussed in the main text. γ_expected is the exact value expected in the limit ξ/l → 0, where ξ is the correlation length and l is the subsystem size.

State	Expected γ	γcalculated/γexpected
Unprojected (ν=1)	0	−0.0008±0.0059*
CSL LA=3	$\log \sqrt {2}$	0.99±0.03
CSL LA=4	$\log \sqrt {2}$	0.99±0.12
Lattice ν=1/3	$\log \sqrt {3}$	1.07±0.05
Lattice ν=1/4	$\log \sqrt {4}$	1.06±0.11

Equation (5) serves as a definitive test to see if a particular ground state corresponds to topological ordered state or not. If the answer is indeed in the affirmative and if the state also possesses a global SU(2) spin–rotation symmetry, then it is a gapped QSL. Having established that the wavefunction corresponds to QSL, it is imperative to further characterize the state since the total quantum dimension D is only a partial characterization of topological order. For example, for D = 2, there exist two distinct QSLs, the Z₂ QSL, which has the same topological order as Kitaev's Toric code [32], and the doubled semionic theory [33], which can be thought of as a composite of the Laughlin ν = 1/2 quantum Hall state and its time-reversed partner. To achieve a detailed characterization of topological order, TEE again turns out to be a very useful concept. As we will show in section 5.2.1, TEE can be used to extract not only the individual quantum dimensions d_i for the quasi-particles, but their braiding statistics as well.

An elegant, different approach to a more complete identification of topological order is through the study of the entanglement spectrum [34]. This method is extremely useful for phases that have a gapless boundary such as quantum Hall and topological insulators/superconductors, although it may not be applicable for topological phases such as the Z₂ spin liquid that do not have edge states.

A QSL can be viewed as a state where each spin forms a singlet with a near neighbor, but the arrangement of singlets fluctuates quantum mechanically, so it is a liquid of singlets [1, 2]. Theoretical models of this singlet liquid fall roughly into two categories. In the first, the singlets are represented as microscopic variables as in quantum dimer and related models [35–39], and are suggested by large N calculations [5]. Topological order can then be established by a variety of techniques including exact solutions. In contrast, there has been less progress in establishing topological order in the second category, which are SU(2) symmetric spin systems where valence bonds are emergent degrees of freedom. More than two decades ago, Anderson [3] proposed constructing an SU(2) symmetric spin liquid wavefunction by starting with a Bardeen–Cooper–Schrieffer (BCS) state, derived from the mean-field Hamiltonian:

$$\begin{equation} H = -\sum_{rr'} \{t_{rr'}f^{\dagger}_{\sigma,r} f_{\sigma,r'} + \Delta_{r r'} f^{\dagger}_{\uparrow,r} f^{\dagger}_{\downarrow,r'}\} + \mathrm{h.c.} \end{equation} \tag{ 6 }$$

and Gutzwiller projecting [40] it so that there is exactly one fermion per site, hence a spin wavefunction. Variants of these are known to be good variational ground states for local spin Hamiltonians (see, e.g., [41, 42]) and serve as a compact description of most experimental and SU(2) symmetric liquids. Approximate analytical treatments of projection, which include small fluctuations about the above mean-field state, indicate that at least two kinds of gapped spin liquids can arise: chiral spin liquids (CSLs) [43] and Z₂ [5, 44, 45] spin liquids. The low-energy theory of CSL corresponds to a Chern–Simons theory, while in the Z₂ spin liquid, it corresponds to a deconfined phase of the Z₂ gauge theory. However, given the drastic nature of projection, it is unclear if the actual wavefunctions obtained from this procedure are in the same phase as the ground state of a Chern–Simons theory/Z₂ gauge theory. In fact, detecting deconfinement in a gauge theory in the presence of matter field is considered notoriously hard [46]. As we will see below, quantum entanglement provides an elegant solution to this problem.

In contrast to gapped QSLs, the low-energy description of gapless QSLs consist of strongly interacting matter-gauge theories which are harder to analyze [17]. Remarkably, the above procedure of constructing mean-field spin-liquid states and using Gutzwiller projection to obtain wavefunctions can be applied to obtain gapless QSLs as well. To do so, one sets the BCS gap Δ_rr' in equation (6) to zero. This Gutzwiller projection technique is known to work well in one dimension where the projected Fermi sea spin wavefunction captures long distance properties of the Heisenberg chain, and is even known to be the exact ground state of the Haldane–Shastry [48] Hamiltonian. Although similar rigorous results are not available in two dimensions, a critical spin liquid, the spinon Fermi sea (SFS) state, has been invoked [41, 49] to account for the intriguing phenomenology of aforementioned triangular lattice organic compounds [50]. In the SFS state, the spinons f hop on the triangular lattice sites giving rise to a Fermi sea, while strongly interacting with an emergent 'electromagnetic' U(1) gauge field. The metal-like specific heat and thermal conductivity seen in these materials is potentially an indication of the spinon Fermi surface. In this regard, the Gutzwiller projected Fermi sea is known to have excellent variational energy for the J₂ − J₄ spin model on the triangular lattice, which is believed to be appropriate for the aforementioned triangular lattice compounds [41].

As mentioned in the introduction, the main task of this review is to illustrate the extraction of the universal properties of a QSL from the ground-state wavefunction. Therefore, we now turn to the technical challenge of calculating EE of Gutzwiller projected wavefunctions that are candidates for QSL.

In this section, using a Monte Carlo scheme, we outline how to obtain the EE of a ground-state function that has the form of a Gutzwiller projected Slater determinant. As mentioned in the previous sections, such wavefunctions are a natural candidate for QSLs. We will focus on the second Renyi entropy S₂ (see equation (1)) since it suffices to establish topological order and is easiest to calculate within our scheme.

Let ϕ(a,b) be the wavefunction of interest, where a (b) is the configuration of subsystem A (B). Then the reduced density matrix is $\rho _A(a,a')=\sum _b\phi ^*(a,b)\phi (a',b)/N$ where the normalization $N=\sum _{a,b}\phi ^*(a,b)\phi (a,b)$. The second Renyi entropy S₂ is

$$\begin{eqnarray} \mathrm{e}^{-S_2} & = & {\rm tr} \rho_A^2 \nonumber \\ & = & \frac1{N^2}\sum_{aa'bb'} \phi^*(a,b)\phi(a',b) \phi^*(a',b')\phi(a,b'). \end{eqnarray} \tag{ 7 }$$

One way to interpret this expression is to consider two copies of the system, in configuration |a,b〉|a',b'〉. Then, the product wavefunction $|\Phi \rangle = \sum _{a,b,a',b'} \phi (a,b)\phi (a',b')|a,b\rangle |a',b'\rangle$ appears in the expression above. In addition, if we define the Swap_A operator, following [51], which swaps the configurations of the A subsystem in the two copies: Swap_A|a,b〉|a',b'〉 = |a',b〉|a,b'〉, then Renyi entropy can be written simply as

$$\begin{equation} \mathrm{e}^{-S_2} =\frac{ \langle \Phi|{\rm Swap}_A |\Phi\rangle}{\langle \Phi|\Phi\rangle}. \end{equation} \tag{ 8 }$$

The above expression suggests a variational Monte Carlo (VMC) algorithm for the Renyi entropy.

Defining configurations α₁ = a, b, α₂ = a', b' and β₁ = a', b, β₂ = a, b', and ϕ_α1 = ϕ(a,b), etc, equation (8) can be rewritten in the suggestive form

$$\begin{equation} \left\langle {\rm Swap}_{A}\right\rangle = {\sum_{\alpha_{1}\alpha_{2}}}\rho_{\alpha_{1}}\rho_{\alpha_{2}}f\left(\alpha_{1},\alpha_{2}\right). \end{equation} \tag{ 9 }$$

Here the weights $\rho _{\alpha _{i}}= |\phi _{\alpha _i}|^2/{\sum _{\alpha _i}|\phi _{\alpha _i}|^2}$ are non-negative and normalized. The quantity $f(\alpha _{1},\alpha _{2})=\frac {\phi _{\beta _1}\phi _{\beta _2}}{\phi _{\alpha _1}\phi _{\alpha _2}}$ when averaged over Markov chains generated with the probability distribution ρ(α₁) ρ(α₂) yields $\langle {\mathrm { Swap}}_{A}\rangle$. The Gutzwiller projected wavefunctions to describe the QSLs we consider in this review can be written as products of determinants, as explained in the previous section. Such wavefunctions can be efficiently evaluated, which forms the basis for the VMC technique [52].

As shown in [53, 54], the above algorithm, with slight modifications that we will discuss briefly in section 6, can be used to calculate S₂ within a few per cent of error bar for subsystems with linear size L_A ≲ 10, embedded in systems with linear size L ≲ 20.

Model systems such as the Toric code [31] and quantum dimer models [55] provide a testing ground for the validity of equation (5). More realistic interacting quantum systems pose a greater challenge since the analytical methods are limited. Some of the earlier numerical work to extract TEE included exact diagonalization studies on small systems, looking at quantum Hall Laughlin states [56] and perturbed Kitaev Toric code models [57]. In the context of QSLs, a recent quantum Monte Carlo study [58] on a sign problem free Hamiltonian used TEE to detect Z₂ topological order. This was a positive definite wavefunction, with U(1) rather than SU(2) spin symmetry. In the rest of this section, we will focus on SU(2) symmetric spin liquids [54], using the technique explained in the previous section.

5.1. Establishing topological order in gapped QSLs

In this subsection, we discuss a few examples where by using the algorithm to calculate the EE discussed in the previous section, one can establish topological order in SU(2) symmetric QSLs [54]. As a first example, let us consider an SU(2) symmetric lattice wavefunction, the Kalmeyer–Laughlin CSL.

The CSL is a spin SU(2) singlet ground state that breaks time reversal and parity symmetry [43, 59]. A wavefunction in this phase wavefunction may be obtained using the slave-particle formalism of Gutzwiller projecting a d + id BCS state [43]. Alternatively, it can be obtained by Gutzwiller projection of a hopping model on the square lattice (figure 2). This model has fermions hopping on the square lattice with a π flux through every plaquette and imaginary hoppings across the square lattice diagonals:

$$\begin{eqnarray} &&H={\sum_{\left\langle ij\right\rangle}}t_{ij}f_{i}^{\dagger}f_{j}+ \mathrm{i}{\sum_{\left\langle\left\langle ik\right\rangle \right\rangle}}\Delta_{ik}f_{i}^{\dagger}f_{k}. \end{eqnarray} \tag{ 10 }$$

Here i and j are nearest neighbors and the hopping amplitude t_ij is t along the $\hat y$ direction and alternating between t and −t in the $\hat x$-direction from row to row; and i and k are second nearest neighbors connected by hoppings along the square lattice diagonals, with amplitude Δ_ik = iΔ along the arrows and Δ_ik = −iΔ against the arrows, see figure 1. The unit cell contains two sublattices A and B. This model leads to a gapped state at half filling and the resulting valence band has unit Chern number. This hopping model is equivalent to a d + id BCS state by an $SU(2)$ gauge transformation [60]. We use periodic boundary conditions throughout this section.

Zoom In Zoom Out Reset image size

Due to the fact that this Hamiltonian contains only real bipartite hoppings and imaginary hoppings between the same sublattices and preserves the particle–hole symmetry, this wavefunction $\left \langle \alpha |\Phi \right \rangle$ can be written as a product of two Slater determinants $\mathcal M {\mathrm { Det}}(M_{ij})^2$, where $\mathcal M$ is just an unimportant Marshall sign factor, and

$$\begin{equation} M_{ij} = \left\{\left[\psi_A(k_i)+\psi_B(k_i)\right] +(-1)^{y_j}\left[\psi_A(k_i)-\psi_B(k_i)\right]\right\}\mathrm{e}^{\mathrm{i}k_i \cdot r_j}. \end{equation} \tag{ 11 }$$

Here ψ_A(k) (ψ_B(k)) is the wavefunction amplitude on sublattice A (B), r_j is the coordinates of the up spins in configuration α and k_i is the momentums in the momentum space. The Renyi entropy S₂ of this wavefunction can be calculated by the VMC method detailed in the last section.

Since TEE γ is a subleading term in EE, naively it is rather difficult to extract since the leading boundary law term dominates the total entropy. However, as shown by Kitaev and Preskill [18], one can divide the total system into four regions A,B,C,D and consider the following linear combination of entropies that completely cancels out the leading non-universal term:

$$\begin{equation} -\gamma = S_{2,A}+S_{2,B}+S_{2,C}-S_{2,AB}-S_{2,AC}-S_{2,BC} +S_{2,ABC}, \end{equation} \noindent \tag{ 12 }$$

where $S_{2,\cal {R}}$ is the Renyi entropy S₂ corresponding to the region $\cal {R}$. Note that the above formula neglects terms of O(e^−l'/ξ), where ξ is the correlation length and l' = min(l_A,l_B,l_C). This is because the above combination of entropies cancels out all local contributions to the EE which are expected to form a Taylor series expansion in the variable ξ/l', following the arguments in [20]. One could still have terms that are non-perturbative in ξ/l' which do not cancel out, hence the exponential dependence of the error made. This has also been confirmed by a direct calculation on select models in [47]. Therefore, to accurately determine TEE γ, it is important that ξ ≪ l'. Since we have a variational parameter Δ/t at our disposal, one may tune it to minimize the correlation length ξ. Choosing an optimal value Δ = 0.5t and dividing the total system with dimensions 12 × 12 lattice spacings into L_A × L_A squares A and B, and an L_A × 2L_A rectangle C (see figure 3), one finds γ = 0.343 ± 0.012 for an L_A = 3 system and γ = 0.344 ± 0.043 for an L_A = 4 system [54]. Both are in excellent agreement with the expectation of $\gamma =\log (\sqrt {2}) = 0.347$ for its ground states' two-fold degeneracy.

Zoom In Zoom Out Reset image size

To confirm the role of Gutzwiller projection in obtaining the QSL, one may also study the BCS wavefunction used to obtain CSL without Gutzwiller projection. This is an exactly solvable problem and one indeed finds a value ≈0 for TEE. This also serves as a benchmark for the VMC calculation. For an L_A = 3 system, the VMC calculation gives γ = −0.0008 ± 0.0059, in agreement with vanishing TEE [54].

One may also study the effect of increasing correlation length on the TEE results. As mentioned above, as correlation length increases, the finite size effects from subleading terms become more important. Indeed, as seen in figure 4, the TEE calculated using VMC deviates monotonically from the universal value of $\gamma =\log (\sqrt {2})$ as we lower the gap size Δ/t for a system with the typical length scale L_A = 3.

Zoom In Zoom Out Reset image size

As mentioned above, the CSL ground state can be written as ∼Det(M_ij)², where M_ij is the Slater determinant corresponding to a Chern insulator. This leads to a generalization to wavefunctions that are cube or the fourth power of a Chern insulator wavefunction. For example, consider the wavefunction

$$\begin{equation} \Psi_{1/3}(r_1,\,r_2\ldots,\,r_N)=\Phi^3(r_1,\ldots, \,r_N), \end{equation} \tag{ 13 }$$

where Φ is the Chern insulator Slater determinant defined above. Clearly, the product is a fermionic wavefunction, since exchanging a pair of particles leads to a change of sign. This is similar in spirit to constructing the corresponding Laughlin liquid of m = 3 of fermions, by taking the cube of the Slater determinant wavefunction in the lowest Landau level

$$\begin{equation} \psi(z_1,\ldots,\,z_N)=\prod_{i<j}(z_i-z_j)\mathrm{e}^{-\sum_i\frac{|z_i|^2}{4l^2_B}}. \end{equation} \tag{ 14 }$$

However, unlike the canonical Laughlin state, composed of the lowest Landau level states, these are lattice wavefunctions. An interesting question is whether the lowest Landau level structure is important in constructing states with the topological order of the Laughlin state or whether bands with identical Chern number are sufficient, as suggested by field theoretic arguments.

Calculating the TEE using VMC with L_A = 3, one finds γ = 0.5894 ± 0.0272 for the m = 3 wavefunction, consistent with the ideal value $\gamma =\log(\sqrt{3})=0.549$. Similarly, for the fourth power of the Chern insulator Slater determinant, one finds γ = 0.732 ± 0.076, again consistent with the ideal value in the thermodynamic limit: $\gamma =\log(\sqrt{4})=0.693$.

These results offer direct support for the TEE formula γ = logD as well as their validity as topological ground-state wavefunctions carrying fractional charge and statistics. The lattice fractional quantum Hall wavefunctions discussed here may be relevant to the recent studies of flat band Hamiltonians with fractional quantum Hall states [61–64].

5.2. Beyond topological entanglement entropy: extracting quasi-particle statistics using entanglement

At the fundamental level, topologically ordered phases in two dimensions display a number of unique properties such as non-trivial statistics of emergent excitations. For example, in Abelian topological phases, exchange of identical excitations or taking one excitation around another (braiding) leads to characteristic phase factors, which implies that these excitations are neither bosonic nor fermionic. A further remarkable generalization of statistics occurs in non-Abelian phases where excitations introduce a degeneracy. An important and interesting question is whether the ground state directly encodes this information and if so how one may access it? It is well known that topologically ordered phases feature a ground-state degeneracy that depends on the topology of the space on which they are defined. Also, as discussed in the previous sections, the ground-state wavefunctions of such states contain a topological contribution to the quantum entanglement, the TEE. In this section, we show that combined together, these two ground-state properties can be used to extract the generalized statistics associated with excitations in these states.

The generalized statistics of quasi-particles is formally captured by the modular $\mathcal S$ and $\mathcal {U}$ matrices, in both Abelian and non-Abelian states [65–69]. The element $\mathcal {S}_{ij}$ of the modular $\mathcal S$ matrix determines the mutual statistics of the ith quasi-particle with respect to the jth quasi-particle while the element $\mathcal {U}_{ii}$ of (diagonal) $\mathcal {U}$ matrix determines the self-statistics ('topological spin') of the ith quasi-particle. Note that these provide a nearly complete description of a topologically ordered phase—for instance, fusion rules that dictate the outcome of bringing together a pair of quasi-particles are determined from the modular $\mathcal S$ matrix, by the Verlinde formula [70]. Previously, Wen proposed [65] to use the non-Abelian Berry phase to extract statistics of quasi-particles. However, the idea in  [65] requires one to have access to an infinite set of ground states labeled by a continuous parameter, and is difficult to implement. Recently, Bais et al [71] also discussed extracting $\mathcal S$ matrix in numerical simulations, by explicit braiding of excitations. In contrast, here we will just use the set of ground states on a torus, to determine the braiding and fusing of gapped excitations.

It is sometimes stated without qualification that TEE is a quantity solely determined by the total quantum dimension D of the underlying topological theory. However, this holds true only when the boundary of the region A consists of topologically trivial closed loops. If the boundary of region A is non-contractible, for example if one divides the torus into a pair of cylinders, generically the EE is different for different ground states (see figure 5). Indeed as shown in [72] for a class of topological states, the TEE depends on the particular linear combination of the ground states when the boundary of region A contains non-contractible loops. We will exploit this dependence to extract information about the topological phase beyond the total quantum dimension D.

Zoom In Zoom Out Reset image size

Let us briefly summarize the key ideas involved in extracting braiding statistics using quantum entanglement [73] and then discuss the details. We recall that the number of ground states N on a torus corresponds to the number of distinct quasi-particle types [44, 65].⁴ Intuitively, different ground states are generated by inserting appropriate fluxes 'inside' the cycle of the torus, which is only detected by loops circling the torus. We would like to express these quasi-particle states as a linear combination of ground states. A critical insight is that this can be done using the TEE for a region A that wraps around the relevant cycle of the torus. To achieve this, we introduce the notion of minimum entropy states (MESs) Ξ_j, j = 1 − N, namely the ground states with minimal EE (or maximal TEE, since the TEE always reduces the entropy) for a given bipartition. These states are generated by the insertion of a definite quasi-particle into the cycle enclosed by region A. With this in hand, one can readily access the modular $\mathcal S$ and $\mathcal U$ matrices as we now describe.

5.2.1. Extracting statistics from TEE

As mentioned above, for Abelian phases, the ijth entry of the $\mathcal {S}$ matrix corresponds to the phase the ith quasi-particle acquires when it encircles the jth quasi-particle. The $\mathcal {U}$ matrix is diagonal and the iith entry corresponds to the phase the ith quasi-particle acquires when it is exchanged with an identical one. Since the MESs are the eigenstates of the non-local operators defined on the entanglement cut, the MESs are the canonical basis for defining $\mathcal S$ and $\mathcal U$. The modular matrices are just some unitary transformations of the MES basis. For example,

$$\begin{equation} \mathcal{S}_{\alpha\beta}=\frac{1}{D}\left\langle \Xi_{\alpha}^{\hat{x}}|\Xi_{\beta}^{\hat{y}}\right\rangle. \end{equation} \tag{ 15 }$$

Here D is the total quantum dimension and $\hat {x}$ and $\hat {y}$ are two directions on a torus. Equation (15) is just a unitary transformation between the particle states along different directions. In the case of a system with square geometry, the $\mathcal S$ matrix acts as a π/2 rotation on the MES basis $|\Xi _{\beta }^{\hat {y}} \rangle$. In general, however, $\hat {x}$ and $\hat {y}$ do not need to be geometrically orthogonal, and the system does not need to be rotationally symmetric, as long as the loops defining $|\Xi _{\alpha }^{\hat {x}} \rangle$ and $|\Xi _{\beta }^{\hat {y}} \rangle$ interwind with each other. We should mention one caveat with equation (15). When we write equation (15), we assume that the quasi-particles transform trivially under rotations. However, there exist topologically ordered states where rotations can induce a change in the quasi-particle type. For example, in the Wen–plaquette model [74], the rotations convert a charge-type excitation to a flux-type excitation. In such topological states, the procedure outlined above can give incorrect results, and we exclude such states from our discussion.

More concretely, let us denote the set of ground-state wavefunctions as |ξ_α〉, where α = 1 − N and N is the total ground-state degeneracy. To obtain the modular S matrix, one calculates the MESs |Ξ_1,α〉 and |Ξ_2,α〉 for entanglement cuts along the directions $\vec {w}_1$ and $\vec {w}_2$, respectively. Doing so leads one to the matrices U₁ and U₂ that relates the states |ξ_α〉 to |Ξ_1,α〉 and |Ξ_2,α〉, respectively. The modular $\mathcal {S}$ matrix is given by U⁻¹₂U₁ up to an undetermined phase for each MES. The existence of an identity particle that obtains a trivial phase encircling any quasi-particle helps to fix the relative phase between different MESs, requiring the entries of the first row and column to be real and positive. This completely defines the modular $\mathcal {S}$ matrix. Therefore, the modular $\mathcal {S}$ matrix can be derived even without any presumed symmetry of the given wavefunctions.

The above algorithm is able to extract the modular transformation matrix $\mathcal {S}$ and hence braiding and mutual statistics of quasi-particle excitations just using the ground-state wavefunctions as an input. Further, there is no loss of generality for non-Abelian phases, which can be deal with by enforcing the orthogonality condition in step 2 which guarantees that one obtains states with quantum dimensions d_α in increasing order. As shown in [73], if the lattice has 2π/3 rotation symmetry, then the operation of rotating MES states by an angle of 2π/3 corresponds to the action of $\mathcal {US}$ on MES. Therefore, in the presence of such a symmetry, one can extract the $\mathcal {U}$ matrix as well.

5.2.2. Example: semionic statistics in chiral spin liquid from EE

Let us revisit CSL, first discussed in section 5.1. Even though topological properties of CSL are well established using field-theoretic methods [65], from a wavefunction point of view, CSL cannot be dealt with analytically unlike continuum Laughlin states. For example, the low-energy theory of CSL predicts a particle with semionic statistics (that is, the particle picks up a factor of −1 when it goes 2π around an identical one), which is difficult to verify in the lattice wavefunction directly. We will show that the entanglement-based algorithm discussed in the last subsection readily demonstrates the existence of a semion as an excitation in CSL.

As discussed in section 5.1, CSL corresponds to a d + id-wave state for chargeless spinons. The corresponding ground-state wavefunction is obtained by Gutzwiller projecting the mean-field BCS ground-state wavefunction. CSL has a topological ground-state degeneracy of two on torus. One of the ground states can be obtained by enforcing periodic boundary (figure 6) condition along the $\hat {x}$-direction and anti-periodic along $\hat {y}$. We denote it by |0,π〉. The orthogonal ground state is given by |π,0〉 with the same notation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our task is to obtain the modular $\mathcal {S}$-matrix since it encodes the mutual statistics of quasi-particles. To do so, as outlined in the previous subsection, we consider linear combinations of $\left |0,\pi \right \rangle$ and $\left | \pi ,0\right \rangle$:

$$\begin{equation} \left|\Phi(\phi)\right\rangle =\cos\phi\left|0,\pi\right\rangle +\sin\phi\left|\pi,0\right\rangle \end{equation} \tag{ 16 }$$

and numerically minimize the EE with respect to the parameter ϕ so as to obtain the MESs [73]. For a system with square geometry the $\mathcal {S}$ matrix describes the action of π/2 rotation on the MESs. Since the two states |0,π〉 and |π,0〉 transform into each other under π/2 rotation, this implies that the $\mathcal {S}$ matrix is given by V^†σ_xV , where V is the unitary matrix that rotates the {|0,π〉 |π,0〉} basis to the MESs basis. From figure 8, one finds that

$$\begin{eqnarray*} V \approx \left(\begin{array}{@{}cc@{}} \cos(0.14 \pi) & -\sin(0.14 \pi )\mathrm{e}^{\mathrm{i}\varphi} \\ \sin(0.14\pi) & \cos(0.14\pi )\mathrm{e}^{\mathrm{i}\varphi} \end{array} \right). \end{eqnarray*}$$

The existence of an identity particle requires positive real entries in the first row and column and implies φ = 0, which gives

$$\begin{eqnarray*} \mathcal{S} \approx \left(\begin{array}{@{}cc@{}} 0.77 & 0.63 \\ 0.63 & -0.77 \end{array}\right) = \frac{1}{\sqrt{2}}\left(\begin{array}{cc} 1.08 & 0.89 \\ 0.89 & -1.08 \end{array}\right), \end{eqnarray*}$$

which is in good agreement with the exact result

$$\begin{equation} \mathcal{S}_{{\rm exact}} = \frac{1}{\sqrt{2}}\left(\begin{array}{@{}cc@{}} 1 & 1 \\ 1 & -1 \end{array}\right). \end{equation} \tag{ 17 }$$

Even though the $\mathcal {S}$ matrix obtained using our method is approximate, the quasi-particle statistics can be extracted exactly. In particular, the above $\mathcal {S}$ matrix tells us that the quasi-particle corresponding to d₀ = 1 does not acquire any phase when it goes around any other particle and corresponds to an identity particle as expected, while the quasi-particle corresponding to d_1/2 = 1 has semion statistics since it acquires a phase of π when it encircles another identical particle. Numerical improvements can further reduce the error in pinpointing the MES and thereby lead to a more accurate value of the $\mathcal S$ matrix.

Zoom In Zoom Out Reset image size

It is also interesting to compare the full dependence of TEE on the angle ϕ. Theoretically, one finds the following expression for TEE [72]:

$$\begin{eqnarray} 2\gamma-\gamma' &=&\log\frac{4}{3+\sin\left(4\phi\right)}, \end{eqnarray} \tag{ 18 }$$

where γ is the TEE for a region with contractible boundary, while γ' is that for a region with non-contractible boundary. From figure 8, we see that the numerically evaluated TEE is in rather good agreement with the exact expression.

In this section, we consider Gutzwiller projected Fermi liquid wavefunctions, which are considered good ansätze for ground states of critical spin liquids [41]. We analyze two different classes of critical spin liquids: states that at the slave-particle mean-field level have a full Fermi surface of spinons and those with only nodal spinons. The basic idea we use here is that for free fermions, the boundary law is violated by a multiplicative logarithm, that is, S(l) ∼ l^d−1log(l) [23, 24]. This suggests that for critical spin liquids, whose low-energy theory contains a spinon Fermi surface [17], a similar result might be expected. Indeed, one finds that for such Gutzwiller projected wavefunctions the EE scales as S(l) ∼ l log(l) for small system sizes [53]. On the other hand, for wavefunctions with nodal spinons, the low-energy theory is believed to be a 2D CFT for which the result in equation (4) must hold. Thus, in this case, the VMC algorithm can be used to extract the universal shape-dependent subleading constant.

The wavefunctions for these spin liquids are constructed by starting with a system of spin-1/2 fermionic spinons f_rα hopping on a finite lattice at half-filling with a mean-field Hamiltonian: $H_{\mathrm { MF}} = \sum _{r r'} \left [ -t_{r r'}f^{\dagger }_{r\sigma } f_{r'\sigma } + \mathrm {h.c.} \right ]$. The spin wavefunction is given by |ϕ〉 = P_G|ϕ〉_MF, where |ϕ〉_MF is the ground state of H_MF and the Gutzwiller projector $P_{\mathrm { G}} = \prod _i \left ( 1- n_{i \uparrow } n_{i \downarrow }\right )$ ensures exactly one fermion per site. The 1D case was discussed in [75].

6.1. Critical QSL with the spinon Fermi surface

Consider the mean-field ansatz that describes spinons hopping on a triangular lattice. In this case, one might expect that the projected wavefunction describes the ground state of a QSL whose low-energy theory contains a spinon Fermi surface [17, 41]. Indeed, one finds [53] that for a triangular lattice of total size 18 × 18 on a torus and a subregion A of square geometry with linear size L_A up to eight sites, the Renyi entropy scales as L_A logL_A (figure 9). This is rather striking since the wavefunction is a spin wavefunction and therefore could also be written in terms of hard-core bosons. This result strongly suggests the presence of an underlying spinon Fermi surface. In fact, the coefficient of the L_A logL_A term is rather similar before and after projection. This observation may be rationalized as follows. Consider 1D fermions with N flavors, where the central charge (and hence the coefficient of the logL_A entropy term) is proportional to N. However, Gutzwiller projection to single occupancy reduces the central charge to N − 1. This is a negligible change for large N. Similarly, a 2D Fermi surface may be considered as a collection of many independent 1D systems in momentum space, and Gutzwiller projection removes only one degree of freedom.

Zoom In Zoom Out Reset image size

The area-law violation of the Renyi entropy for Gutzwiller projected wavefunctions substantiates the theoretical expectation that an underlying Fermi surface is present in the spin wavefunction.

6.2. Critical QSL with nodal spinons

6.3. Effect of projection on the nested Fermi surface

Finally, consider a mean-field ansatz as the fermions on a square lattice with nearest-neighbor hopping. The unprojected Fermi surface is nested. Since projection amounts to taking correlations into account, one might wonder whether the Fermi surface undergoes a magnetic instability after the projection. Indeed, one finds non-zero magnetic order in the projected wavefunction [53], consistent with an independent study [78]. This was verified by calculating spin–spin correlations on a 42 × 42 lattice. Renyi entropy results for a lattice of total size 24 × 24 with region A being a square up to size 12 × 12 are shown in figure 10. Although it is difficult to rule out the presence of a partial Fermi surface from Renyi entropy, there is a significant reduction in the Renyi entropy as well as the coefficient of the L_A logL_A term as compared to the unprojected Fermi sea.

Zoom In Zoom Out Reset image size

In this paper, we reviewed the perspective provided by the concept of EE on QSLs. In particular, we discussed the practical uses of the concept of TEE to detect topological order in realistic QSLs and a detailed characterization of topological order using entanglement as well. Returning to the questions posed in the introduction, let us discuss whether EE is of any value in finding new QSLs? Remarkably, one of the most accurate methods to study strongly correlated systems, namely, the density matrix renormalization group (DMRG) [79], employs EE as a resource to find ground states of particular Hamiltonians. In brief, DMRG implements a real-space renormalization group that successively coarse-grains the system on the basis of the entanglement between a sub-region and the rest of the system. Although DMRG is most accurate for 1D quantum systems, recent works are now beginning to explore the realm of 2D strongly correlated systems. As an output, DMRG yields an approximate ground state of the system along with the EE for a bipartition of the system into two halves. Having access to the EE implies that one can use the results reviewed in the previous sections, thus allowing one to establish (or rule out) a QSL ground state. As an example, using DMRG, Yan et al [10] found that the ground state of kagome lattice Heisenberg Hamiltonian is a QSL. In a more recent study, Jiang et al [11] established the topological nature of this spin liquid by calculating the TEE of the ground state for a region A with non-contractible boundary and found it to be close to log(2), a value that is apparently consistent with that of a Z₂ QSL. Similar results are found for other models which are known to have topologically ordered ground states [11]. However, as discussed in section 5.2.1, TEE for a region with non-contractible boundary can be any number between 0 and log(2), the maximum value of TEE (hence the minimum value for total entropy), being attained when the ground state corresponds to quasi-particle states. These results indicate that DMRG might be biased towards a linear combination of ground states that minimize total entropy. As discussed in section 5.2, entanglement can be used to determine quasi-particle statistics from a topologically degenerate set of ground states. Recently, using such an approach, Cincio and Vidal [80] extracted the modular S and U matrices for a microscopic lattice Hamiltonian that has a CSL ground state. Entanglement-based measures have therefore emerged as a powerful conceptual and numerical tool in the context of QSLs. One may expect that, in future, the confluence of ideas of many-body physics and quantum information will throw up other equally fruitful directions.

We thank Ari Turner and Masaki Oshikawa for useful discussions and collaboration on the topics in this review. Support from the National Science Foundation (NSF) under grant no. NSF DMR-0645691 and hospitality at KITP, Santa Barbara, where part of this work was carried out, are gratefully acknowledged. This research was supported in part by the NSF under grant no. NSF PHY11-25915.