Gapless Kitaev Spin Liquid to Classical String Gas through Tensor Networks

TL;DR

This paper introduces a tensor network framework for the gapless Kitaev spin liquid, revealing its hidden string gas structure and accurately modeling its properties without Majorana fermions.

Contribution

It presents a novel tensor network approach that captures key features of the Kitaev spin liquid, including symmetries, gauge structure, and phase transitions, with minimal parameters.

Findings

Achieves a highly accurate KSL ansatz with bond dimension D=8

Reveals the hidden string gas structure of the KSL

Provides insights into gap opening and non-Abelian phases under magnetic field

Abstract

We provide a framework for understanding the gapless Kitaev spin liquid (KSL) in the language of tensor network(TN). Without introducing Majorana fermion, most of the features of the KSL including the symmetries, gauge structure, criticality and vortex-freeness are explained in a compact TN representation. Our construction reveals a hidden string gas structure of the KSL. With only two variational parameters to adjust, we obtain an accurate KSL ansatz with the bond dimension D = 8 in a compact form, where the energy is about 0.007% higher than the exact one. In addition, the opening of gap and non-Abelian phase driven by a magnetic field are naturally understood in our construction.

Figures

Equations

$$\hat{\mathcal{H}}=-\sum _{⟨\alpha \beta {⟩}_{\gamma }}{J}_{\gamma }{\hat{\sigma }}_{\alpha }^{\gamma }{\hat{\sigma }}_{\beta }^{\gamma }$$

$$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 }},$$

$$(C_6{\hat{U}}_{C_6}){\hat{Q}}_{ijk}(C_6{\hat{U}}_{C_6}{)}^{†}={\hat{Q}}_{ijk}.$$

$${\hat{\sigma }}^x{\hat{Q}}_{ijk}$$

$${\hat{\sigma }}^y{\hat{Q}}_{ijk}$$

$${\hat{\sigma }}^z{\hat{Q}}_{ijk}$$

$$\text{includegraphics},$$

$$|T_{ijk}^0⟩\equiv {\hat{Q}}_{ijk}|(111)⟩.$$

$$\text{includegraphics}.$$

$$⟨{\psi }_0|{\psi }_0⟩=N_{\Gamma }\sum _{G\in \Gamma }{\left(\frac{1}{\sqrt{3}}\right)}^{l_G}=N_{\Gamma }\times Z_{O(1)}\left(\frac{1}{\sqrt{3}}\right),$$

$$g_{i{i}^{\prime }}{g}_{j{j}^{\prime }}{g}_{k{k}^{\prime }}|T_{i^{\prime }{j}^{\prime }{k}^{\prime }}^0⟩=|T_{ijk}^0⟩.$$

$$\text{includegraphics},$$

$${\hat{R}}_{ijk}={\zeta }_{ijk}({\hat{\sigma }}^x{)}^i({\hat{\sigma }}^y{)}^j({\hat{\sigma }}^z{)}^k.$$

$$\text{includegraphics}.$$

$$|T_{i_1{j}_1{k}_1}^1(\varphi )⟩={\hat{R}}_{ijk}(\varphi )|T_{i_0{j}_0{k}_0}^0⟩,$$

$$|T_{i_2{j}_2{k}_2}^2(\varphi ,\theta )⟩={\hat{R}}_{ijk}(\theta )|T_{i_1{j}_1{k}_1}^1(\varphi )⟩.$$

$${\hat{Q}}_{ijk}={\tau }_{ijk}({\hat{\sigma }}^x{)}^{1-i}({\hat{\sigma }}^y{)}^{1-j}({\hat{\sigma }}^z{)}^{1-k},$$

$${\tau }_{000}=-i,{\tau }_{011}={\tau }_{101}={\tau }_{110}=1.$$

$${\hat{U}}_{C_6}=\frac{1}{2}\left({\hat{\sigma }}^0+i{\hat{\sigma }}^x+i{\hat{\sigma }}^y+i{\hat{\sigma }}^z\right),{\hat{U}}_{\sigma }=\frac{i}{\sqrt{2}}\left({\hat{\sigma }}^x-{\hat{\sigma }}^y\right),$$

$${\hat{U}}_{C_6}{\hat{\sigma }}^x{\hat{U}}_{C_6}^{†}={\hat{\sigma }}^z,{\hat{U}}_{C_6}{\hat{\sigma }}^y{\hat{U}}_{C_6}^{†}={\hat{\sigma }}^x,{\hat{U}}_{C_6}{\hat{\sigma }}^z{\hat{U}}_{C_6}^{†}={\hat{\sigma }}^y,$$

$${\hat{U}}_{\sigma }{\hat{\sigma }}^x{\hat{U}}_{\sigma }^{†}=-{\hat{\sigma }}^y,{\hat{U}}_{\sigma }{\hat{\sigma }}^y{\hat{U}}_{\sigma }^{†}=-{\hat{\sigma }}^x,{\hat{U}}_{\sigma }{\hat{\sigma }}^z{\hat{U}}_{\sigma }^{†}=-{\hat{\sigma }}^z.$$

$${\hat{U}}_{C_6}{\hat{Q}}_{ijk}{\hat{U}}_{C_6}^{†}={\tau }_{ijk}{\hat{U}}_{C_6}({\hat{\sigma }}^x{)}^{1-i}({\hat{\sigma }}^y{)}^{1-j}({\hat{\sigma }}^z{)}^{1-k}{\hat{U}}_{C_6}^{†}={\tau }_{ijk}({\hat{\sigma }}^z{)}^{1-i}({\hat{\sigma }}^x{)}^{1-j}({\hat{\sigma }}^y{)}^{1-k},$$

$${\hat{U}}_{\sigma }{\hat{Q}}_{ijk}{\hat{U}}_{\sigma }^{†}={\tau }_{ijk}{\hat{U}}_{\sigma }({\hat{\sigma }}^z{)}^{1-i}({\hat{\sigma }}^x{)}^{1-j}({\hat{\sigma }}^y{)}^{1-k}{\hat{U}}_{\sigma }^{†}={\tau }_{ijk}(-{\hat{\sigma }}^y{)}^{1-i}(-{\hat{\sigma }}^x{)}^{1-j}(-{\hat{\sigma }}^z{)}^{1-k},$$

$$C_6({\hat{U}}_{C_6}{\hat{Q}}_{ijk}{\hat{U}}_{C_6}^{†})C_6^{-1}=C_6\left[{\tau }_{ijk}({\hat{\sigma }}^z{)}^{1-i}({\hat{\sigma }}^x{)}^{1-j}({\hat{\sigma }}^y{)}^{1-k}\right]C_6^{-1}={\tau }_{ijk}({\hat{\sigma }}^x{)}^{1-i}({\hat{\sigma }}^y{)}^{1-j}({\hat{\sigma }}^z{)}^{1-k}={\hat{Q}}_{ijk},$$

$$\sigma ({\hat{U}}_{\sigma }{\hat{Q}}_{ijk}{\hat{U}}_{\sigma }^{†}){\sigma }^{-1}=\sigma \left[{\tau }_{ijk}(-{\hat{\sigma }}^y{)}^{1-i}(-{\hat{\sigma }}^x{)}^{1-j}(-{\hat{\sigma }}^z{)}^{1-k}\right]{\sigma }^{-1}={\tau }_{ijk}(-{\hat{\sigma }}^x{)}^{1-i}(-{\hat{\sigma }}^y{)}^{1-j}(-{\hat{\sigma }}^z{)}^{1-k},$$

$$\mathcal{T}{\hat{Q}}_{ijk}\mathcal{T}=\{\begin{array}{ll}{\hat{Q}}_{ijk}& \mathrm{if}i+j+k=0\\ -{\hat{Q}}_{ijk}& \mathrm{if}i+j+k=2\end{array},$$

$$\mathcal{T}{\hat{Q}}_{ijk}\mathcal{T}=W_{i{i}^{\prime }}{W}_{j{j}^{\prime }}{W}_{k{k}^{\prime }}{\hat{Q}}_{i^{\prime }{j}^{\prime }{k}^{\prime }}\mathrm{with}W=\left(\begin{array}{cc}1& 0\\ 0& i\end{array}\right).$$

$${\hat{R}}_{ijk}={\zeta }_{ijk}({\hat{\sigma }}^x{)}^i({\hat{\sigma }}^y{)}^j({\hat{\sigma }}^z{)}^k,$$

$${\zeta }_{000}=1,{\zeta }_{100}={\zeta }_{010}={\zeta }_{001}=c.$$

$$(C_6{\hat{U}}_{C_6}){\hat{R}}_{ijk}({\hat{U}}_{C_6}^{†}{C}_6^{-1})={\zeta }_{ijk}(C_6{\hat{U}}_{C_6})\left[({\hat{\sigma }}_x{)}^i({\hat{\sigma }}_y{)}^j({\hat{\sigma }}_z{)}^k\right]({\hat{U}}_{C_6}^{†}{C}_6^{-1})={\zeta }_{ijk}{C}_6\left[({\hat{\sigma }}_z{)}^i({\hat{\sigma }}_x{)}^j({\hat{\sigma }}_y{)}^k\right]C_6^{-1}={\hat{R}}_{ijk},$$

$$(\sigma {\hat{U}}_{\sigma }){\hat{R}}_{ijk}({\hat{U}}_{\sigma }^{†}{\sigma }^{-1})$$

$$={\zeta }_{ijk}(-{\hat{\sigma }}_x{)}^i(-{\hat{\sigma }}_y{)}^j(-{\hat{\sigma }}_z{)}^k.$$

$$\displaystyle\mathcal{T}\hat{R}_{ijk}\mathcal{T}=\begin{cases}W_{ii^{\prime}}W_{jj^{\prime}}W_{kk^{\prime}}\hat{R}_{i^{\prime}j^{\prime}k^{\prime}}&{\rm if}\quad$c${\rm\,\,is\,\,real}\\ \hat{R}_{ijk}&{\rm if}\quad$c${\rm\,\,is\,\,imaginary}\end{cases}$$

Full text

Introduction- Quantum spin liquids (QSL) represent a state of quantum matter which is not characterized by any local order parameters even at zero temperature. These novel states are expected to exhibit long-range entanglement leading to the topological order and fractionalized excitationsSavary and Balents (2017). For example, the nearest-neighbor resonating valence bond (nnRVB) statesAnderson (1987) have been extensively studied as variational wavefunctions for the ground states of frustrated quantum magnetsMoessner and Sondhi (2001); Wen (2002). Indeed, the nnRVB states are topologically orderedMoessner and Sondhi (2001); Wen (2002); Zhou et al. (2017); Poilblanc et al. (2012); Lee et al. (2017) and support spinon excitations carrying the fractionalized quantum numberWen (2002); Zhou et al. (2017). However, since they are not exact ground states of the antiferromagnetic Heisenberg model, variational methods with the nnRVB states have been employed to search for true ground statesKaneko et al. (2014); Jiang et al. (2016); Iqbal et al. (2016); Mei et al. (2017). The Haldane phase, which is also known as a symmetry-protected topological phase, is another fascinating phase one can find in the $S=1$ quantum spin chain. The novel character that discriminates the Haldane phase from trivial gapped states was most clearly revealed by the discovery of Affleck–Kennedy–Lieb–Tasaki (AKLT) model and its exact ground state or AKLT stateAffleck et al. (1987). The compact representation of AKLT stateKlümper et al. (1993) provided a new insight into the Haldane phase. In addition, it was subsequently used in a variety of contexts for variational calculations on the quantum spin systemsKlümper et al. (1993); Kolezhuk and Mikeska (1998).

Kitaev honeycomb model (KHM) is an exactly soluble model which exhibits gapless and gapped KSL ground states with fractionalized excitationsKitaev (2006). Recent successful realizations of Kitaev materialsKhaliullin (2005); Jackeli and Khaliullin (2009); Plumb et al. (2014); Zhou et al. (2016); Banerjee et al. (2016); Trebst (2017); Singh and Gegenwart (2010) triggered a burst of theoretical investigations on KHM and its extensionsKimchi and You (2011); Catuneanu et al. (2018); Gohlke et al. (2018). In addition, due to the non-Abelian phase of the KSL driven by the magnetic fieldJiang et al. (2011); Zhu et al. (2018) and its potential application to quantum computation, it has been a subject of active research for the last decade. We refer the readers to Ref.Nayak et al. (2008) for an exhaustive list of relevant literature. The TN methods have been also employed to represent the KSLSchmoll and Orús (2017); Osorio Iregui et al. (2014). However, the Majorana basis TN requires a three-dimensional structure which makes it impractical as a tool for the numerical optimizationSchmoll and Orús (2017). On the other hand, the spin basis TN study, which was done with the computationally expensive optimization, suffers from undesirable breaking of symmetriesOsorio Iregui et al. (2014). In this Letter, we provide a compact TN representation for KHM that is defined with the spin basis and retains various symmetries.

Model. The KHM is defined asKitaev (2006)

[TABLE]

where $\langle\alpha\beta\rangle_{\gamma}$ stands for a pair on the $\gamma\,(=x,y,z)$ links connecting sites $\alpha$ and $\beta$ as depicted in Fig. 1 (a). As demonstrated in Kitaev’s seminal workKitaev (2006), the Hamiltonian commutes with the so-called flux operators defined on all hexagonal plaquette ($p$): $[\hat{\mathcal{H}},\hat{W}_{p}]=0$ with $\hat{W}_{p}=\hat{\sigma}^{x}_{1}\hat{\sigma}^{y}_{2}\hat{\sigma}^{z}_{3}\hat{\sigma}^{x}_{4}\hat{\sigma}^{y}_{5}\hat{\sigma}^{z}_{6}$ where the site indices 1-6 are defined in Fig. 1 (a). Therefore, the Hilbert space is sectorized by each flux number $\{\hat{W}_{p}=\pm 1\}$. Even further, in each sector, the KHM becomes a noninteracting Majorana fermion hopping model in the background of static $Z_{2}$ gauge fields. The ground states live in the vortex-free sector ($\hat{W}_{p}=+1$ for all $p$), which form a critical KSL phase around the isotropic point ($J_{x}=J_{y}=J_{z}=\pm 1$). In this Letter, we only consider the isotropic point at which the KHM is invariant under the following symmetry transformations: $C_{6}\hat{U}_{C_{6}}$ and $\sigma\hat{U}_{\sigma}$, where $\hat{U}_{C_{6}}=\bigotimes_{\alpha}\left(\hat{\sigma}^{0}_{\alpha}+i\hat{\sigma}^{x}_{\alpha}+i\hat{\sigma}^{y}_{\alpha}+i\hat{\sigma}^{z}_{\alpha}\right)/2$, $\hat{U}_{\sigma}=\bigotimes_{\alpha}\left(\hat{\sigma}^{x}_{\alpha}-\hat{\sigma}^{y}_{\alpha}\right)/\sqrt{2}$, and $C_{6}$, $\sigma$ respectively denote the $60^{\circ}$ spatial rotation and inversion as depicted in Fig. 1 (a). One can easily verify $[C_{6}\hat{U}_{C_{6}},\hat{\mathcal{H}}]=0=[\sigma\hat{U}_{\sigma},\hat{\mathcal{H}}]$. Also, the KHM is invariant under the time-reversal and translational symmetries.

Tensor network representation. We employ the tensor product state (TPS) representation Verstraete et al. (2008). Since the KSL is a zero-flux state You et al. (2012), we reasonably assume the TPS to be translationally symmetric. Additionally, we assume that the tensor does not depend on the sublattice, and therefore our ansatz is rewritten as $|\psi\rangle={\rm tTr}\prod_{\alpha}|T_{i_{\alpha}j_{\alpha}k_{\alpha}}\rangle,$ where tTr stands for the tensor trace or contraction of all virtual indices $\{i_{\alpha},j_{\alpha},k_{\alpha}\}$, $\alpha$ labels the site, $|T_{ijk}\rangle\equiv\sum_{s}T_{ijk}^{s}|s\rangle$ with $s$ being the local quantum number. Its graphical illustration is presented in Fig. 1, where the black open leg denotes the physical degrees of freedom. In what follows, we will construct the local tensor $|T_{ijk}\rangle$ (we identify $T_{ijk}^{s}$ and $|T_{ijk}\rangle$ and call both “tensor” hereafter) with consideration for the symmetries, the vortex-free condition and gauge structure. In the main text, we only discuss the ferromagnetic model $(J_{\gamma}=1)$, since the antiferromagnetic one is a trivial generalization and discussed in the supplementary material (SM) Lee et al. (2019).

Zeroth order tensor. We begin with introducing a bond dimension $D=2$ tensor product operator (TPO) referred to as the loop gas (LG) operator, $\hat{Q}_{\rm LG}={\rm tTr}\prod_{\alpha}Q_{i_{\alpha}j_{\alpha}k_{\alpha}}^{ss^{\prime}}|s\rangle\langle s^{\prime}|$ with a building block tensor

[TABLE]

which is depicted in Fig. 1 (c). The virtual indices $i,j$ and $k$ range from 0 to 1 ($D=2$), and non-zero elements of $\tau$-tensor are $\tau_{000}=-i$ and $\tau_{011}=\tau_{101}=\tau_{110}=1$. To simplify the notation, we define a local operator $\hat{Q}_{ijk}=\sum_{s,s^{\prime}}Q_{ijk}^{ss^{\prime}}|s\rangle\langle s^{\prime}|$. One can verify in the local tensor level that the LG operator respects the symmetries of KHM. For instance, applying $C_{6}\hat{U}_{C_{6}}$ on the $\hat{Q}$-tensor leaves it intactly, i.e.,

[TABLE]

Here, we use the facts that the $\hat{U}_{C_{6}}$-transformation rotates the spin, i.e., $U_{C_{6}}\hat{\sigma}^{x,y,z}_{\alpha}U_{C_{6}}^{\dagger}=\hat{\sigma}^{z,x,y}_{\alpha}$, while the $C_{6}$-rotation permutes the virtual indices as followse $C_{6}\circ(ijk)=(kij)$. Therefore, the resulting LG operator is invariant under the $(C_{6}\hat{U}_{C_{6}})$-transformation, and other symmetries of KHM can be shown in a similar wayLee et al. (2019). Note that the $Q$-operator satisfies the following relation

[TABLE]

with a matrix $v$, of which non-zero elements are $v_{01}=i$ and $v_{10}=1$, acting on the virtual bonds. Repeated indices are summed over, except where explicitly stated otherwise. Using Eq. (3), one can verify a relation $\hat{W}_{p}\hat{Q}_{\rm LG}=\hat{Q}_{\rm LG}\hat{W}_{p}=\hat{Q}_{\rm LG}$. To be more specific, the invariance of a patch of $\hat{Q}_{\rm LG}$ under the action of $\hat{W}_{p}$ can be shown as follows:

[TABLE]

where the connected green squares denote $\hat{W}_{p}$, and a physical leg is omitted for simplicity. Here, Eq. (3) and $v^{\dagger}v=1$ are used in the first and second equalities, respectively. This remarkable relation guarantees a quantum state $|\psi\rangle=\hat{Q}_{\rm LG}|\phi\rangle$, where $|\phi\rangle$ is an arbitrary state, being vortex-free and thus non-magnetic. Notice that the LG operator is identical to the projector $\prod_{p}(1+\hat{W}_{p})/2$ up to a normalization factor.

Regarding the $(C_{6}\hat{U}_{C_{6}})$-symmetry, let us apply $\hat{Q}_{\rm LG}$ on a product state $|\phi_{0}\rangle=\otimes_{\alpha}|(111)\rangle_{\alpha}$, where $|(111)\rangle$ denotes a spin aligned along (1,1,1) direction: $\langle(111)|\,\vec{\sigma}\,|(111)\rangle=(1,1,1)/\sqrt{3}$. Note that the ansatz $|\phi_{0}\rangle$ is a classical ground state respecting the $(C_{6}\hat{U}_{C_{6}})$-symmetry. Now, we define a quantum state $|\psi_{0}\rangle\equiv\hat{Q}_{\rm LG}|\phi_{0}\rangle$ which consists of a building block tensor

[TABLE]

We refer to it as zeroth order tensor. By virtue of the $\tau$-tensor in $\hat{Q}_{ijk}$, one can visualize the ansatz $|\psi_{0}\rangle$ as follows

[TABLE]

Here, the empty site stands for $|(111)\rangle$ state while the red loops denote the product of $\hat{\sigma}^{x}|(111)\rangle,\hat{\sigma}^{y}|(111)\rangle$ and $\hat{\sigma}^{z}|(111)\rangle$ states depending on the direction of loop on each site.

By computing the norm of the LG ansatz, we can show its criticality. To this end, we first note that the LG operator is hermitian as well as idempotentLee et al. (2019): $\hat{Q}_{\rm LG}^{\dagger}\hat{Q}_{\rm LG}=N_{\Gamma}\hat{Q}_{\rm LG}$ where $N_{\Gamma}$ is a total number of the loop configuration in the system. Using such properties and a simple identity $\langle(111)|\hat{\sigma}^{\gamma}|(111)\rangle=1/\sqrt{3}$, it is straightforward to show that the norm of $|\psi_{0}\rangle$ reads

[TABLE]

where $\Gamma$ denotes a set of all possible loop configurations and $l_{G}$ is a total length of loops in a configuration $G$. Also, $Z_{O(1)}(x)$ stands for the partition function of the classical $O(1)$ loop gas model with the fugacity $x$, which is exactly solvable and critical at $x_{c}=1/\sqrt{3}$ Nienhuis (1982). It indicates that the norm of $|\psi_{0}\rangle$ is exactly mapped into the partition function of the critical classical model which guarantees the criticality of $|\psi_{0}\rangle$Ardonne et al. (2004). In addition, the Ising conformal field theory (CFT) with the central charge $c=1/2$ is known to characterize the critical LG modelNienhuis (1982), which is consistent with the KSL of KHMMånsson et al. (2013); Lahtinen et al. (2014); Meichanetzidis et al. (2016).

The LG structure encoded in the $\tau$-tensor is useful in describing the vortex excitation of the KSL. To see this, we first note that the $\tau$-tensor is invariant under a gauge transformation $g=\hat{\sigma}^{z}$, i.e. $g_{ii^{\prime}}g_{jj^{\prime}}g_{kk^{\prime}}\tau_{i^{\prime}j^{\prime}k^{\prime}}=\tau_{ijk}$, and thus

[TABLE]

With a trivial gauge transformation $I_{2}$ being a two-dimensional identity matrix, they form a $Z_{2}$ invariant gauge group (IGG). String-like action of $g$ on links would twist the gauge fieldsKitaev (2006) along the string and hence create two vortices $\hat{W}_{p}=-1$ at both ends as demonstrated below:

[TABLE]

where $\pm 1$ in the hexagon denotes $\hat{W}_{p}$. One can explicitly showLee et al. (2019) such creation and move of fluxes using Eq. (3).

Finally, we measure the KHM energy (per bond) of $|\psi_{0}\rangle$ and obtain $E=-0.16349$ which is rather higher than the exact one $E_{\rm Kitaev}\simeq-0.19682$Kitaev (2006). Details in numerics will be discussed later. By construction, the LG ansatz $|\psi_{0}\rangle$ made of zeroth order tensor satisfies most of the physical constraints respected in the KSL [see SM for the time-reversal and $\sigma\hat{U}_{\sigma}$ symmetries] but is energetically far away from the exact solution. In what follows, we present a simple but effective TPO ($D=2$) applied to the LG ansatz which reduces the energy greatly without violating the constraints. We refer to it as the dimer gas (DG) operator.

Higher order tensors. The DG operator is defined by $\hat{R}_{\rm DG}={\rm tTr}\prod_{\alpha}\hat{R}_{i_{\alpha}j_{\alpha}k_{\alpha}}$ with

[TABLE]

Here, non-zero elements of $\zeta$-tensor are $\zeta_{000}=1$ and $\zeta_{100}=\zeta_{010}=\zeta_{001}=z$ with $i,j,k=0,1$, and $z$ is a real(or pure imaginary) variational parameter fixing the fugacity of a dimer. In this context, the dimer denotes the operator $\hat{\mathcal{H}}_{\alpha\beta}^{\gamma}$. Then, the DG operator can be interpreted as a sum of all possible dimer configurations, i.e., $\hat{R}_{\rm DG}(z)=\sum_{G\in\Gamma_{\rm D}}\hat{R}_{G}(z)$ where $\hat{R}_{G}(z)=\bigotimes_{\langle\alpha\beta\rangle_{\gamma}\in G}(z\hat{\mathcal{H}}_{\alpha\beta}^{\gamma})$ is defined for each dimer configuration $G$, and $\Gamma_{\rm D}$ is the set of all dimer configurations:

[TABLE]

Due to $[\hat{\mathcal{H}}_{\alpha\beta}^{\gamma},\hat{W}_{p}]=0$, it is obvious that $\hat{R}_{G}$ commutes with $\hat{W}_{p}$ for any $G$, and hence $\hat{R}_{\rm DG}$ does; $[\hat{R}_{\rm DG},\hat{W}_{p}]=0$. In fact, we can easily prove $[\hat{R}_{\rm DG},\hat{Q}_{\rm LG}]=0$ and that the DG operator respects all symmetries of the KSL Lee et al. (2019). Therefore, its multiplication to $|\psi_{0}\rangle$ does not contaminate the features of the KSL regardless of $z$. Moreover, it can be expressed as the polynomial function of the KHM Hamiltonian, which may be the reason why it improves the energy of the ansatz quite efficiently. The first key observation is that we can graphically represent Eq. (1) raised to the $n$-th power as the linear combination of elements of $\Gamma_{\rm D}$ as

Here, the number of sites in the system is assumed to be $2N$. The terms grouped with a coefficient $\alpha_{0}$ are the fully-packed configurations while the second ones are configuration with $N-2$ dimers. The terms on the second line have $q$-mers longer than dimer, e.g. trimer $\hat{\mathcal{H}}_{\alpha\beta}^{x}\hat{\mathcal{H}}_{\beta\gamma}^{y}$. All those terms with $q$-mers are canceled by the anticommutativity of Pauli matrices, and thus $\beta=0$. Note that the configurations with the same number of dimers share the coefficient which resembles the $\hat{R}_{\rm DG}$. Then, one can recast it as $\hat{R}_{\rm DG}=\sum_{M}h_{M}\mathcal{\hat{H}}^{N-2M}$ with proper coefficients $h_{M}$. Note that our approach is not a perturbative oneVanderstraeten et al. (2017).

Now, we define the $n$-th order ansatz as $|\psi_{n}(\{z_{i}\})\rangle=\prod_{i=1}^{n}\hat{R}_{\rm DG}(z_{i})|\psi_{0}\rangle$ having $n$ complex variational parameters. Due to the application of the DG operator, the ansatz $|\psi_{n}\rangle$ can be interpreted as a string gas state which is a linear superposition of string configurations. The string configuration consists of open and closed strings, connected loops and string-connected loops as depicted in Fig. 2. The building block tensor of $|\psi_{n}\rangle$, referred to as the $n$-th order tensor, is obtained by applying the $\hat{R}_{ijk}$-operator $n$-times on the zeroth order tensor in Eq. (14). The bond dimension scales as $D=2^{n+1}$. Note that the LG feature or $\tau$-tensor in the zeroth tensor is inherited by all higher order tensors. Furthermore, the $\hat{R}_{ijk}$-operator is invariant only under the trivial gauge transformation, and thus its action does not enlarge the $Z_{2}$ IGG of which the non-trivial element is simply $g_{n}=I_{2^{n}}\otimes\hat{\sigma}^{z}$. In contrast to the zeroth order case, the norm of $|\psi_{n}\rangle$ does not map to the LG model. However, by employing the loop TN renormalizationYang et al. (2017), we numerically prove that the $n$-th order ansatz are also critical and characterized by the Ising CFT as summarized in Fig. 2. We also present the best variational energies at each order in Fig. 2, and the details are given below.

Variational ansatz. Now, we turn on and tune variational parameters to obtain a better ansatz than the zeroth one. We parametrize the $\zeta$-tensor as follows: $\zeta_{000}=\cos\phi$, $\zeta_{100}=\zeta_{010}=\zeta_{001}=\sin\phi$, and hence $\hat{R}_{ijk}(z)\rightarrow\hat{R}_{ijk}(\phi)$. For measuring the energy, we employ the corner transfer matrix renormalization group method (CTMRG)Nishino and Okunishi (1996); Orús and Vidal (2009); Corboz et al. (2010) of which accuracy is controlled by the dimension $\chi$ of CTM. The parallel C++ library mptensorMorita et al. (2016–) is utilized to perform CTMRG.

Let us begin with the first order ansatz $|\psi_{1}(\phi)\rangle$ and its building block tensor

[TABLE]

where $i_{n}=2^{n}i+i_{n-1}$ and $i=0,1$. The energy of $|\psi_{1}(\phi)\rangle$ is presented in Fig. 3 (a) as a function of $\phi$, of which the lowest $E=-0.19644$ is found at $\phi=0.24\pi$. Here, we fix $\chi=64$. It is remarkable that the first order tensor ($D=4$) already attains such a small error of $0.2\%$. Furthermore, we perform the loop TN renormalization to evaluate the norm of $|\psi_{1}(\phi)\rangle$ and extract the central charge and scaling dimensions presented in Fig. 3 (b)[see SM for more details]. All those are in excellent agreement with the ones of Ising CFT, and therefore our ansatze are critical and belong to the same universality class.

To obtain an ansatz even closer to the KSL, we consider the second order ansatz $|\psi_{2}\rangle$ and tensor ($D=8$):

[TABLE]

Its overall energy landscape is shown in Fig. 4 (a) as functions of $(\phi,\theta)$ and minimized at $(\phi,\theta)=(0.342\pi,0.176\pi)$. After an additional scaling with respect to $\chi$Lee et al. (2019), we obtain the best variational ansatz with $E=-0.19681$ which is only $0.007\%$ higher than the exact one. Also, using the environment tensorsLee and Kawashima (2018), the five largest correlation lengths ($\xi_{i}$) are extracted and shown in Fig. 4 (b), which are diverging with $\chi$. Analogous figure is shown for $\psi_{1}$ in SM. Therefore, we reasonably conclude that the ansatz made of higher order tensors form a family of gapless states which we believe are smoothly connected to the exact KSL and, as a series, converge to it.

Further, we foundLee et al. (2019) that applying the (111)-direction magnetic field drives the ansatze into the gapped phaseKitaev (2006). We speculate that these gapped ansatze host non-Abelian anyonic excitations. The description of the non-Abelian and Abelian topological phases under the LG and SG schemes is an interesting question, and now further study is in progressLee et al. .

Conclusion. Based on the physical and gauge symmetries and the vortex-free condition, we have constructed the compact TN representation, which generates a family of KSL-like states sharing the features of the KHM ground state. In this sense, the ansatze given in this study are analogous to the AKLT state as a member of the Haldane states or the RVB state as an ansatz of frustrated quantum magnetsAffleck et al. (1987); Wen (2002). Under this scheme, the string gas structure of the KSL comes in sight clearly which offers a novel viewpoint for the KSL and its physics. It also provides an intuitive picture for the KSL in the spin language without referring to the Majorana fermion, which has never been provided before. There are many generalizations that one can envision as well as concrete open questions involving the LG and SG ansatze, e.g., general LGs having larger internal degrees of freedom and their parent Hamiltonians. The relation between the general LGs and the string-net statesFendley (2008) is another interesting question to ask. We also find that the ansatz discussed in the present Letter provides a good initial state for variational method for the KHM with the magnetic fieldKaneko et al. . Further, for the anisotropic KHM, one can choose the initial magnetic state which differs from the state $|(111)\rangle$ and introduce a bond-dependent dimer fugacity as additional variational parameters to optimize the modelLee et al. . Therefore, we expect our work could furnish a better understanding of KSL and its neighboring phases observed in the Kitaev quantum magnets such as $\alpha$-RuCl3 Banerjee et al. (2016, 2018) and studied theoretically in extended KHMsJackeli and Khaliullin (2009); Chaloupka et al. (2010); Singh and Gegenwart (2010); Kimchi and You (2011); Gohlke et al. (2018). Using two variational parameters, the accuracy of $0.007\%$ in energy is obtained, which has never been achieved by other numerical optimizationsPhien et al. (2015); Corboz (2016); Vanderstraeten et al. (2016). This high accuracy, together with the observed systematic convergence, leads us to believe that the present scheme not only correctly captures the essence of KSL physics, but also provides a new direction for quantitatively accurate description of quantum spin liquids.

I Symmetries of the Loop Gas Operator

In the main text, we define the loop gas (LG) operator $\hat{Q}_{\rm LG}={\rm tTr}\prod_{\alpha}\hat{Q}_{i_{\alpha}j_{\alpha}k_{\alpha}}$ with the building block tensor

[TABLE]

where $i,j,k=0,1$, and the non-zero elements of $\tau$-tensor are

[TABLE]

We consider unitary operators

[TABLE]

which transform the Pauli matrices in the following way:

[TABLE]

Let us see how the LG operator transforms locally under the $\hat{U}_{C_{6}}$ and $\hat{U}_{\sigma}$ rotations:

[TABLE]

where Eq. (4) is used. Now, we consider the $C_{6}$ spatial rotation and $\sigma$ reflection transformations defined in Fig. (1) in the main text. Such lattice symmetry transformations permute the virtual indices: $C_{6}\circ(ijk)=(kij)$ and $\sigma\circ(ijk)=(jik)$. Let us apply those transformations on Eq. (5),

[TABLE]

where we use the fact that the $\tau$-tensor is fully symmetric under any permutation. The first equation in Eq. (6) directly indicates that the operator $\hat{Q}_{ijk}$, and thus the LG operator, is invariant under the $(C_{6}\hat{U}_{C_{6}})$-transformation, $(C_{6}\hat{U}_{C_{6}})\hat{Q}_{\rm LG}(C_{6}\hat{U}_{C_{6}})^{-1}=\hat{Q}_{\rm LG}$. Note that the length of any loop on the honeycomb lattice is even, which indicates that the extra minus signs in the second equation in Eq. (6) are redundant. Therefore, the LG operator remains invariant under the $(\sigma\hat{U}_{\sigma})$-transformation: $(\sigma\hat{U}_{\sigma})\hat{Q}_{\rm LG}(\sigma\hat{U}_{\sigma})^{-1}=\hat{Q}_{\rm LG}$. Now, we consider the time-reversal transformation $\mathcal{T}$ which transforms the operator $\hat{Q}_{ijk}$ as follows:

[TABLE]

where $\mathcal{T}\hat{\sigma}^{\gamma}\mathcal{T}=-\hat{\sigma}^{\gamma}$ is used. Even though the $Q$-tensor is not symmetric under the time-reversal transformation, an additional gauge transformation can restore its original form, i.e.,

[TABLE]

Therefore, the LG operator is time-reversal symmetric: $\mathcal{T}\hat{Q}_{\rm LG}\mathcal{T}=\hat{Q}_{\rm LG}$. By construction, the translational symmetry is respected in the LG operator, and therefore it keeps all symmetries of the isotropic Kitaev honeycomb model.

II Details on the Dimer gas operator

As demonstrated in the main text, in order to reduce the energy, we define the DG operator $\hat{R}_{\rm DG}={\rm tTr}\prod_{\alpha}\hat{R}_{i_{\alpha}j_{\alpha}k_{\alpha}}$ with

[TABLE]

and $i,j,k=0,1$, and the non-zero elements of $\zeta$-tensor are

[TABLE]

The constant $c$ is a variational parameter. As shown in the previous section, let us first consider the $C_{6}\hat{U}_{C_{6}}$-symmetry, i.e., $\hat{R}_{ijk}\rightarrow(C_{6}\hat{U}_{C_{6}})\hat{R}_{ijk}(C_{6}\hat{U}_{C_{6}})^{-1}$:

[TABLE]

where we used the fact that the $\zeta$-tensor is invariant under any permutation in the indices. Above relation implies $(C_{6}\hat{U}_{C_{6}})\hat{R}_{\rm DG}(C_{6}\hat{U}_{C_{6}})^{-1}=\hat{R}_{\rm DG}$. Next, under the $(\sigma\hat{U}_{\sigma})$-symmetry, the $\hat{R}_{ijk}$ transforms as

[TABLE]

Here, although an extra minus sign appears, it will be canceled after after contraction since the dimer is a two-site object. Consequently, the DG operator is symmetric under the $(\sigma\hat{U}_{\sigma})$-transformation: $(\sigma\hat{U}_{\sigma})\hat{R}_{\rm DG}(\sigma\hat{U}_{\sigma})^{-1}=\hat{R}_{\rm DG}$. As for the time-reversal symmetry, one can show

[TABLE]

Therefore, in any case, the DG operator is invariant under the time-reversal transformation: $\mathcal{T}\hat{R}_{\rm DG}\mathcal{T}=\hat{R}_{\rm DG}$. Consequently, if one can apply the PEPO on some ansatz respecting the symmetries, then the resulting state is also guaranteed to satisfy those symmetries.

III The time-reversal symmetry and $(\sigma\hat{U}_{\sigma})$-symmetry of ansatze

The zeroth order ansatz $|\psi_{0}\rangle$ is obtained by contracting the zeroth order tensor

[TABLE]

Let us see how it transforms under the $(\sigma\hat{U}_{\sigma})$-symmetry and time-reversal symmetry $\mathcal{T}$,

[TABLE]

where $|(-1,-1,-1)\rangle$ denotes a spin aligned along $(-1,-1,-1)$ direction: $\langle(-1,-1,-1)|\,\vec{\sigma}\,|(-1,-1,-1)\rangle=(-1,-1,-1)/\sqrt{3}$. Here, relations $\mathcal{T}|(111)\rangle=e^{i\pi}|(-1,-1,-1)\rangle$, $\hat{U}_{\sigma}|(111)\rangle=e^{i\pi/4}|(-1,-1,-1)\rangle$ and $(\mathcal{T})^{2}=-1$ are used. Note that one cannot restore $\mathcal{T}|T^{0}_{ijk}\rangle$ and $\sigma\hat{U}_{\sigma}|T^{0}_{ijk}\rangle$ to $|T^{0}_{ijk}\rangle$ by applying a gauge transformation. Consequently, the zeroth order tensor does not ensure the resulting state $|\psi_{0}\rangle$ to be time-reversal symmetric and $(\sigma\hat{U}_{\sigma})$-symmetric. However, it is invariant under the combination of $(\sigma\hat{U}_{\sigma})$ and $\mathcal{T}$ transformations, i.e., $(\sigma\hat{U}_{\sigma})\,\mathcal{T}|\psi^{0}\rangle=e^{i\theta}|\psi^{0}\rangle$. More precisely, the zeroth order tensor is transformed as follows

[TABLE]

where Eqs. (6) and (7) are used. The overall phase does not affect the resulting state, and therefore the zeroth order ansatz is invariant under the $(\sigma\hat{U}_{\sigma})\mathcal{T}$-transformation.

Even though the building block tensor does not guarantee the $\mathcal{T}$ and $\sigma\hat{U}_{\sigma}$ symmetries, those symmetries might be restored in a larger unit-cell. In order to carve this out, we measure the fidelity between the states $|\psi_{0}\rangle$ and its transformed one $O|\psi_{0}\rangle$ on a torus geometry with size $(L_{x},L_{y})$, where $O=\mathcal{T},\sigma\hat{U}_{\sigma}$ and $(\sigma\hat{U}_{\sigma})\mathcal{T}$, and $(L_{x},L_{y})$ denotes the number of unit-cell on horizontal and vertical directions. The results are shown in Fig. 1. As one can see, the transformed state becomes orthogonal to each other, and thus not symmetric under the transformations. Then, how can we make them symmetric? In fact, Eq. (15) provides a quick cure to construct an ansatz respecting both the $\mathcal{T}$ and $\sigma\hat{U}_{\sigma}$ symmetries by doubling the bond dimension:

[TABLE]

where $i=2i_{1}+i_{0}$ and

[TABLE]

The resulting state from above tensor is simply $|\widetilde{\psi}_{0}\rangle=|\psi_{0}\rangle+|\mathcal{T}\psi_{0}\rangle$ and thus, by construction, invariant under the time-reversal operation. The $\sigma\hat{U}_{\sigma}$ symmetry is also guaranteed as shown below

[TABLE]

which is restored to its original form by flipping $0\leftrightarrow 1$ in $i_{1},j_{1},k_{1}$ with a proper coefficient to eliminate the minus sign in the second line. It can be simply done by applying a gauge trasformaion

[TABLE]

such that

[TABLE]

This procedure also applies to the higher order tensors, and therefore one can always construct a higher order ansatz, which is time-reversal and $(\sigma\hat{U}_{\sigma})$-symmetric, by doubling the bond dimension.

IV Manipulation of $Z_{2}$ vortices

As mentioned in the main text, one can create the $Z_{2}$ vortices by acting the non-trivial elements of $Z_{2}$ invariant guage group (IGG) $g$ on the virtual bonds. Let us see how it happens. For simplicity, we only consider the zeroth order tensor with $g_{0}=\hat{\sigma}^{z}$, but its generalization to the general case is straightforward. In the main text, we explicitly showed that a plaquette patch of the LG operator is invariant under the action of flux operator using Eq. (3) in the main text. Let us consider the same tensor network except that the $g_{0}$ is inserted between two sites as depicted below:

[TABLE]

where the connected green squares denote the flux operator while

[TABLE]

In the last equality, we use $v^{\dagger}v=1$ and $v^{\dagger}\hat{\sigma}^{z}v=-\hat{\sigma}^{z}$. Therefore, the resulting state is an eigenstate of the flux operator with the eigenvalue $-1$. The $Z_{2}$ vortex is created. Similarly, another vortex is created on the opposite plaquette covering the $g_{0}$-inserted bond as depicted in Fig. 2 (a). Inserting another $g_{0}$ in the plaquette, the minus sign will be canceled such that the vortex is removed. But, on the opposite plaqutte of newly $g_{0}$-inserted bond, another $Z_{2}$ vortex is created. In other words, the vortex can be moved from a plaquette to another one as demonstrated in Fig. 2 (b). The $\zeta$-tensor embedded in the DG operator is invariant only under the trivial gauge transformation $e$, and thus the IGG of the DG operator is the trivial IGG. It indicates that the multiplication of the DG operator does not enlarge the IGG. Therefore, one can create and move the vortices with $g_{n}=I_{2^{n}}\otimes\hat{\sigma}^{z}$ with the higher order tensors and ansatze.

V Loop Configurations

V.1 Deformation of loops

As mentioned in the main text, the zeroth ansatz $|\psi_{0}\rangle$ has the quantum loop gas structure. In other words, the ansatz are represented by linear superpositions of all possible closed loop configurations with an equal weight. Here, the loop denotes the product of $\hat{\sigma}^{x}|(111)\rangle,\hat{\sigma}^{y}|(111)\rangle$ and $\hat{\sigma}^{z}|(111)\rangle$ states along the loop as depicted in Fig. 3 (a), while an empty state is simply $|(111)\rangle$. By applying the flux operator $W_{p}$, one can deform the loop configurations following a simple rule. First, one regards the flux operator $W_{p}$ on $p$ as a loop along the boundary of $p$ as demonstrated in Fig. 3 (b). Then, we put $W_{p}$ on a loop gas configuration and draw the loop along the plaquette. If some part of loops are overlapped, then we eliminate the overlapped fragments. This procedure does not break any loop but just deform or detour the loops. Also, it does not give any extra phase after the deformation. For example, applying the flux operator on $p$ in Fig. 3 (a), then the $\hat{\sigma}^{x}$ and $\hat{\sigma}^{z}$ apply on $3$ and $4$, respectively. Therefore, the local states on the sites $3$ and $4$ rotates as follows:

[TABLE]

The phases are canceled each other, and a part of loop on the bond between the sites 3 and 4 is deleted and extended to wrap the plaquette $p$. In Fig. 3 (c)-(f), we present some exemplary deformations of some loops by applying the flux operator. Using these local deformations, one can completely remove some loop configurations by applying the flux operators or create them from empty configurations.

V.2 Norm of $|\psi_{0}\rangle$

The norm of wavefunction contains important informations on the low-lying excitations in the system. In this subsection, we compute the norm of $|\psi_{0}\rangle$ state. It is easy to see that an inner product between the loop free and an arbitrary loop configurations is simply

[TABLE]

where $|G\rangle$ and $|0\rangle=\prod|(111)\rangle$ denote respectively the arbitrary loop configuration and loop free configuration and $l_{G}$ is the total length of loops in the configuration $G$. Here, we used $\langle(111)|\sigma^{\gamma}|(111)\rangle=1/\sqrt{3}$. As explained in the previous subsection, any loop configuration can be obtained by a product of flux operators, i.e.,

[TABLE]

where $p\in G$ denotes a proper choice of plaquettes to create the $|G\rangle$ from the loop free configuration. Then, the norm of $|\psi_{0}\rangle$ is rewritten as

[TABLE]

We use the fact that $(W_{p})^{2}=1$ in the third equality, and Eq. (24) is substituted in the last equality. The norm of $|\psi_{0}\rangle$ turns out to be the partition function of the $O(n)$ loop gas model at a critical point as presented in the main text.

V.3 Norm of $|\psi_{1}\rangle$

The norm of $|\psi_{1}(\phi)\rangle$ does not map to the exactly solvable point of the $O(n)$ loop gas model. Therefore, we employ the loop tensor network renormalization (LTNR) to numerically obtain the norm of $|\psi_{1}(\phi)\rangle$ and extract the central charge $c$ and scaling dimensions $\Delta_{i}$. Results for $\phi=0,0.125\pi$ and $0.25\pi$ are shown in Fig. 4 as a function of the real space renormalization step (RG step). Here, the bond dimension of LTNR is fixed to $\chi=32$. The number of iteration for loop optimization varies up to 20 to find the best ansatz at each RG step. As shown in Fig. 4 (a), at $\phi=0$ where the tensor becomes zeroth order one, the conformal data match nicely with the exact values from Ising universality class and shows a stable behavior up to about 20 RG step. Although the accuracy of LTNR becomes less as increasing the fugacity of dimer operator $(\phi>0)$, we could obtain reasonable and consistent results up to around $\phi=0.25\pi$ as presented in Fig. 4 (b) and (c).

VI General dimer fugacity and energy land scapes

Generally, the variational parameter in the DG operator is allowed to be complex though it breaks the time-reversal symmetry. The results with only real parameters are shown in the main text. Here, the energy dependence on the complex coefficients is presented and briefly discussed. The first order ansatz $|\psi_{1}\rangle$ is in general defined as

[TABLE]

where the $\zeta$-tensor in the tensor $\hat{R}_{ijk}$ is parameterized as follows: $\zeta_{000}=\cos\phi$ and $\zeta_{100}=\zeta_{010}=\zeta_{001}=e^{i\theta}\sin\phi$. Here, the parameter $\theta$ is additionally introduced to give an arbitrary phase. Now, we should fix two parameters to find the energy minimum point, and the energy landscape is shown in Fig. 5 (a). The lowest energy is obtained at $(\phi,\theta)=(0.24\pi,0)$. Though we could not mathematically prove whether or not the lowest energy is found by real coefficients, we observe that real coefficients give the lowest energy at several parameter points. Now, let us consider the second order ansatz

[TABLE]

Two phase variables $\alpha$ and $\beta$ are introduced. Therefore, we have to fix four independent parameters. Therefore, we present, in Fig. 5 (b), the $\alpha$ and $\beta$ dependence of energy only at $(\phi,\theta)=(0.342\pi,0.176\pi)$ at which the lowest energy is measured.

VII Antiferromagnetic Kitaev Honeycomb Model

In the main text, we construct the ansatze for the ferromagnetic Kitaev model by applying the LG and DG operators on the classical ground state, $|\phi^{\rm fm}\rangle=\otimes_{\alpha}|(111)\rangle_{\alpha}$ where $\alpha$ runs over all sites. Following the same strategy, we prepare the classical ground state of the antiferromagnetic model:

[TABLE]

where $\alpha$ labels the unit-cell, $a$ and $b$ denote two different sublattices, respectively. Then, we apply the LG and DG operators on the state $|\phi^{\rm afm}\rangle$ as we did in the ferromagnetic model and fix the variational parameters in the DG operator to find the lowest energy ansatz. The zeroth order ansatz $|\psi_{0}^{\rm afm}\rangle=\hat{Q}_{\rm LG}|\phi^{\rm afm}\rangle$ gives the energy $E=-0.14746$ which is rather higher than the one obtained by the zeroth ansatz of the ferromagnetic model (see the main text). Now, let us see how the DG operator reduces the energy. Here, we only consider the first order ansatz

[TABLE]

where the $\zeta$-tensor in the tensor $\hat{R}_{ijk}$ is parameterized as follows: $\zeta_{000}=\cos\phi$ and $\zeta_{100}=\zeta_{010}=\zeta_{001}=e^{i\theta}\sin\phi$. The energy as functions of $\phi$ and $\theta$ is shown in Fig. 6 (a). In contrast to the ferromagnetic case, the lowest energy is found with a negative dimer fugacity or $\theta=\pi/2$. In Fig. 6 (b), the energy is presented as a function of $\phi$ with $\theta=0.5\pi$. Here, the lowest energy $E=-0.195356$ is obtained at $(\phi,\theta)=(0.25\pi,0.5\pi)$, and it is only $0.75\%$ higher than the exact one. Again, this energy is slightly higher than the one obtained for the ferromagnetic model but still surprisingly close to the exact one with only $D=4$. We also confirmed that our ansatz for the antiferromagnetic model exhibit critical behavior. We believe that one could obtain much better ansatz by applying another DG operator ($D=8$) as shown in the main text.

VIII Effect of the (111)-direction magnetic field

In the current representation of the KSL, it is naturally anticipated that the (111)-direction magnetic field opens the excitation gap and drives the KSL into the non-Abelian phaseKitaev (2006). Applying the field, the vortex-free condition is not required anymore, and therefore introducing a parameter $\alpha$ into the $\tau$-tensor in $\hat{Q}_{ijk}$ is allowedLee et al. (2019), such that $\tau_{011}=\tau_{101}=\tau_{110}=\alpha$ and $\tau_{000}=-i$. It is obvious that $\tau_{000}$ element generates the local magnetic state $|(111)\rangle$ while the others do $\hat{\sigma}^{\gamma}|(111)\rangle$. In the weak-field limit, one can control $\alpha$ to modify the weight of LG and reasonably choose $\alpha<1$, since the (111)-field prefers the state $|(111)\rangle$ rather than $\hat{\sigma}^{\gamma}|(111)\rangle$ states. Then, the norm of wavefunction maps to $Z_{O(1)}(\alpha/\sqrt{3})$ [Eq. (6) in the main text], where the model enters into a massive phaseDomb (2000). Consequently, the gap is opened by the magnetic field.

In order to find a better ansatz with $D=4$, we introduce two additional parameters in $|\psi_{1}(\phi)\rangle$, which tune the weight of local $\hat{\sigma}^{\gamma}|(111)\rangle$ states. To be more specific, let us explicitly write down the non-zero elements of the first order tensor $|T^{1}_{ijk}(\phi)\rangle=\hat{R}_{i_{1}j_{1}k_{1}}(\phi)\hat{Q}_{i_{0}j_{0}k_{0}}|(111)\rangle$ where $i=2i_{1}+i_{0}$. For simplicity, we redefine the first order tensor as $|T^{1}_{ijk}(\phi)\rangle=\cos\phi|A_{ijk}\rangle+\sin\phi|B_{ijk}\rangle$ with

[TABLE]

Note that the $A$-tensor is the same as the zeroth order tensor. Now, we assign parameters tuning the weight of $\hat{\sigma}^{\gamma}|(111)\rangle$ states in the $A$-tensor and $B$-tensor, respectively. In other words, the ansatz becomes dependent on three parameters: $|T_{ijk}^{1}(\phi,\alpha_{0},\alpha_{1})\rangle=\cos\phi|A_{ijk}(\alpha_{0})\rangle+\sin\phi|B_{ijk}(\alpha_{1})\rangle$ where the parameter $\alpha_{i}$ multiplied to the $\hat{\sigma}^{\gamma}|(111)\rangle$ states. For example, the non-zero elements of $A$-tensor are

[TABLE]

Note that the $(C_{6}\hat{U}_{C_{6}})$-symmetry is still valid even in the presence of the (111)-field. Therefore, one allows to introduce only a single parameter $\alpha_{1}$ in the $B$-tensor. Since the $(111)$-field prefers the local state $|(111)\rangle$ rather than $\hat{\sigma}^{\gamma}|(111)\rangle$, one may naively expect that reducing the parameters $\alpha_{i}$ from 1 helps lowering the energy. Indeed, we found that the energy is optimized at $(\phi,\alpha_{0},\alpha_{1})=(0.225\pi,0.825,0.95)$ with the energy $E=-0.19688$ which is competitive to the one obtained by the numerical optimization Kaneko et al. . The norm of the ansatz is not mapped into the LG model similar to $|\psi_{1}(\phi)\rangle$. Therefore, in order to show its gapped nature, we directly measure the most dominant correlation lengths of the ansatz using the environment tensor in CTMRGTakahashi (2005); Lee and Kawashima (2018). The result for ansatz with $B=0.01J$ is presented in Fig. 7 (a) as a function of $\chi$. As one can see, the correlation lengths converge to finite values with increasing $\chi$ indicating a finite gap in the ansatz. For comparison, those in the critical state $|\psi_{1}(\phi=0.24\pi)\rangle$ are shown in Fig. 7 (b), which exhibit diverging behavior with $\chi$. Therefore, the gapped ansatz in the presence of (111)-field can be reasonably obtained by giving some fugacity to the tensor element generating $\hat{\sigma}^{\gamma}|(111)\rangle$ state. In addition, it has been recently shownFendley (2008) that the gapped LG having non-trivial inner products between two configurations can be systematically mapped into string-net states describing non-Abelian anyonic excitationsFendley (2008). It strongly suggests that our ansatz belongs to a non-Abelian phase in the presence of (111)-field, which is consistent with the perturbation calculation using the Majorana fermion in Kitaev’s original work Kitaev (2006).

IX $\chi$-scaling of the variational energies

We provide the bond dimension of CTMRG, $\chi$, dependence of the variational ansatze shown in the main manuscript. In Figure 8, the scaling behavior of the variational energies of the LG state $|\psi_{0}\rangle$, the first order SG ansatz $|\psi_{1}(\phi=0.24\pi)\rangle$, and the second order SG ansatz $|\psi_{2}(\phi=0.342\pi,\theta=0.176\pi)\rangle$ are shown from left to right, respectively.