Abelian and non-Abelian chiral spin liquids in a compact tensor network representation

TL;DR

This paper introduces tensor network-based variational ansatze for Abelian and non-Abelian chiral spin liquids in the Kitaev model on the star lattice, revealing their properties and anyonic excitations with high accuracy.

Contribution

It develops compact tensor network representations of chiral spin liquids, accurately captures phase diagram features, and characterizes anyons and edge modes in the non-Abelian phase.

Findings

Variational energies are highly accurate, with LG states exact at certain limits.

Short-range LG effectively describes Abelian CSLs.

Non-Abelian CSLs connect adiabatically to critical LG with macroscopic loops.

Abstract

We provide new insights into the Abelian and non-Abelian chiral Kitaev spin liquids on the star lattice using the recently proposed loop gas (LG) and string gas (SG) states [H.-Y. Lee, R. Kaneko, T. Okubo, N. Kawashima, Phys. Rev. Lett. 123, 087203 (2019)]. Those are compactly represented in the language of tensor network. By optimizing only one or two variational parameters, accurate ansatze are found in the whole phase diagram of the Kitaev model on the star lattice. In particular, the variational energy of the LG state becomes exact(within machine precision) at two limits in the model, and the criticality at one of those is analytically derived from the LG feature. It reveals that the Abelian CSLs are well demonstrated by the short-ranged LG while the non-Abelian CSLs are adiabatically connected to the critical LG where the macroscopic loops appear. Furthermore, by constructing the minimally entangled states and exploiting their entanglement spectrum and entropy, we identify the nature of anyons and the chiral edge modes in the non-Abelian phase with the Ising conformal field theory.

Equations

$$\hat{\mathcal{H}}=-\frac{J}{4}\sum _{⟨ij⟩\in \gamma }{\hat{\sigma }}_i^{\gamma }{\hat{\sigma }}_j^{\gamma }-\frac{J^{\prime }}{4}\sum _{⟨ij⟩\in {\gamma }^{\prime }}{\hat{\sigma }}_i^{{\gamma }^{\prime }}{\hat{\sigma }}_j^{{\gamma }^{\prime }},$$

$$⟨\theta ,\gamma |{\hat{\sigma }}^{{\gamma }^{\prime }}|\theta ,\gamma ⟩={\delta }_{\gamma {\gamma }^{\prime }}\cos \theta +(1-{\delta }_{\gamma {\gamma }^{\prime }})\frac{\sin \theta }{\sqrt{2}},$$

$$\text{includegraphics}.$$

$$⟨{\psi }_{\mathrm{LG}}(\theta )|{\psi }_{\mathrm{LG}}(\theta )⟩=c{Z}_{O(1)}\left(\frac{q^2}{1-r+r^2}\right),$$

$$\text{includegraphics},$$

$$\text{includegraphics}.$$

$$Q_{ijk}^{s{s}^{\prime }}={\tau }_{ijk}[({\hat{\sigma }}^x{)}^{1-i}({\hat{\sigma }}^y{)}^{1-j}({\hat{\sigma }}^z{)}^{1-k}{]}_{s{s}^{\prime }},{\tau }_{ijk}=\{\begin{array}{ll}-i& \mathrm{if}i+j+k=0\\ 1& \mathrm{if}i+j+k=2\end{array},$$

$$⟨{\psi }_0(\theta )|{\psi }_0(\theta )⟩=⟨\Psi (\theta )|{\hat{Q}}_{\mathrm{LG}}^{†}{\hat{Q}}_{\mathrm{LG}}|\Psi (\theta )⟩=N_{\Gamma }⟨\Psi (\theta )|{\hat{Q}}_{\mathrm{LG}}|\Psi (\theta )⟩=N_{\Gamma }\sum _{G\in \Gamma }{x}_G(\theta ),$$

$$\text{includegraphics},$$

$$⟨{\psi }_0(\theta )|{\psi }_0(\theta )⟩=N_{\Gamma }\sum _{G^{\prime }\in {\Gamma }^{\prime }}(1+r^3{)}^{n-n_{G^{\prime }}}[q^2(1+r){]}^{n_{G^{\prime }}}=N_{\Gamma }(1+r^3{)}^n\sum _{G^{\prime }\in {\Gamma }^{\prime }}{\left(\frac{q^2(1+r)}{1+r^3}\right)}^{n_{G^{\prime }}}.$$

$$R_{ijk}^{s{s}^{\prime }}={\zeta }_{ijk}{\left[({\hat{\sigma }}^x{)}^i({\hat{\sigma }}^y{)}^j({\hat{\sigma }}^z{)}^k\right]}_{s{s}^{\prime }},{\zeta }_{ijk}=\{\begin{array}{ll}1& \mathrm{if}i+j+k=0\\ c_1/c_2& \mathrm{if}i+j+k=1\end{array},$$

$${\zeta }_{ijk}^{\mathrm{site}8}=\{\begin{array}{ll}1& \mathrm{if}i+j+k=0\\ c_1& \mathrm{if}j=0\mathrm{and}k+i=1\\ c_2& \mathrm{if}j=1\mathrm{and}k+i=0\end{array}.$$

$${\zeta }_{ijk}=\{\begin{array}{ll}\cos \beta & \mathrm{if}i+j+k=0\\ \sin \beta \sqrt{\cos \alpha }/\sin \beta \sqrt{\sin \alpha }& \mathrm{if}i+j+k=1\end{array}.$$

$$v=\left(\begin{array}{ccc}0& & i\\ 1& & 0\end{array}\right),$$

$${\sigma }_{s{s}^{\prime \prime }}^x{Q}_{ijk}^{s^{\prime \prime }{s}^{\prime }}=v_{j{j}^{\prime }}{v}_{k{k}^{\prime }}^{\ast }{Q}_{i{j}^{\prime }{k}^{\prime }}^{s{s}^{\prime }},{\sigma }_{s{s}^{\prime \prime }}^y{Q}_{ijk}^{s^{\prime \prime }{s}^{\prime }}=v_{k{k}^{\prime }}{v}_{i{i}^{\prime }}^{\ast }{Q}_{i^{\prime }j{k}^{\prime }}^{s{s}^{\prime }},{\sigma }_{s{s}^{\prime \prime }}^z{Q}_{ijk}^{s^{\prime \prime }{s}^{\prime }}=v_{i{i}^{\prime }}{v}_{j{j}^{\prime }}^{\ast }{Q}_{i{j}^{\prime }{k}^{\prime }}^{s{s}^{\prime }}.$$

$$\text{includegraphics},$$

Full text

I introduction

Discovery of the fractional quantum Hall (FQH) effectTsui et al. (1982) had brought a paradigm shift in understanding of condensed phases of matter. Exotic quantum liquid states, i.e., chiral spin liquids (CSL), were proposed as the ground states of the FQH systemKalmeyer and Laughlin (1987); Wen et al. (1989), which cannot be featured by Landau’s symmetry breaking theory but the so-called topological order. The topological order can be interpreted as the pattern of long-range entanglement which leads to the ground state degeneracy depending only on the topology of systemWen (1990). Those CSL ansatze successfully explained the nature of the FQH fluids such as the fractional statistics of quasiparticles (or anyons)Arovas et al. (1984); Halperin (1984). Furthermore, the anyons obeying the non-Abelian braiding statistics were theoretically realized in the FQH systemWen (1991a); Moore and Read (1991); Greiter and Thomale (2009). Due to the robust topological degeneracy against the local perturbations and exotic statistics of anyons, the non-Abelian topological states have been proposed as a promising platform for fault-tolerant quantum computingKitaev (2003) and thus attracted lots of attention in the field of quantum information for the last decadeNayak et al. (2008). Another interesting feature of the FQH fluids and CSLs is that the chiral gapless edge modes appear at the boundary of the system, and it leads to perfect heat conduction at the edgeWen (1991b). The edge states are described by the conformal field theories (CFT) which also characterize and hence have been employed to identify the topological orderLi and Haldane (2008); Qi et al. (2012); Zhu et al. (2015); Yang et al. (2015); Poilblanc et al. (2015, 2016); Poilblanc (2017); Chen et al. (2018). By solving the Kitaev modelKitaev (2006) on the star lattice (KSM), Yao and Kivelson showed the existence of the CSL as an exact ground state of local Hamiltonian and found Abelian and non-Abelian phases characterized by the topological degeneracies four and three on the torus, respectivelyYao and Kivelson (2007); Chung et al. (2010).

The aim of this Letter is to understand the Abelian and non-Abelian Kitaev CSLs without referring to the Majorana fermion. Recently, a particular LG state and its extension, which is referred to as SG state, have been proposed as ansatze for the Kitaev spin liquid (KSL) on the honeycomb latticeKitaev (2006) in a compact tensor product state (TPS) representationLee et al. (2019). The LG ansatz was found to reflect most qualitative features of the KSL, and SG provides a quantitatively accurate approximation to the KSL while keeping the qualitative features intact. In what follows, we reinterpret the CSLs as the LG and SG states and provide direct evidences identifying the topological order in each phase.

II model

The KSM is defined asYao and Kivelson (2007)

[TABLE]

where $\hat{\sigma}_{i}^{\gamma}$ stands for the Pauli matrix with $\gamma,\gamma^{\prime}=x,y,z$, while $\langle ij\rangle_{\gamma}$ and $\langle ij\rangle_{\gamma^{\prime}}$ denote the nearest-neighbor pair respectively on the intra-triangle ($\gamma$) and inter-triangle ($\gamma^{\prime}$) bonds connecting sites $i$ and $j$ as defined in Fig. 1 (a). Note that the Hamiltonian commutes with two types of flux operators defined on the triangle plaquette $\hat{V}_{p}=\hat{\sigma}_{1}^{x}\hat{\sigma}_{2}^{y}\hat{\sigma}_{3}^{z}$ and dodecagon plaquette $\hat{W}_{p}=\hat{\sigma}_{1}^{x}\hat{\sigma}_{2}^{z}\hat{\sigma}_{3}^{y}\hat{\sigma}_{4}^{x}\cdots\hat{\sigma}_{12}^{y}$Yao and Kivelson (2007), where the site indices are defined in Fig. 1 (a). Therefore, the Hamiltonian is block-diagonalized, and each block is characterized by the set of the flux numbers or eigenvalues of flux operators $\{V_{p}=\pm 1,W_{p}=\pm 1\}$. Since the operator $\hat{V}_{p}$ consists of three Pauli matrices, the time-reversal transformation flips its flux number, i.e., $\mathcal{T}\hat{V}_{p}\mathcal{T}^{-1}=-\hat{V}_{p}$, and hence the TR-symmetry is spontaneously broken in eigenstates of $\hat{\mathcal{H}}$. It was foundYao and Kivelson (2007) that the ground states do not break any lattice symmetry (that is, the CSLs) and live in the vortex-free sector, i.e., $\{\hat{W}_{p}=1,\hat{V}_{p}=1\}$. In addition, the model exhibits a topological phase transition between the non-Abelian and Abelian CSLs at $J^{\prime}/J=\sqrt{3}$Yao and Kivelson (2007).

III variational ansatze

III.1 Loop gas states

We begin with defining a product state $|\Psi(\theta)\rangle=\otimes_{\alpha}|\theta,\gamma_{\alpha}\rangle_{\alpha}$ with a local magnetic state given by

[TABLE]

and $\gamma_{\alpha}$ being the inter-triangle bond at site $\alpha$. For instance, at site-8 in Fig. 1 (a), we assign the state $|\theta,y\rangle$. Applying the loop gas (LG) operator $\hat{Q}_{\rm LG}$ defined in Ref. Lee et al. (2019) on top of $|\Psi(\theta)\rangle$ results in a LG state $|\psi_{\rm LG}(\theta)\rangle\equiv\hat{Q}_{\rm LG}|\Psi(\theta)\rangle$, which is simply a superposition of all possible loop configurations of magnetic states, i.e.,

[TABLE]

Here, the empty site denotes the state $|\theta,\gamma_{\alpha}\rangle$ while the loop-occupied site stands for $\hat{\sigma}^{\gamma}|\theta,\gamma_{\alpha}\rangle$ depending on the direction of loop as depicted in Fig. 1 (b). Due to the LG operator, the desired symmetries, $Z_{2}$ gauge redundancy and the vortex-freeness are guaranteed in the LG stateLee et al. (2019). Notice that the norm of $|\psi_{\rm LG}(\theta)\rangle$ maps into the partition function of the classical $O(1)$ LG modelNienhuis (1982) with the local weights of loops being $r=\cos\theta$ and $q=\sin\theta/\sqrt{2}$ along the triangle and dodecagon plaquttes, respectively, on the star latticeLee et al. (2019). After simple algebra, the partition function can even be mapped into that on the honeycomb lattice, i.e., $Z_{O(1)}(x)$ with $x$ being the fugacity per site:

[TABLE]

where $c$ is a constant(see Appendix A). Now, one can optimize the variational parameter $\theta$ to minimize the energy for a given $(J,J^{\prime})$. For simplicity, let us parameterize the exchange couplings as $J^{\prime}/J=\tan\phi$. However, at $\phi=0$, we consider $J\rightarrow\infty$ while $J^{\prime}=1$ being finite to avoid the trivial solution at $J^{\prime}=0$, and vice versa at $\phi=\pi/2$. We employ the corner transfer matrix renormalization group (CTMRG) methodNishino and Okunishi (1996); Orús and Vidal (2009); Corboz et al. (2010) to measure the energy $E=\langle\psi_{\rm LG}|\hat{\mathcal{H}}|\psi_{\rm LG}\rangle$ and find $\theta^{*}(\phi)$ minimizing the energy at a given $\phi$. The resulting $E$ and $\theta^{*}$ are presented in Fig. 2 (a) and (b), respectively, as a function of $\phi$. Here, we present the exact energy $E_{\rm ex}$. The optimal local weights $(q^{*},r^{*})=(\sin\theta^{*}/\sqrt{2},\cos\theta^{*})$ are also presented in Fig. 2 (b), which provide new insights into the nature of each phase. As one can see, the variational energy of the LG ansatz is quite accurate in $0.4\pi<\phi\leq 0.5\pi$, where the energy deviation $dE=1-E/E_{\rm exact}$ is less than $0.1\%$ as shown in the inset of Fig. 2 (a). In particular, the energy becomes exact (within machine precision) at $\phi=0.5\pi$ at which $r$ is maximized while $q$ vanishes as shown in Fig. 2 (b). It indicates that the configurations with only triangle loops and holes survive. In this sense, the Abelian CSL phase can be understood as the triangle loop gas in which longer loops are suppressed.

Interestingly, the variational energy becomes also exact (within machine precision) as approaching the opposite limit $\phi=0$ where the excitation gap closes as $J^{\prime 2}/J$Yao and Kivelson (2007); Dusuel et al. (2008). Notice that $Z_{O(1)}(x)$ in Eq. (4) becomes critical at $x_{c}=1/\sqrt{3}$Nienhuis (1982), which leads to the critical LG state at $(q_{c},r_{c})=(\sqrt{2\sqrt{3}-3},2-\sqrt{3})$. Remarkably, the LG state is optimized with $(q_{c},r_{c})$ at $\phi=0$ which implies that the ground state is the critical LG state exhibiting macroscopic loops. Furthermore, its low-energy physics is described by the Ising CFTNienhuis (1982) which is consistent with the expected oneYao and Kivelson (2007); Månsson et al. (2013); Lahtinen et al. (2014); Meichanetzidis et al. (2016). Therefore, according to these circumstantial evidences, we may conclude that the LG ansatz at $\phi=0$ is the exact ground state. It tells us that the non-Abelian CSL around $\phi=0$ are well described by the long-ranged LG states which is qualitatively distinct from the short-ranged feature of the Abelian CSL. Excluing $\phi=0$, all LG ansatze map into the gapped phase of $Z_{O(1)}(x)$, i.e., $x<1/\sqrt{3}$, ensuring its gapped nature.

III.2 String gas states

In the non-Abelian phase and around the phase transition point, the variational energies of the LG ansatz are away from the exact ones. Therefore, in order to lower the energy and find better ansatze, we apply the dimer gas (DG) operator $\hat{R}_{\rm DG}$ as suggested in Ref. Lee et al. (2019),

[TABLE]

where the dimer on a bond $\langle ij\rangle$ stands for $\hat{\sigma}^{\gamma}_{i}\hat{\sigma}^{\gamma}_{j}$ depending on the bond, and different colors for the intra-bond and inter-bond. Note that one can introduce variational parameters in the DG operator, which determine the fugacities of the dimersLee et al. (2019), and optimize them to lower energy. Here, we introduce two independent fugacites such that the one on the intra-triangle bonds ($c_{1}$) and another one on the inter-triangle bond ($c_{2}$), i.e., $\hat{R}_{\rm DG}=\hat{R}_{\rm DG}(c_{1},c_{2})$. This operator does not spoil the symmetries and gauge structure of the ansatzeLee et al. (2019). Then, we employ the state $|\psi_{\rm SG}(c_{1},c_{2})\rangle=\hat{Q}_{\rm LG}\hat{R}_{\rm DG}(c_{1},c_{2})|\Psi_{(111)}\rangle$ as our ansatz which can be regarded as the SG stateLee et al. (2019). Here, the state $|\Psi_{(111)}\rangle\equiv|\Psi(\tan^{-1}\sqrt{2})\rangle$ is the the (111)-magnetic state where all spins point to $(111)$-direction. In fact, one can apply the DG operator on the general product state $|\Psi(\theta)\rangle$ to have an additional variational parameter $\theta$. Instead, for brevity, we fix the initial product state as $|\Psi_{(111)}\rangle$ and optimize $c_{1}$ and $c_{2}$ for a given $\phi$. The obtained variational energy is presented in Fig. 2 (a), of which the inset is the energy deviation from the exact one (see Appendix B for details on the optimized parameters). As one can see, the DG operator drastically reduces the energy and provides reasonably good ansatz even around the transition point. Furthermore, the second derivative of the energy allows us to estimate the transition point correctly(see Appendix B). Note that, in the non-Abelian phase, the variational energies are particularly good around $\phi=0.25\pi$. This is because $|\Psi_{(111)}\rangle$ is used as the initial state, which is optimal at $\phi=0.25\pi$. We thus believe that one can obtain even better ansatze ($dE\sim O(10^{-4})$) throughout the non-Abelian phase by choosing the initial state $|\Psi(\theta)\rangle$ properly.

IV minimally entangled states

So far we have considered the ansatz only on the infinite system. Now we discuss the ansatz on the compact manifold (e.g., torus), where the topological sectors allows us to distinguish the Abelian and non-Abelian phases. With periodic boundary conditions (PBC), one should also consider the so-called global flux measured by the flux operator $\hat{\Phi}_{\Gamma}=\prod_{i\in\Gamma}\hat{\sigma}_{i}^{\gamma_{i}}$ defined on a non-contractible closed path $\Gamma$Kitaev (2006). Its eigenvalues $\pm 1$ determine the topological sector, say even (odd) sector for $+1$ ($-1$). However, it turns out that our ansatze with PBC are not eigenstates of the global flux operators. To be more specific, let us consider ansatz on the torus (or an infinitely long cylinder) and then apply $\hat{\Phi}_{y}$ wrapping the inner tube of torus (say $y$-direction). One can verify that multiplying $\hat{\Phi}_{y}$ is equivalent to the gauge twisting along a closed path encircling the tube as illustrated below:

[TABLE]

where the five-leg tensor is composed of six onsite tensors in the unit-cell, and

[TABLE]

Here, the red squares denote $\hat{\Phi}_{y}$ whereas the yellow one stands for the non-trivial element of the $Z_{2}$ invariant gauge group of our ansatz, $g=\hat{\sigma}^{z}$Lee et al. (2019). The detailed derivation is presented in Appendix C. Note that RHS of Eq. (5) is equivalent to the procedure of creating a vortex pair, moving one of those around the loop $\Gamma_{y}$, and then annihilating the pair. We define $G_{y}=\otimes_{i=1}^{L_{y}}g$ as the string of $g$ wrapping the inner tube where $L_{y}$ is the circumference in units of the unit-cell, i.e., the ring of yellow tensors in the right-hand side of Eq. (5). Therefore, acting $\hat{\Phi}_{y}$ changes our ansatz $|\psi\rangle$ (regardless of LG or SG) to a different state $|\psi_{y}\rangle$ that is the $G$-inserted $|\psi\rangle$. However, since the square of $\hat{\Phi}_{y}$ is identity, its eigenstates are simply obtained by $|\psi_{\pm}\rangle=|\psi\rangle\pm|\psi_{y}\rangle$, and the subscript $\pm$ labels the global flux number, i.e., $\hat{\Phi}_{y}|\psi_{\pm}\rangle=\pm|\psi_{\pm}\rangle$. In a similar way, one can set the simultaneous eigenstates of both $\hat{\Phi}_{x}$ and $\hat{\Phi}_{y}$, i.e., $|\psi_{(\pm,\pm)}\rangle$ living in one of four topological sectors specified by $(\hat{\Phi}_{x},\hat{\Phi}_{y})=(\pm 1,\pm 1)$. Interestingly, in case of the LG ansatze, those in distinct sector are characterized by the parity of the number of non-contractible loop configurations in each directionKitaev (2003); Poilblanc et al. (2012). It can be easily seen that the action of $G_{y}$ gives a minus sign to all the configurations with odd number of non-contractible loops enclosing the hole of torus (say $x$-direction). Now, using those topologically degenerate ansatze, we construct the so-called minimally entangled states (MES) Zhang et al. (2012), e.g., $|\mathbb{I}\rangle=|\psi_{(+,+)}\rangle+|\psi_{(-,+)}\rangle$ and $|\mathfrak{m}\rangle=|\psi_{(+,-)}\rangle+|\psi_{(-,-)}\rangle$ characterized by each anyon ($\mathbb{I}$: trivial, $\mathfrak{m}$: vortex) flux threading the inner tube of torus. In this basis, one can read off the quantum dimension ($d_{i}$) of each anyon from the topological entanglement entropy (TEE), $\gamma_{i}=\log(D/d_{i})$, where $i=\mathbb{I},\mathfrak{m},\cdots$. The total quantum dimension $D=\sqrt{\sum_{i}d_{i}^{2}}$Kitaev and Preskill (2006); Levin and Wen (2006); Dong et al. (2008) is known to be four for the Abelian and non-Abelian KSL phasesYao and Qi (2010); Chung et al. (2010). To this end, we employ the bulk-edge correspondence in TPS Cirac et al. (2011) to evaluate the entanglement spectrum (ES) and entanglement entropy (EE) on the infinitely long cylinder. Here, we impose PBC in the $y$-direction.

Firstly, the results of TEE obtained from the LG and SG ansatze are presented in Fig. 3. Here, (a) and (b) show the EEs in each sector obtained from $|\psi_{\rm LG}\rangle$ at $\phi=0.02\pi$ and $\phi=0.4\pi$, respectively, as a function of the circumference $L_{y}$. As expected from the geometry of TPSVerstraete et al. (2006), all EEs follow the area lawSrednicki (1993): $S=\alpha L_{y}-\gamma_{i}$ where $\alpha$ is a non-universal prefactor, and $\gamma_{i}$ is extracted by fitting the data with linear functions (black solid lines). At $\phi=0.02\pi$, we obtained $(\alpha,\gamma_{i})=(0.5502,0.6786)$ and $(0.5535,0.3544)$ in each $\mathbb{I}$ and $\mathfrak{m}$ sector, respectively [Fig. 3 (a)]. Those TEE are remarkably close to the TEE in the vacuum sector (i.e., $\log 2$) and the $\sigma$-anyon (vortex) sector (i.e., $\log\sqrt{2}$) of Ising anyon modelKitaev (2006); Nayak et al. (2008). On the other hand, at $\phi=0.4\pi$ [Fig. 3 (a)], both EEs almost perfectly fit to $(\alpha,\gamma_{i})=(\log 2,\log 2)$, which is consistent with the one from the toric codeKitaev (2003, 2006). Similarly, we have extracted $\gamma_{i}$ at each $\phi$, and the results obtained from $|\psi_{\rm LG}\rangle$ and $|\psi_{\rm SG}\rangle$ are shown in Fig. 3 (c) and (d), respectively. Those of $|\psi_{\rm LG}\rangle$ are in an excellent agreement with the ones of Ising anyon model around $\phi=0$ and with the ones of the toric code mostly in the Abelian phase [Fig. 3 (c)]. Meanwhile, the SG ansatz gives almost consistent TEEs even in the non-Abelian phase agreeing with the ones of Ising anyon model and predicts the transition point rather correctly [Fig. 3 (d)].

Furthermore, the identification of the topological excitations becomes even clearer from characteristic structures in the ESLi and Haldane (2008). Figure 4 (a) and (b) present the ESs of $|\mathbb{I}\rangle$ and $|\mathfrak{m}\rangle$ obtained by the SG ansatz at $\phi=0.25\pi$. Here, the circumference is $L_{y}=6$, and the horizontal axis $k_{y}$ denotes the momentum. There are four branches of two distinct chiral modes in the $\mathbb{I}$-sector, which linearly disperse in one direction. Those are highlighted by the red and blue solid lines in Fig. 4 (a). Assuming the close ESs (dashed boxes) as degenerate levels, the degeneracy pattern is consistent with the ones of the primary fields ${\bf 1}$ (blue) and $\psi$ (red) and their descendants in the Ising CFTFriedan et al. (1984); Henkel (1999), respectively. On the other hand, in the $\mathfrak{m}$-sector, we find six branches of a single chiral mode of which the degeneracy counting obeys $\{1,1,1,2,2,3,4,5,7,\cdots\}$Henkel (1999), i.e., the characteristic of the primary field $\sigma$ and its descendants in the Ising CFT. In addition, the level spacings in the low-lying spectrum are in excellent agreement with the exact ones (see Appendix D). From our MES setup, the state $|\mathfrak{m}\rangle$ is expected to accommodate the vortex at each boundary, and thus the vortex is identified with the $\sigma$-anyon exhibiting non-Abelian braiding statisticsHenkel (1999); Nayak et al. (2008).

V conclusion

In this Letter, we show that the Abelian and non-Abelian CSL ground states of the KSM are well represented by the LG and SG states. In particualr, at both limits $\phi=0$ and $\pi/2$, the LG states become exact. Further, the gap closing at $\phi=0$ is understood by mapping the norm of ansatz into the partition function of the critical LG model. In addition, the fate of long-ranged loops is found to determine the Abelianess and non-Abelianess of CSL. By constructing the MES and measuring its TEE, we directly show that our ansatze host indeed the non-Abelian vortex with the quantum dimension $d_{\mathfrak{m}}=\sqrt{2}$. On the other hand, it becomes trivial, i.e., $d_{\mathfrak{m}}=1$, as the ansatz enters into the Abelian phase. We also identify the chiral edge modes in the non-Abelian phase with the Ising CFT, not the SU(2)2 Wess-Zumino-Witten theory conjectured in Ref. Yao and Kivelson (2007) by exploiting the level spacing and their degeneracy patternsFrancesco et al. (2012); Chen et al. (2018). We believe that the LG ansatze are the simplest CSLs in a compact representation, and therefore it could provide a platform bridging the quantum loop modelsFendley (2008) with the Abelian and non-Abelian topological states. It is also worth noting that our ansatze are the rigorous example explicitly revealing that general TPSs can represent the chiral gapped states. In the case of the fermionic Gaussian TPS, there exists a no-go theoremDubail and Read (2015) prohibiting the chiral Gaussian TPS to be gapped. However, it was not so clear whether the theorem applies to generic TPSs or notPoilblanc et al. (2015); Yang et al. (2015); Poilblanc et al. (2016); Poilblanc (2017); Chen et al. (2018). We believe that our ansatze are the counter evidence against the generalization of the theorem. In technical aspects, two independent optimization schemes, which can be combined together, are introduced for the LG and SG ansatz. Those can be employed to study the anisotropic Kitaev model on the honeycomb lattice and its extensions which are relevant to Kitaev materials such as $\alpha$-RuCl3Jackeli and Khaliullin (2009); Chaloupka et al. (2010); Singh and Gegenwart (2010); Kimchi and You (2011); Gohlke et al. (2018); Banerjee et al. (2016, 2018). Also, the variational ansatze are of interest of the deformed topological wavefunctionXu and Zhang (2019).

Acknowledgements.

The computation in the present work was executed on computers at the Supercomputer Center, ISSP, University of Tokyo, and also on K-computer (project-ID: hp190196). N.K.’s work is funded by MEXT KAKENHI No.19H01809. H.-Y.L. was supported by MEXT as “Exploratory Challenge on Post-K computer” (Frontiers of Basic Science: Challenging the Limits). T.O. was supported by JSPS KAKENHI No.15K17701 and 19K03740. R.K. was supported by MEXT as ”Priority Issue on Post-K computer” (Creation of New Functional Devices and High-Performance Materials to Support Next-Generation Industries).

Appendix A Norm of the zeroth ansatz

In this section, we show that the norm of the loop gas (LG) ansatz $|\psi_{0}(\theta)\rangle=\hat{Q}_{\rm LG}|\Psi(\theta)\rangle$ can be exactly mapped into the partition function of the $O(1)$ loop gas model on the honeycomb lattice. To this end, we first note that the LG operator is an hermitian projector: $\hat{Q}_{\rm LG}^{\dagger}=\hat{Q}_{\rm LG}$ and $(\hat{Q}_{\rm LG})^{2}=N_{\Gamma}\hat{Q}_{\rm LG}$ where $N_{\Gamma}$ is the total number of loop configurations on the star lattice. The LG operator is efficiently represented by the tensor product operatorLee et al. (2019), i.e., $\hat{Q}_{\rm LG}={\rm tTr}\prod_{\alpha}Q_{i_{\alpha}j_{\alpha}k_{\alpha}}^{ss^{\prime}}|s\rangle\langle s^{\prime}|$ where tTr stands for the tensor trace, $\alpha$ labels the site index, the building block tensor

[TABLE]

and the indices of physical and virtual legs are $s,s^{\prime}=0,1$ and $i,j,k=0,1$, respectively. Also, the virtual legs $i,j$ and $k$ lie on the $x,y$ and $z$ bonds in the model (see Fig. 1 (a) in the main text), respectively. Since the LG operator is obtained by summing over all possible loop operators (that is, product of $\hat{\sigma}^{x}$, $\hat{\sigma}^{y}$ and $\hat{\sigma}^{z}$ along the loops), it is straightforward to verify $(\hat{Q}_{\rm LG})^{2}=N_{\Gamma}\hat{Q}_{\rm LG}$ using the manipulation rules of loop defined in Ref. Lee et al. (2019). One can also easily show its hermiticity using the $Q$-tensor in Eq. (6). Now, let us compute the norm of the LG ansatz:

[TABLE]

where we used the hermiticity and idempotence (up to overall $N_{\Gamma}$), and $G$ denotes each loop configuration, $x_{G}(\theta)=\langle\Psi(\theta)|\hat{Q}_{G}|\Psi(\theta)\rangle$ is the weight of a loop operator $\hat{Q}_{G}$. With the definition of $|\Psi(\theta)\rangle$ and Eq. (2) in the main text, one can easily verity $x_{G}(\theta)=(\sin\theta)^{l_{G}^{t}}(\cos\theta/\sqrt{2})^{l_{G}^{d}}$ where $l_{G}^{t}$ and $l_{G}^{d}$ are the total lengths of partial loops along the triangle and dodecagon plaquettes, respectively, in the configuration $G$. Now, let us consider local loop configurations on a single triangle plaquetteBatchelor (1998). There are eight configurations on the single triangle plaquette which are depicted in the left hand side of graphical equations below:

[TABLE]

where $r=\cos\theta$ and $q=\sin\theta/\sqrt{2}$ as defined in the main text. As shown equation above, summing two of them can be regarded as a configuration on the honeycomb lattice by treating the triangle plaquette as a single site. Then, our task falls into the calculation of the partition function of the $O(1)$ LG model on the honeycomb lattice, where the weights of hole and loop per site are $1+r^{3}$ and $q^{2}(1+r)$, respectively, and it is simply given by

[TABLE]

Here, $\Gamma^{\prime}$ denotes all possible loop configurations on the honeycomb lattice, and $n$ is the total number of sites on the honeycomb lattice while $n_{G^{\prime}}$ stands for the total length of loops in a configuration $G^{\prime}$. The RHS is identical to the partition function of the $O(1)$ LG model on the honeycomb lattice, i.e., $Z_{O(1)}(x)$, with the loop fugacity $x=q^{2}/(1-r+r^{2})$, of which the critical point is $x_{c}=1/\sqrt{3}$Nienhuis (1982). Consequently, the norm of the LG ansatz maps to $Z_{O(1)}(q^{2}/(1-r+r^{2}))$ which becomes critical at $q_{c}=\sqrt{(1-r+r^{2})/\sqrt{3}}$ for a given $r$.

Appendix B optimal variational parameters in the first order ansatz

In this section, we present optimal variational parameters in the first order ansatz $|\psi_{1}(c_{1},c_{2})\rangle=\hat{Q}_{\rm LG}\hat{R}_{\rm DG}(c_{1},c_{2})|\Psi_{111}\rangle$, where $|\Psi_{111}\rangle$ is the product state of local magnetic (111)-state. In a similar way to the LG operator, the dimer gas (DG) operatorLee et al. (2019) $\hat{R}_{\rm DG}(c_{1},c_{2})$ is also efficiently represented in the tensor network and by the following building block tensor

[TABLE]

where we assign $c_{1}\,(c_{2})$ if the non-zero element comes from the intra-triangle (inter-triangle) bond. The dimension and direction of the virtual indices $i,j,$ and $k$ are the same as the ones of the $Q$-tensor in Eq. (6). To be more specific, on the site 8 in Fig. 1 (a) in the main text on which the inter-triangle bond is the $y$-bond, we put the building block tensor with

[TABLE]

Note that the variational parameters $c_{1}$ and $c_{2}$ determine the fugacity of the dimer on the intra-triangle bond and inter-triangle bond, respectively, whereas the fugacity of hole is set to unity. Let us reparametrize the variational parameters as follows: $\hat{R}_{\rm DG}(c_{1},c_{2})\rightarrow\hat{R}_{\rm DG}(\alpha,\beta)$ with

[TABLE]

Now, one can vary two variational parameters $\alpha$ and $\beta$ to minimize the energy expectation value of the Kitaev model on the star lattice. Using the corner transfer matrix renormalization group method, we measured the expectation values and then found the optimal $\alpha$ and $\beta$ at a given $\phi$ (see the main text) which are presented in Fig. 5 (a). The energy expectation values shown in Fig. 2(a) in the main text are obtained with the optimal parameters presented in Fig. 5 (a). Interestingly, the optimal $\alpha$, which determines the relative weight between dimers on the intra-triangle and inter-triangle bonds, increases linearly with the model parameter $\phi$, which determines the relative strength between the exchange couplings on the intra-triangle ($J$) and inter-triangle ($J^{\prime}$) bonds. Meanwhile, the parameter $\beta$, which determines the relative weight between the hole and dimers, is optimized in $0.2\pi<\beta<0.3\pi$. Note that one can also introduce complex fugacities, which lead to two more variational parameters, i.e., $(c_{1},c_{2})\rightarrow(c_{1}e^{i\eta_{1}},c_{2}e^{i\eta_{2}})$. However, we found that real fugacities always give the lowest energy throughout the model parameter $\phi$. Figure 5 (b) shows the second derivative of the energy expectation value of $|\psi_{1}\rangle$ in terms of $\phi$. Here, the transition point expected from $|\psi_{1}\rangle$ is a bit different from the exact one $\phi_{c}=\pi/3$. However, considering the fact that the curve is obtained by numerical differenciation twice, its accuracy and smoothness are quite remarkable.

Appendix C global flux

Here, we show that the LG and SG ansatze are not the eigenstates of the global flux operators, i.e., the Wilson loop operator $W_{\Gamma}=\prod_{i\in\Gamma}\sigma_{i}^{\alpha_{i}}$ along a non-contractible loop $\Gamma$ on a compact manifold where $\alpha=x,y$ and $z$ depending on the site $i$Kitaev (2006). To this end, we first note that, as shown in Ref. Lee et al. (2019), the multiplication of the Pauli matrices on the physical leg of local tensor of the LG operator is identical to the multiplication of the matrix

[TABLE]

and its conjugate on two virtual legs:

[TABLE]

Here, the tensor $Q_{ijk}^{ss^{\prime}}$ denotes the local tensor of the LG operator, and repeated indices are implicitly summed over. Above relation can be simply described by the following graphical representation:

[TABLE]

where the gray circle denotes the $Q$-tensor, black solid line stands for the physical leg (index $s$) and red, blue and yellow solid lines are virtual legs (indices $i,j$ and $k$) on the $x,y$ and $z$ bonds of the Kitaev model. Also, the red square attached on the physical leg denotes the Pauli matrix. Note that a physical leg $s^{\prime}$ is omitted in the graphical representation for simplicity. In what follows, using the above relation, we show how the LG and SG ansatze react on the action of the global flux operator $W_{\Gamma}$ and how to construct the eigenstates of $W_{\Gamma}$ using them. Let us consider the periodic boundary condition along the $y$-direction as defined in the main text and then act the operator $W_{\Gamma}$ along the $y$-direction as depicted in Fig. 6 (a). Now, using the relation in Eq. (14), one can show that the ring of the tensor network, where the operator $W_{\Gamma}$ is applied, has the following equalities,

[TABLE]

where the matrix $g=\sigma^{z}$ is the non-trivial element of the $Z_{2}$ invariant gauge group of the LG operatorLee et al. (2019). In the first equality, the relation in Eq. (14) is applied, and we use relations $vv^{\dagger}=1$ and $vv^{T}=-\sigma^{z}$. Finally, in the last equality, the invariant gauge symmetry is used, i.e., $g_{ii^{\prime}}g_{jj^{\prime}}g_{kk^{\prime}}Q_{i^{\prime}j^{\prime}k^{\prime}}^{ss^{\prime}}=Q_{ijk}^{ss^{\prime}}$. Therefore, applying the operator $W_{\Gamma}$ on our ansatze results in a different tensor network where a string of $g$, $G=\prod_{i=1}^{L_{y}}g$, along the $y$-direction is inserted in the original state as illustrated in Fig. 6 (b). One can easily notice that such a $G$ on the non-contractible loop can not be eliminated by a gauge transformation, and therefore the resulting state is not identical to the original state, i.e. Fig. 6 (a) = Fig. 6 (b). That is, our ansatze are not the eigenstate of the global flux operator. However, since $(W_{\Gamma})=1$ and $g^{2}=1$, it is easy to construct the eigenstate of $W_{\Gamma}$ using our ansatz, i.e., $|\psi_{\pm}\rangle=|\psi\rangle\pm|\psi_{G}\rangle$ where $|\psi\rangle$ denotes the LG or SG ansatz and $|\psi_{G}\rangle$ the $G$-inserted $|\psi\rangle$ along the $y$-direction as shown in Fig. 6 (b). Then, the state $|\psi_{\pm}\rangle$ is the eigenstate of $W_{\Gamma}$ with the eigenvalue or global flux number $\pm 1$.

Appendix D entanglement spectrums and conformal towers

In the main text, the entanglement spectrums are presented as a function of the momentum $k_{y}$ in each topological sector, and their degeneracy patterns are discussed. Here, we directly compare the entanglement spectrums and the Virasoro towers of the Ising conformal field theory. In Fig. 7 (a), we compare the entanglement spectrums of the string gas (SG) ansatz $|\psi_{1}\rangle$ optimized at $\phi=0.25\pi$ with the Virasoro characters of the primary operators having the conformal weights $\Delta=0,1/2$ and $1/16$Henkel (1999). As one can see, the spectrums in the sector $|\mathbb{I}\rangle$ (left panel) can be regarded as the sum of two Virasoro towers of $\Delta=0$ and $1/2$, while the ones in the sector $|\mathfrak{m}\rangle$ (right panel) match with the Virasoro tower of $\Delta=1/16$. Furthermore, their spacings and degeneracy patterns are in excellent agreement with the exact ones up to the seventh level. In Fig. 7 (b), the entanglement spectrums of the loop gas ansatz $|\psi_{0}\rangle$ optimized at $\phi=0.02\pi$ are presented. However, we could not identify their characteristics and relationship to the Ising conformal field theory, though the ansatz provides very accurate variational energy and entanglement entropy as shown in the main text.