Kagome model for a ${\mathbb Z}_2$ quantum spin liquid

TL;DR

This paper introduces a sign-problem-free model on the Kagome lattice that exhibits a gapped ${\mathbb Z}_2$ topological quantum spin liquid, supported by numerical evidence including entanglement entropy analysis.

Contribution

The study presents a novel, sign-problem-free Kagome lattice model that demonstrates a gapped ${\mathbb Z}_2$ quantum spin liquid phase, with insights from large-$N$ analysis and entanglement entropy.

Findings

Identification of a quantum spin liquid phase in the model

Evidence of a gapped ${\mathbb Z}_2$ topological phase

Numerical confirmation via entanglement entropy

Abstract

We present a study of a simple model antiferromagnet consisting of a sum of nearest neighbor SO($N$) singlet projectors on the Kagome lattice. Our model shares some features with the popular $S=1/2$ Kagome antiferromagnet but is specifically designed to be free of the sign-problem of quantum Monte Carlo. In our numerical analysis, we find as a function of $N$ a quadrupolar magnetic state and a wide range of a quantum spin liquid. A solvable large-$N$ generalization suggests that the quantum spin liquid in our original model is a gapped ${\mathbb Z}_2$ topological phase. Supporting this assertion, a numerical study of the entanglement entropy in the sign free model shows a quantized topological contribution.

Equations

$$H$$

$$|S_{ij}⟩$$

$$H_b$$

$$Z$$

$${\mathcal{L}}_b$$

$$+$$

$$2\gamma =S_{\text{leavevmode}to5.98pt\text{vbox}to5.41pt\text{pgfpicture}\text{makeatletter}\text{lower}-0.56905ptto0.0pt\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{definecolor}pgfstrokecolorrgb\mathrm{0,0,0}\text{pgfsys@color@rgb@stroke}000\text{pgfsys@invoke}\text{pgfsys@color@rgb@fill}000\text{pgfsys@invoke}\text{pgfsys@setlinewidth}0.4pt\text{pgfsys@invoke}\text{nullfont}to0.0pt\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{pgfsys@setlinewidth}1.1381pt\text{pgfsys@invoke}\text{pgfsys@moveto}4.83698pt4.26791pt\text{pgfsys@lineto}0.0pt4.26791pt\text{pgfsys@lineto}0.0pt0.0pt\text{pgfsys@lineto}4.83698pt0.0pt\text{pgfsys@stroke}\text{pgfsys@invoke}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{hss}\text{pgfsys@discardpath}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{hss}\text{lxSVG@closescope}\text{endpgfpicture}}^{(2)}-S_{\text{leavevmode}to5.41pt\text{vbox}to5.98pt\text{pgfpicture}\text{makeatletter}\text{lower}-1.1381ptto0.0pt\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{definecolor}pgfstrokecolorrgb\mathrm{0,0,0}\text{pgfsys@color@rgb@stroke}000\text{pgfsys@invoke}\text{pgfsys@color@rgb@fill}000\text{pgfsys@invoke}\text{pgfsys@setlinewidth}0.4pt\text{pgfsys@invoke}\text{nullfont}to0.0pt\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{pgfsys@setlinewidth}1.1381pt\text{pgfsys@invoke}\text{pgfsys@moveto}0.0pt0.0pt\text{pgfsys@lineto}4.26791pt0.0pt\text{pgfsys@lineto}4.26791pt4.26791pt\text{pgfsys@lineto}0.0pt4.26791pt\text{pgfsys@lineto}0.0pt-0.56906pt\text{pgfsys@stroke}\text{pgfsys@invoke}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{hss}\text{pgfsys@discardpath}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{hss}\text{lxSVG@closescope}\text{endpgfpicture}}^{(2)}+S_{\text{leavevmode}to5.98pt\text{vbox}to5.41pt\text{pgfpicture}\text{makeatletter}\text{lower}-0.56905ptto0.0pt\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{definecolor}pgfstrokecolorrgb\mathrm{0,0,0}\text{pgfsys@color@rgb@stroke}000\text{pgfsys@invoke}\text{pgfsys@color@rgb@fill}000\text{pgfsys@invoke}\text{pgfsys@setlinewidth}0.4pt\text{pgfsys@invoke}\text{nullfont}to0.0pt\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{pgfsys@setlinewidth}1.1381pt\text{pgfsys@invoke}\text{pgfsys@moveto}-0.56906pt0.0pt\text{pgfsys@lineto}4.26791pt0.0pt\text{pgfsys@lineto}4.26791pt4.26791pt\text{pgfsys@lineto}-0.56906pt4.26791pt\text{pgfsys@stroke}\text{pgfsys@invoke}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{hss}\text{pgfsys@discardpath}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{hss}\text{lxSVG@closescope}\text{endpgfpicture}}^{(2)}-S_{\text{leavevmode}to6.54pt\text{vbox}to5.41pt\text{pgfpicture}\text{makeatletter}\text{lower}-0.56905ptto0.0pt\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{definecolor}pgfstrokecolorrgb\mathrm{0,0,0}\text{pgfsys@color@rgb@stroke}000\text{pgfsys@invoke}\text{pgfsys@color@rgb@fill}000\text{pgfsys@invoke}\text{pgfsys@setlinewidth}0.4pt\text{pgfsys@invoke}\text{nullfont}to0.0pt\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{pgfsys@setlinewidth}1.1381pt\text{pgfsys@invoke}\text{pgfsys@moveto}-0.56906pt0.0pt\text{pgfsys@lineto}4.83698pt0.0pt\text{pgfsys@stroke}\text{pgfsys@invoke}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{pgfsys@setlinewidth}1.1381pt\text{pgfsys@invoke}\text{pgfsys@moveto}-0.56906pt4.26791pt\text{pgfsys@lineto}4.83698pt4.26791pt\text{pgfsys@stroke}\text{pgfsys@invoke}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{hss}\text{pgfsys@discardpath}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{hss}\text{lxSVG@closescope}\text{endpgfpicture}}^{(2)}.$$

$$2\gamma =\log (Z_{\text{leavevmode}to5.41pt\text{vbox}to5.98pt\text{pgfpicture}\text{makeatletter}\text{lower}-1.1381ptto0.0pt\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{definecolor}pgfstrokecolorrgb\mathrm{0,0,0}\text{pgfsys@color@rgb@stroke}000\text{pgfsys@invoke}\text{pgfsys@color@rgb@fill}000\text{pgfsys@invoke}\text{pgfsys@setlinewidth}0.4pt\text{pgfsys@invoke}\text{nullfont}to0.0pt\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{pgfsys@setlinewidth}1.1381pt\text{pgfsys@invoke}\text{pgfsys@moveto}0.0pt0.0pt\text{pgfsys@lineto}4.26791pt0.0pt\text{pgfsys@lineto}4.26791pt4.26791pt\text{pgfsys@lineto}0.0pt4.26791pt\text{pgfsys@lineto}0.0pt-0.56906pt\text{pgfsys@stroke}\text{pgfsys@invoke}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{hss}\text{pgfsys@discardpath}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{hss}\text{lxSVG@closescope}\text{endpgfpicture}}^{(2)}/Z_{\text{leavevmode}to5.98pt\text{vbox}to5.41pt\text{pgfpicture}\text{makeatletter}\text{lower}-0.56905ptto0.0pt\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{definecolor}pgfstrokecolorrgb\mathrm{0,0,0}\text{pgfsys@color@rgb@stroke}000\text{pgfsys@invoke}\text{pgfsys@color@rgb@fill}000\text{pgfsys@invoke}\text{pgfsys@setlinewidth}0.4pt\text{pgfsys@invoke}\text{nullfont}to0.0pt\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{pgfsys@setlinewidth}1.1381pt\text{pgfsys@invoke}\text{pgfsys@moveto}4.83698pt4.26791pt\text{pgfsys@lineto}0.0pt4.26791pt\text{pgfsys@lineto}0.0pt0.0pt\text{pgfsys@lineto}4.83698pt0.0pt\text{pgfsys@stroke}\text{pgfsys@invoke}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{hss}\text{pgfsys@discardpath}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{hss}\text{lxSVG@closescope}\text{endpgfpicture}}^{(2)})-\log (Z_{\text{leavevmode}to5.98pt\text{vbox}to5.41pt\text{pgfpicture}\text{makeatletter}\text{lower}-0.56905ptto0.0pt\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{definecolor}pgfstrokecolorrgb\mathrm{0,0,0}\text{pgfsys@color@rgb@stroke}000\text{pgfsys@invoke}\text{pgfsys@color@rgb@fill}000\text{pgfsys@invoke}\text{pgfsys@setlinewidth}0.4pt\text{pgfsys@invoke}\text{nullfont}to0.0pt\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{pgfsys@setlinewidth}1.1381pt\text{pgfsys@invoke}\text{pgfsys@moveto}-0.56906pt0.0pt\text{pgfsys@lineto}4.26791pt0.0pt\text{pgfsys@lineto}4.26791pt4.26791pt\text{pgfsys@lineto}-0.56906pt4.26791pt\text{pgfsys@stroke}\text{pgfsys@invoke}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{hss}\text{pgfsys@discardpath}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{hss}\text{lxSVG@closescope}\text{endpgfpicture}}^{(2)}/Z_{\text{leavevmode}to6.54pt\text{vbox}to5.41pt\text{pgfpicture}\text{makeatletter}\text{lower}-0.56905ptto0.0pt\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{definecolor}pgfstrokecolorrgb\mathrm{0,0,0}\text{pgfsys@color@rgb@stroke}000\text{pgfsys@invoke}\text{pgfsys@color@rgb@fill}000\text{pgfsys@invoke}\text{pgfsys@setlinewidth}0.4pt\text{pgfsys@invoke}\text{nullfont}to0.0pt\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{pgfsys@setlinewidth}1.1381pt\text{pgfsys@invoke}\text{pgfsys@moveto}-0.56906pt0.0pt\text{pgfsys@lineto}4.83698pt0.0pt\text{pgfsys@stroke}\text{pgfsys@invoke}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{pgfsys@setlinewidth}1.1381pt\text{pgfsys@invoke}\text{pgfsys@moveto}-0.56906pt4.26791pt\text{pgfsys@lineto}4.83698pt4.26791pt\text{pgfsys@stroke}\text{pgfsys@invoke}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{hss}\text{pgfsys@discardpath}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{hss}\text{lxSVG@closescope}\text{endpgfpicture}}^{(2)}).$$

$${\mathcal{Z}}_{AB}^{(2)}(\lambda )=\sum _{C\subseteq A-B}{\lambda }^{N_C}(1-\lambda {)}^{N_A-N_B-N_C}{Z}_{B+C}^{(2)}.$$

$$\log (Z_A^{(2)}/Z_B^{(2)})={\int }_0^{1}d\lambda \frac{\partial \log {\mathcal{Z}}_{AB}^{(2)}(\lambda )}{\partial \lambda }.$$

$$2\gamma ={\int }_0^{1}d\lambda \frac{\Delta ⟨N_C{⟩}_{\lambda }}{\lambda (1-\lambda )},$$

Full text

Quantum antiferromagnetism on the Kagome lattice is an important playground in the study of quantum spin liquids emerging from frustrated magnetism. The most popular model in this family is the $S=1/2$ Kagome anti-ferromagnet $H=J\sum_{\braket{ij}}\vec{S}_{i}\cdot\vec{S}_{j}$. Despite a quarter of a century of intense research using an array of numerical and analytic methods on this important model the ground state of the $S=1/2$ Kagome anti-ferromagnet remains hotly contested. While the absence of magnetic order is uncontroversial Marston and Zeng (1995); Chalker and Eastmond (1992); Singh and Huse (1992); Leung and Elser (1993); Lecheminant et al. (1997), various nonmagnetic ground states have been proposed including, e.g., an array of quantum spin liquids Ran et al. (2007); Yan et al. (2011); Iqbal et al. (2013); He et al. (2017) and valence bond solid ordering Marston and Zeng (1995); Nikolic and Senthil (2003); Singh and Huse (2007). In parallel to the theoretical work, a number of synthetic quantum materials have been identified that provide venues where the interplay of quantum fluctuations and frustration on the Kagome lattice give rise to novel unexplained behavior Balents (2010).

Of all the proposed phases of matter on the Kagome, the so-called gapped ${\mathbb{Z}}_{2}$ quantum spin liquid Sachdev (1992) is the simplest example of an exotic state with long range entanglement Savary and Balents (2016), a prototypical quantum state that cannot be deformed into a simple product or mean field state. In its simplest incarnation, the excitations above the ground state come in two basic varieties, an $e$ particle and an $m$ particle which by themselves are bosons but are mutual semions Wen (1991); Kitaev and Laumann (2008). Remarkably it has been shown that the presence of these excitations can be detected in the entanglement of the ground state wavefunction itself, giving rise to a contribution called the “topological entanglement entropy” Levin and Wen (2006); Kitaev and Preskill (2006). Although this state is not yet experimentally accesible, we now have a few model Hamiltonians that realize this topological order, including the toric code Kitaev (2003), the honeycomb Kitaev model Kitaev and Laumann (2008), non-bipartite quantum dimer models Moessner and Sondhi (2001); Misguich et al. (2002) and models of frustrated bosons Balents et al. (2002); Isakov et al. (2006); Dang et al. (2011). It is clearly of great interest to extend this family of models with an eye to finding simple models that could find realizations in physical systems.

Model: A number of variations on the basic $S=1/2$ Heisenberg model have been introduced and studied on the Kagome lattice, including Sp($N$) Sachdev (1992), SU($N$) Corboz et al. (2012), larger spin versions of the two spin Heisenberg exchange Changlani and Läuchli (2015), as well as certain multi-spin interactions Bauer et al. (2014). In this work we present and study a new variant of the Kagome anti-ferromagnet. Our model is constructed from spins which have a local Hilbert space of $N$ states, denoted for site $j$ as $\Ket{\alpha}_{j}$ where $\alpha=1,\ldots,N$. The Hamiltonian can be written simply as a sum of singlet projectors on the nearest neighbors of the Kagome lattice,

[TABLE]

Physically, the Hamiltonian Eq. (1) can be viewed as lowering the energy of singlet formation locally between nearest neighbors. Since all pairs of neighbors cannot simultaneously form singlets, quantum fluctuations play an important role in stabilizing the ground state. We note here that the usual $S=1/2$ Heisenberg model is also a sum of singlet projectors, of the form Eq. (1) but with $\Ket{S_{ij}}\rightarrow\frac{\Ket{\uparrow\downarrow}-\Ket{\downarrow\uparrow}}{\sqrt{2}}$, which aside from the crucial relative minus sign, is identical to our singlet Eq. (2) at $N=2$. It is this discrepancy of sign that allows us to sidestep the infamous sign problem and carry out large volume numerical studies that are so far impossible for the $S=1/2$ Kagome Heisenberg model.

The model Eqs. (1,2) has a global SO($N$) symmetry in which each site transforms in the fundamental representation, $|\alpha\rangle\rightarrow O_{\alpha\beta}|\beta\rangle$ and in the path integral can be interpreted as a statistical mechanics model of tightly packed unoriented loops Kaul (2012). A previous study Kaul (2015) on the triangular lattice found a $\sqrt{12}\times\sqrt{12}$ valence bond solid order at large values of $N$. Here by introducing a solvable large-$N$ limit and a numerical study of the entanglement entropy at finite-$N$, we show that the increased geometric frustration of the Kagome lattice realizes a ${\mathbb{Z}}_{2}$ topological quantum spin liquid.

We simulate the model Hamiltonian Eqs. (1-2) using the stochastic series expansion Sandvik (2010) with loop updates on $3\times L\times L$ lattices at an inverse temperature $\beta$. To characterize the breaking of SO($N$) symmetry we introduce the operator $\hat{Q}_{\alpha\beta}=\Ket{\alpha}\Bra{\beta}-\frac{\delta_{\alpha\beta}}{N}$ which because of its tensorial nature we will call the “quadrupolar” order parameter. The Fourier transformed susceptibility, $\chi_{Q}(\mathbf{k})=\frac{1}{\beta N_{\text{site}}}\int_{0}^{\beta}d\tau\sum_{\mathbf{r}}e^{i\mathbf{k}\cdot\mathbf{r}}\langle{Q}_{\alpha\alpha}(\mathbf{r},\tau){Q}_{\alpha\alpha}(\mathbf{0},0)\rangle,$ is used to diagnose quadrupolar order. We define $Q^{2}_{\rm zf}=\chi_{Q}({\bf 0})$ as the order parameter, a quantity which scales to a finite value in the thermodynamic limit in the quadrupolar phase and to zero otherwise. To facilitate detection of the long range order, we also study a correlation ratio $R_{\text{S}}=1-\frac{\chi_{Q}(\mathbf{G}/L)}{\chi_{Q}(\mathbf{0})}$ where ${\bf G}$ is the shortest reciprocal lattice vector, which scales to 1(0) in the symmetry broken (unbroken) phase. As shown in Fig. 1, the quadrupolar order decreases as $N$ is increased. Finite size scaling shows that for $N\leq 9$ there is quadrupolar order that breaks the SO($N$) symmetry and for $N>10$ the quadrupolar order vanishes. The $N=9$ case is on the verge of transition but a careful finite size scaling indicates that it is quadrupolar ordered. We have searched extensively for translational symmetry breaking at the $N$ for which quadrupolar order is absent (as was found in the triangular lattice Kaul (2015)), but we find no evidence for this order, indicating the possibility of a liquid like state. We now present field theoretic arguments and numerical evidence that this phase is a ${\mathbb{Z}}_{2}$ quantum spin liquid.

*Large-$N$ limit – * To introduce a solvable large-$N$ limit that can capture both the quadrupolar as well as non-magnetic phase, we generalize the spins in our model to transform under larger representations than the fundamental SO($N$), using Schwinger bosons in which each local spin state is associated with one of $N$ flavors of boson $b_{i\alpha}$ ( with $[b_{i\alpha},b^{\dagger}_{j\beta}]=\delta_{ij}\delta_{\alpha\beta}$). The generalized spin model is then,

[TABLE]

with the constraint $\sum_{\alpha}b^{\dagger}_{i\alpha}b_{i\alpha}=n_{b}$, which fixes the representation of the spin. Thus the family of models Eq. (3) has two parameters $n_{b}$ and $N$. Clearly $n_{b}=1$ corresponds to Eq. (1). Increasing $n_{b}$ is a generalization of Eq. (1), with different representations of SO($N$). These are SO($N$) analogues of the well known Schwinger boson method of implementing higher representations of SU($N$)  Arovas and Auerbach (1988); Read and Sachdev (1989a).

Using boson coherent states we obtain a functional integral representation of the partition function $Z={\rm Tr}[e^{-\beta H_{b}}]$, with a field $\lambda_{i}(\tau)$ that enforces the on-site constraint and a Hubbard-Stratonovich field $Q_{ij}(\tau)$ that decouples the quartic interaction Arovas and Auerbach (1988),

[TABLE]

where we have set $n_{b}=\kappa N$. Integrating out the $b$ fields we obtain an effective action proportional to $N$. By fixing $\kappa$ a large-$N$ limit can be accessed simply by a saddle point evaluation. Assuming space and time independent $Q$ and $\lambda$ we find evaluating the trace over bosons: $Z=e^{-\beta VNf}$ (where $V$ is the total number of spatial unit cells), where, $f=\frac{1}{V}\sum_{{\bf k}\alpha}\left(\frac{zQ^{2}}{2J}-\lambda\left(\kappa+\frac{1}{2}\right)+\frac{1}{\beta}\log 2\sinh\left(\frac{\beta\omega_{{\bf k}\alpha}}{2}\right)\right)$, and $\omega_{{\bf k}\alpha}=\sqrt{\lambda^{2}-4Q^{2}\gamma^{2}_{{\bf k}\alpha}}$ the dispersion of bosons and $\gamma_{{\bf k}\alpha}$ are the three modes $\alpha=1,2,3$ of the adjacency matrix on the Kagome, $\{1,\frac{1}{2}\left(-1\pm\sqrt{3+2{\rm cos}(\vec{k}\cdot\vec{a}_{i})}\right)\}$, where $\vec{a}_{i}$ are the three shortest lattice vectors on the triangular lattice Bravais lattice. At the saddle point (obtained by extremizing $Q$ and $\lambda$), there are two phases one where the $b_{\alpha}$ have a gap and the other where they are condensed. The condensed phase of the $b_{\alpha}$ breaks the SO($N$) symmetry and corresponds to the quadrupolar order for the spin model, the implications of the gapped phase for the spin model are more subtle – we address them below. Numerically we find the transition between these two phases where the gap goes to zero is at $\kappa_{c}\approx 0.148\dots$ This gives a phase boundary $n_{b}=\kappa_{c}N$ that we show in Fig. 2 as a solid line. Also shown in solid circles are the phases determined from QMC for $n_{b}=1$. From this figure it is plausible by continuity that the quadruplar phase found in QMC corresponds to the condensation of $b_{\alpha}$ ($\kappa>\kappa_{c}$) and the liquid like phase in the QMC corresponds to the state in which the $b_{\alpha}$ are gapped ($\kappa<\kappa_{c}$). We now ask what non-magnetic state the original spin model goes into when the $b_{\alpha}$ acquire a gap? Following previous work Read and Sachdev (1989a), the state is determined by $1/N$ fluctuations beyond mean field, which take the structure of a U(1) gauge theory. The unique aspect here is that because of the non-bipartite lattice all the $b_{i}$ carry the same sign of gauge charge (as opposed to the staggered signs on bipartite lattices), and thus because of the structure of the saddle point there is a charge-2 Higgs field coupled to the U(1) gauge theory. As originally discussed in seminal work such a Higgs phase leaves behind a $\mathbb{Z}_{2}$ gauge theory and a topological phase Fradkin and Shenker (1979). This line of argument was used previously to establish emergent $\mathbb{Z}_{2}$ gauge structures in large-$N$ expansions Read and Sachdev (1991). We thus conclude that for $\kappa<\kappa_{c}$ the spin model will be in a $\mathbb{Z}_{2}$ quantum spin liquid phase. This suggests by continuity that the liquid like phase observed in our original model, Eq. (1) [the $n_{b}=1$ limit of Eq. (3)] is also in this interesting phase. We test this conjecture below.

It is also possible to take a direct large-$N$ limit of Eq. (1) (i.e. holding $n_{b}=1$ fixed). Analogous to work on SU($N$) models on bipartite lattices Read and Sachdev (1989b) we obtain as an effective theory, a quantum dimer model on the Kagome lattice, where the spin wavefunction is obtained by replacing each dimer with $|S\rangle$ of Eq. (2). At $N=\infty$ all dimer coverings are degenerate and $1/N$ corrections introduce dynamics into the quantum dimer model. In this limit it is clear that quadrupolar order is absent, consistent with our numerical findings. While the quantum dimer model so obtained is not generally solvable, it is plausible that for $N\gg 1$, our model Eq. (1) ends up in the same ${\mathbb{Z}}_{2}$ spin liquid phase as an exactly solved Kagome quantum dimer model Misguich et al. (2002), since the topological spin liquid is expected to be stable to all small deformations of the Hamiltonian. This limit suggests that at $n_{b}=1$ the model remains in a liquid state for arbitrary large $N$, as shown in Fig. 2.

Entanglement: Having presented circumstantial evidence for a spin liquid in our model Eq. (1) from large-$N$ expansions, we return to numerical simulations to provide direct evidence for the $\mathbb{Z}_{2}$ quantum spin liquid phase. We carry out measurements of the topological entanglement entropy (TEE), which has been a fruitful tool to detect topological order numerically Furukawa and Misguich (2007); Isakov et al. (2011); Grover et al. (2013); Jiang et al. (2012). In phases with topological order the TEE appears as a universal negative contribution to the entanglement entropy Levin and Wen (2006); Kitaev and Preskill (2006). With $L$ the linear size of a smooth simply connected subsystem, for large $L$ in a thermodynamic system: $S_{L}=aL-\gamma+...$, where the first term is the so-called “area law” contribution with $a$ non-universal and the second term is the universal TEE piece. For the $\mathbb{Z}_{2}$ state found in our large-$N$ study and in the Kagome quantum dimer model, it is predicted that $\gamma=\log(2)$ in the ground state. This expectation has been extended to finite temperatures as well, where due to two different excitation gaps associated with $e$ and $m$ particles, $\gamma$ is predicted to show two plateaus at $\log(2)/2$ and $\log(2)$ as a function of inverse temperature Castelnovo and Chamon (2007).

In order to isolate $\gamma$ we compute the difference in entanglement entropy of differently shaped regions Levin and Wen (2006), written as $2\gamma=S^{(2)}_{\leavevmode\hbox to5.98pt{\vbox to5.41pt{\pgfpicture\makeatletter\hbox{\hskip 0.56905pt\lower-0.56905pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{} {{}{}}{} {}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{1.1381pt}\pgfsys@invoke{ }{}\pgfsys@moveto{4.83698pt}{4.26791pt}\pgfsys@lineto{0.0pt}{4.26791pt}\pgfsys@lineto{0.0pt}{0.0pt}\pgfsys@lineto{4.83698pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}-S^{(2)}_{\leavevmode\hbox to5.41pt{\vbox to5.98pt{\pgfpicture\makeatletter\hbox{\hskip 0.56905pt\lower-1.1381pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{} {{}{}}{} {}{} {{}{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{1.1381pt}\pgfsys@invoke{ }{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{4.26791pt}{0.0pt}\pgfsys@lineto{4.26791pt}{4.26791pt}\pgfsys@lineto{0.0pt}{4.26791pt}\pgfsys@lineto{0.0pt}{-0.56906pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}+S^{(2)}_{\leavevmode\hbox to5.98pt{\vbox to5.41pt{\pgfpicture\makeatletter\hbox{\hskip 1.1381pt\lower-0.56905pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{} {}{} {}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{1.1381pt}\pgfsys@invoke{ }{}\pgfsys@moveto{-0.56906pt}{0.0pt}\pgfsys@lineto{4.26791pt}{0.0pt}\pgfsys@lineto{4.26791pt}{4.26791pt}\pgfsys@lineto{-0.56906pt}{4.26791pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}-S^{(2)}_{\leavevmode\hbox to6.54pt{\vbox to5.41pt{\pgfpicture\makeatletter\hbox{\hskip 1.1381pt\lower-0.56905pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{1.1381pt}\pgfsys@invoke{ }{}\pgfsys@moveto{-0.56906pt}{0.0pt}\pgfsys@lineto{4.83698pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{{}}{} {}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{1.1381pt}\pgfsys@invoke{ }{}\pgfsys@moveto{-0.56906pt}{4.26791pt}\pgfsys@lineto{4.83698pt}{4.26791pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}$, where $S^{(2)}_{A}=-\log(\mathrm{Tr}\rho^{2}_{A})$ is the second Rényi entanglement entropy and the subscript $A$ denotes the specific subsystem. Writing the Rényi entanglement entropy in terms of replica partition functions $S^{(2)}_{A}=-\log(Z^{(2)}_{A}/Z^{2})$ Calabrese and Cardy (2004), the Levin-Wen measurement can be expressed as $2\gamma=\log(Z^{(2)}_{\leavevmode\hbox to5.41pt{\vbox to5.98pt{\pgfpicture\makeatletter\hbox{\hskip 0.56905pt\lower-1.1381pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{} {{}{}}{} {}{} {{}{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{1.1381pt}\pgfsys@invoke{ }{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{4.26791pt}{0.0pt}\pgfsys@lineto{4.26791pt}{4.26791pt}\pgfsys@lineto{0.0pt}{4.26791pt}\pgfsys@lineto{0.0pt}{-0.56906pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}/Z^{(2)}_{\leavevmode\hbox to5.98pt{\vbox to5.41pt{\pgfpicture\makeatletter\hbox{\hskip 0.56905pt\lower-0.56905pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{} {{}{}}{} {}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{1.1381pt}\pgfsys@invoke{ }{}\pgfsys@moveto{4.83698pt}{4.26791pt}\pgfsys@lineto{0.0pt}{4.26791pt}\pgfsys@lineto{0.0pt}{0.0pt}\pgfsys@lineto{4.83698pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}})-\log(Z^{(2)}_{\leavevmode\hbox to5.98pt{\vbox to5.41pt{\pgfpicture\makeatletter\hbox{\hskip 1.1381pt\lower-0.56905pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{} {}{} {}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{1.1381pt}\pgfsys@invoke{ }{}\pgfsys@moveto{-0.56906pt}{0.0pt}\pgfsys@lineto{4.26791pt}{0.0pt}\pgfsys@lineto{4.26791pt}{4.26791pt}\pgfsys@lineto{-0.56906pt}{4.26791pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}/Z^{(2)}_{\leavevmode\hbox to6.54pt{\vbox to5.41pt{\pgfpicture\makeatletter\hbox{\hskip 1.1381pt\lower-0.56905pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{1.1381pt}\pgfsys@invoke{ }{}\pgfsys@moveto{-0.56906pt}{0.0pt}\pgfsys@lineto{4.83698pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{{}}{} {}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{1.1381pt}\pgfsys@invoke{ }{}\pgfsys@moveto{-0.56906pt}{4.26791pt}\pgfsys@lineto{4.83698pt}{4.26791pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}})$. To compute these partition function ratios numerically we have adapted a recently introduced algorithm D’Emidio (2019) to the current problem. The method introduces a one parameter family of partition functions $\mathcal{Z}_{AB}^{(2)}(\lambda)$ that interpolates between the two partition functions ($Z^{(2)}_{A}$ and $Z^{(2)}_{B}$) appearing in the ratio. In this extended ensemble the log ratio takes the form of a $\lambda$ integral of a simple Monte Carlo estimator sup . We have also tried other techniques to calculate the EE including the energy integration method Melko et al. (2010) used in Isakov et al. (2011), however we find the current method to be better suited to our problem.

At large values of $N$ and in the quantum spin liquid phase at the low temperatures of interest, it is difficult to efficiently sample our phase space using only traditional QMC loop updates. To improve the quality of our entanglement data we have incorporated annealing and replica exchange methods sup . With these improvements we are able to measure $\gamma$ reliably at moderately low temperatures, after which we encounter difficulties with equilibration and ergodicity. As we shall see, this allows us to observe the first plateau at $\log(2)/2$ but not the second plateau at $\log(2)$. Fig. 3 shows the TEE as a function of inverse temperature $\beta$ for the SO($N$) model with $N=8$ to $N=11$ on an $L=8$ lattice. As $T$ is lowered, we clearly see a pronounced signal in the TEE for $N\geq 10$ near a plateau at $\log(2)/2$. For $N\leq 9$ on the other hand $\gamma$ goes to zero in the low temperature regime, consistent with the study of the quadruplar order parameter show in Fig. 1. Interestingly, even though for $L=8$ the difference region in each ratio contains only $3\times 2\times 2$ we see reasonable quantization at the first plateau. To test how $\gamma$ scales as the system size $L$ (and the subsystem size) is scaled up we present the SO(11) TEE data for $L=8,12,16$ in Fig. 4. The data shows clear convergence to $\log(2)/2$ as $L$ is increased. We have also performed measurements at lower temperatures in an effort to see the second quantized plateau at $\log(2)$ and despite signals that are consistent with this picture, proper equilibration here remains challenging and will be saved for future studies.

In conclusion, we have unambiguously identified in sign-free Monte Carlo simulations, a ${\mathbb{Z}}_{2}$ quantum spin liquid in a simple model of magnetism on a Kagome lattice. Our work paves the way to study various interesting questions, including the theory of phase transitions out of the QSL, the role of isolated impurities, as well as the effect of large scale disorder in QSLs.

The computational work presented here was carried out using the XSEDE awards TG-DMR130040 and TG-DMR140061. Financial support was received through NSF DMR-1611161.

I Supplemental material

I.1 Details on entanglement entropy measurements

As described in the main text, we make use of the Levin-Wen construction Levin and Wen (2006) to extract the topological entanglement entropy, which can be written as

[TABLE]

Here $S^{(2)}_{A}=-\log(\mathrm{Tr}\rho^{2}_{A})$ is the second Rényi entanglement entropy and the subscript $A$ describes the shape of the region on an $L$ x $L$ lattice that is traced only once. We also assume that the region size is large compared to the lattice spacing and small compared to the total system size. Writing the Rényi entanglement entropy in terms of replica partition functions $S^{(2)}_{A}=-\log(Z^{(2)}_{A}/Z^{2})$ Calabrese and Cardy (2004), the Levin-Wen measurement can be expressed as

[TABLE]

In order to compute ratios such as these we employ the equilibrium version of a nonequilibrium method that was very recently introduced D’Emidio (2019). We have tried many different methods for these calculations, including energy integration as was used in Isakov et al. (2011) and described in Melko et al. (2010), but we find the present method to be more efficient for our model.

In [D’Emidio, 2019] it was shown that the $\log$ ratio of partition functions $\log(Z^{(2)}_{A}/Z^{(2)}_{B})$ can be computed by introducing a weighted sum over replica partition functions $\mathcal{Z}_{AB}^{(2)}(\lambda)$ that depends on an external field $\lambda$ such that $\mathcal{Z}_{AB}^{(2)}(0)=Z^{(2)}_{B}$ and $\mathcal{Z}_{AB}^{(2)}(1)=Z^{(2)}_{A}$. In other words $\lambda$ couples to the trace topology of the replica partition functions. If we assume that $B$ is a subset of $A$, we can write $\mathcal{Z}_{AB}^{(2)}(\lambda)$ explicitly as

[TABLE]

Here $C$ is summed over all all proper subsets of the set $A-B$, from the empty set $\o$ up to and including $A-B$ itself. Here $N_{A}$, $N_{B}$, and $N_{C}$ are the number of sites in the sets $A$, $B$, and $C$ respectively. Notice that when $\lambda=0$, only $C=\o$ survives and the sum equals $Z^{(2)}_{B}$. Also, when $\lambda=1$ only $C=A-B$ survives and the sum equals $Z^{(2)}_{A}$, as intended. The log ratio can then be computed as

[TABLE]

This can be efficiently measured in QMC as the equilibrium average $\partial\log\mathcal{Z}^{(2)}_{AB}(\lambda)/\partial\lambda=\frac{\langle N_{C}\rangle_{\lambda}}{\lambda(1-\lambda)}-\frac{N_{A}-N_{B}}{1-\lambda}$, where the average is taken in the $\mathcal{Z}^{(2)}_{AB}(\lambda)$ ensemble at a fixed value of $\lambda$.

In Fig. 5 we show the two different types of $\mathcal{Z}^{(2)}_{AB}(\lambda)$ geometries that are used to compute Eq. (8) for an $L=8$ system. In these pictures black sites are traced once, white sites are traced twice, and grey sites can fluctuate independently between a single and double trace according to the value of $\lambda$. When $\lambda=0(1)$ all of the spins in the grey region are traced twice (once), respectively. The QMC simply measures $\langle N_{C}\rangle_{\lambda}$, or the average number of spins in the grey region that are traced once for a fixed value of $\lambda$ in equilibrium.

In this framework the Levin-Wen measurement takes on a simple and intuitive form, we can write Eq. 8 as

[TABLE]

where $\Delta\langle N_{C}\rangle_{\lambda}$ is the difference in the average number of single trace spins in the grey region for the two different Levin-Wen geometries in Fig. 5.

In order to measure $\gamma$ as a function of temperature, we perform thermal annealing from randomly initialized configurations at high temperature and slowly cool down to low temperature. This is done for 24 different values of $\lambda$ for both $\mathcal{Z}^{(2)}_{\leavevmode\hbox to5.41pt{\vbox to5.98pt{\pgfpicture\makeatletter\hbox{\hskip 0.56905pt\lower-1.1381pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{} {{}{}}{} {}{} {{}{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{1.1381pt}\pgfsys@invoke{ }{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{4.26791pt}{0.0pt}\pgfsys@lineto{4.26791pt}{4.26791pt}\pgfsys@lineto{0.0pt}{4.26791pt}\pgfsys@lineto{0.0pt}{-0.56906pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\,\leavevmode\hbox to5.98pt{\vbox to5.41pt{\pgfpicture\makeatletter\hbox{\hskip 0.56905pt\lower-0.56905pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{} {{}{}}{} {}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{1.1381pt}\pgfsys@invoke{ }{}\pgfsys@moveto{4.83698pt}{4.26791pt}\pgfsys@lineto{0.0pt}{4.26791pt}\pgfsys@lineto{0.0pt}{0.0pt}\pgfsys@lineto{4.83698pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}(\lambda)$ and $\mathcal{Z}^{(2)}_{\leavevmode\hbox to5.98pt{\vbox to5.41pt{\pgfpicture\makeatletter\hbox{\hskip 1.1381pt\lower-0.56905pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{} {}{} {}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{1.1381pt}\pgfsys@invoke{ }{}\pgfsys@moveto{-0.56906pt}{0.0pt}\pgfsys@lineto{4.26791pt}{0.0pt}\pgfsys@lineto{4.26791pt}{4.26791pt}\pgfsys@lineto{-0.56906pt}{4.26791pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\,\leavevmode\hbox to6.54pt{\vbox to5.41pt{\pgfpicture\makeatletter\hbox{\hskip 1.1381pt\lower-0.56905pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}{}}{{}}{} {{}{}}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{1.1381pt}\pgfsys@invoke{ }{}\pgfsys@moveto{-0.56906pt}{0.0pt}\pgfsys@lineto{4.83698pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{{}}{} {}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{1.1381pt}\pgfsys@invoke{ }{}\pgfsys@moveto{-0.56906pt}{4.26791pt}\pgfsys@lineto{4.83698pt}{4.26791pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}(\lambda)$. On regular intervals the thermal annealing schedule is paused and measurements of $\langle N_{C}\rangle_{\lambda}$ are made for each ratio. To improve measurement statistics we employ replica exchange (within each Levin-Wen geometry separately) as a function of $\lambda$, where neighboring configurations can swap $\lambda$ values probabilistically. Since traditional QMC loop updates are extremely inefficient deep in the spin liquid phase, we compute our QMC averages over independent thermal annealing realizations.

I.2 QMC versus exact diagonalization

Here we compare our QMC method against exact results on a 2 x 2 lattice for $N=2$. We have tried to make the comparison as close as possible to the measurement required to compute the topological entanglement entropy. We have therefore chosen the Levin-Wen geometries as depicted in Fig. 6. These are used to compute $\Delta\langle N_{C}\rangle_{\lambda}/\lambda(1-\lambda)$ shown in the main plot, where we see perfect agreement between the QMC and ED.

I.3 TEE Equilibration

Our thermal annealing schedule consists of 2000 equilibration sweeps at each value of $\beta$ on a fine grid of points that are equally spaced on a $\log$ scale (a geometric progression with 1400 points between $\beta=0.1$ and $\beta=100$). The data presented in this work is generated by pausing the thermal annealing schedule at regular intervals to make measurements. At each measured value of $\beta$ we perform 1000 equilibration sweeps, then 20000 measurement sweeps followed by another 1000 equilibration sweeps and another 20000 measurement sweeps. We can separately average measurements from the first and second segments, where we expect statistical agreement in the case of properly equilibrated configurations. This is shown in Fig. 7, where agreement is found for the $L=16$ SO(11) system. The data presented in the main text is the average over both measurement segments.

I.4 TEE density

Finally it is interesting to look at the behavior of $d\gamma/d\lambda=\Delta\langle N_{C}\rangle_{\lambda}/2\lambda(1-\lambda)$ that, when integrated from $\lambda=0$ to $\lambda=1$, gives the Levin-Wen topological entanglement entropy values presented in the main text. This is shown in Fig. 8 for SO(11) $L=8,12,16$ at the largest values of $\beta$ presented in the main text. We can see that the largest contribution to $\gamma$ comes near $\lambda=0.5$ and tapers off to zero near the extremes for sufficiently large system sizes where the quantized value $\log(2)/2$ is observed.