Tensor network wavefunction of $S=1$ Kitaev spin liquids

Recently there has been a great interest in understanding quantum spin liquid phases with varying spin magnitude, partly due to possible material realizations. A number of recent numerical computations suggest that the ground state of the S=1 Kitaev model may be a quantum spin liquid in analogy to the renowned $S$=$1/2$ model. On the other hand, the nature of the ground state remains elusive since the $S$=$1$ model is not exactly solvable in contrast to the $S$=$1/2$ model. In this work, we construct a tensor network ground state wavefunction for the S=1 Kitaev model, which is explicitly written in terms of physical spin operators. We explain how this class of wavefunctions can be successfully used for variational computations and compare the outcomes to known results on finite size systems. We establish the existence of distinct topological sectors on torus by constructing the minimally entangled states in the degenerate ground state manifold and evaluating topological entanglement entropy. Our results suggest that the ground state of the S=1 Kitaev model is a gapped quantum spin liquid with Z2 gauge structure and Abelian quasiparticles.


Introduction -
The exact solution of the S=1/2 Kitaev model on the honeycomb lattice, i.e., Kitaev spin liquid (KSL), is an important milestone in theoretical research on quantum spin liquids [1].It has firmly established the existence of quantum spin liquids in systems with local spin interactions and opened the door to possible quantum spin liquid states arising from the bonddependent competing interactions.Both aspects lead to recent intensive experimental and theoretical explorations of the Kitaev materials , where the Kitaev interaction is dominant and the bond-dependent interactions originate from strong spin-orbit coupling.
On the other hand, the nature may also allow the Kitaev's bond-dependent interactions for higher-spin magnitudes, as recently suggested by electronic structure computations [24].While the large-spin limit would simply correspond to the classical limit, the ground state of the Kitaev model for relatively small spin magnitudes, especially for S=1, has not been fully understood as the higher-spin Kitaev models on honeycomb lattice do not have exact solutions.It was pointed out early on that there exists a flux operator, W p , made of six spin operators on each plaquette for an arbitrary spin S, which commutes with the Hamiltonian [25].Hence the ground state must have a definite eigenvalue of the flux operator, just like the case of the S=1/2 model.Exact diagonalization of the S=1 Kitaev model on small system sizes show that the ground state is non-magnetic [26].Thermal pure quantum state approach at finite temperatures on small system sizes finds an entropy plateau, similar to the S=1/2 model [27].All of these indicate the possibility of a quantum spin liquid in the S=1 Kitaev model while the precise nature of the ground state is not known.
In this paper, we provide the tensor network (TN) representation of the ground state wavefunction for the S=1 Kitaev model, with both ferromagnetic and antiferromagnetic interactions between the S=1 local moments.
Our approach is inspired by the recent construction of the TN ground-state wavefunction of the S=1/2 Kitaev model [28], where the wavefunction is written in terms of physical spin variables, instead of the Majorana fermion operators in the exact solution.We start from the bond dimension D=2 TN representation of the wavefunction, which is an explicit eigenstate of the flux operator via the projector Q LG , or the Loop Gas (LG) operator.The resulting wavefunction can be regarded as the sum over loop gas configurations made of a series of local actions of Q LG along a loop of lattice sites.We show that the Z 2 gauge structure is automatically built in the wavefunction.Using the mapping between the norm of the wavefunction and the partition function of the classical loop gas, it is shown that the wavefunction is in a gapped state, which is corroborated by the direct computation of the correlation length from the transfer matrix.This wavefunction is further improved by applying an additional Dimer Gas (DG) operator [28] and imaginary time evolution (ITE), resulting in excellent variational energy.
Our approach also allows us to understand some partial information about topological properties of the ground state wavefunction on torus via the construction of minimally entangled states in the ground state manifold.Using this construction, we contrast the difference between the degenerate ground states of the S=1 and S=1/2 systems, elucidating the nature of the quantum spin liquid ground state in the S=1 Kitaev model.
Model and symmetry -Hamiltonian of the Kitaev model on the honeycomb lattice [1] is defined as where ij γ stands for the nearest neighboring sites, i and j, on the γ-bond with γ = x, y, z as depicted in Fig. 1 (a), and S γ is the spin-1 operator.One can verify that the Hamiltonian commutes with a flux operator where U γ = e iπS γ is the 180 • spin-rotation operator along the γ-direction, and sites 0-5 are defined in Fig. 1 (a).It indicates that the spin-1 model also exhibits the Z 2 gauge redundancy, and thus the Hilbert space can be sectorized by the combination of the flux number.
Loop gas states.
-In a similar fashion to the spin-1/2 model [28], one can define the LG operator Q LG for the spin-1 model in a bond dimension D=2 TN representation with a local tensor where λ, µ, ν = 0, 1 are the virtual indices.Here, non-zero elements of the tensor τ are only τ λµν = −i for λ=µ=ν=0 and 1 for λ+µ+ν=2.It will be shown later that this operator generates loop configurations made of local actions represented in Fig. 1 (b).One can examine the physical symmetries of Q LG at the local tensor level as done in Ref. [28] and also derive the Z 2 gauge symmetry: λ ,µ ,ν σ z λλ σ z µµ σ z νν Q λ µ ν = Q λµν .Furthermore, it is straightforward to show that the local Q-tensor obeys the following equations where σ x is the Pauli matrix.Note that the Q-tensor for the spin-1/2 model obeys similar relations, but with a non-Hermitian unitary matrix v, defined in Ref. [28], rather than σ x in the above relations.That is because U γ satisfies the commutator relation, [U γ , U γ ] = δ γγ , while the Pauli matrices follow the anti-commutator relation, {σ γ , σ γ } = 2δ γγ .In addition, as will be seen below, it brings about fundamental distinction between the S=1/2 and S=1 LG states including topological property, even though the difference (v ↔ σ x ) may not look critical at first sight.One can also show at the local TN level that Q LG projects any quantum state into the vortex-free sector using Eq. ( 1) or W p Q LG = Q LG .See supplemental material (SM) for more details [29].
Due to the symmetries at the isotropic point, we apply the LG operator onto a classical product state |0 ≡ ⊗ i |0 i where all spins align in the (111)-direction, i.e., 0|S LG model with the fugacity ζ=1/3 which is in the gapped phase of the model [30].It is worth noting that the S = 1/2 LG maps to the critical point ζ c = 1/ √ 3 [28].Utilizing the corner transfer matrix renormalization group (CTMRG), we have computed the energy expectation value (per site) of |ψ LG , E LG = −0.50562,which is far from the previous ones E 24site = −0.648obtained on a 24-site system [31] and E DMRG = −0.644obtained by the density matrix renormalization group method on a (48 × 3)-cylinder system [32].In what follows, we present two different ways to vary the LG state to obtain better variational energies while most physical properties remain intact: (i) applying the so-called dimer gas (DG) operator onto the LG state and optimize its variational parameters as proposed in Ref. [28], (ii) evolving the LG state through the imaginary time evolution (ITE) operator [33].
String gas states.-The DG operator R DG is also represented in a D = 2 TN with the following local tensor LG , which can be optimized to have the lowest energy [28].Here, n is referred to as the order of the SG state, and the bond dimension grows as D = 2 n+1 .The variational energy of the first order ansatz (D=4) is shown in Fig. 2 (a) as a function of the variational parameter φ.The lowest energy is found to be E SG1 =−0.6167 which is still significantly higher than the ones obtained by ED and DMRG.Note that the efficiency of the DG operator is not as good as the case of the S=1/2 model.We believe that this is because the S=1 DG operator cannot be written as the polynomial function of the Hamiltonian in contrast to the spin-1/2 model, where the first order SG state already gives 99.8% accuracy to the exact one energetically [28].We have also optimized the second order ansatz |ψ SG2 (φ 1 , φ 2 ) and obtained the best energy E SG2 =−0.6366 [see the SM [29] for more details].
Imaginary time evolution.-We perform the ITE by applying the two-site gate e τ KS γ i S γ j on every bond and updating the site tensors using the simple update (SU) iteratively [33].Note that it is almost impossible to find a KSL-like spin liquid, e.g., W p =1, with SU from a random initial state.To be more precise, the resulting states are mostly magnetically ordered and not vortexfree, W p <1, and even exhibit the strong initial state dependence.However, the initial LG state with a careful choice of D [34] makes the ITE quite stable and reliable.For instance, the ITE operator commutes with the LG operator, i.e., [e τ KS γ i S γ j , Q LG ] = 0, so that the ITE does not spoil the vortex-freeness.While the truncated singular value decomposition in the SU may destroy the vortex-free condition, it can be avoided by keeping the degenerate singular values.Figure 2 (b) presents the energy flow as a function of the ITE step.The zeroth step denotes the LG state, and the energy drops quickly and converges already in a few hundred steps.Throughout the ITE, the flux expectation value keeps unity, W p =1, while the state remains non-magnetic, S =0 up to the machine precision [see SM [29]].The variational energy decreases monotonically with increasing D, and the lowest energy obtained by the ITE is E D=12 ITE =−0.6453 which is comparable to the ones predicted by ED and DMRG.Comparing to the energies of the SG states, the ITE en- ergies are lower at the same D.However, this is not so surprising since the tensors of the SG states are quite sparse, i.e., many of elements are zero.To estimate the infinite-D variational energy, we have extrapolated the energies E SG 's and E ITE 's, and both seem to converge to almost the same one E D→∞ −0.6464 as shown in Fig. 2 (c).Finally, to see how the gapped nature of the LG state is affected by the ITE, we have directly computed the largest correlation length ξ from the largest and and second largest eigenvalues, λ 0 and λ 1 , of the transfer matrix in the CTMRG algorithm: ξ = −1/ log(λ 1 /λ 0 ).The result from the D=10 ansatz is presented in Fig. 2 (d) as a function of the CTM dimension χ.It clearly indicates a finite and short correlation length in contrast to the case of the spin-1/2 where ξ diverges as increasing χ [28].Additionally, the smooth connection from the SG ansatz to the LG state, which is strictly gapped, is a strong evidence of the gapped nature of the S=1 KSL.
Topological property.-On a compact geometry, one should consider not only the local fluxes but also the global fluxes defined on non-contractable paths in the system.For instance, we are able to define a global flux operator on the cylinder, e.g., W Γ = i∈Γ U y i wrapping the cylinder.Then, one can sectorize the Hilbert space further using the global flux W Γ = ±1.It has been found that, in the case of the spin-1/2, the sector is characterized by the parity of the number of non-contractible loops in the configuration [35].For instance, the parity of the number of loops connecting two edges of the cylinder determines the sector.On the other hand, the S=1 LG state in each sector is characterized in a completely different way.Namely, all loop configurations are allowed regardless of the sector, whereas the sign of the configurations with the odd number of the loops winding the .This qualitative difference between the S=1/2 and S=1 LG states comes from Eq. ( 1) which is different from the S=1/2 case [see the SM [29] for more details].It can be shown that the (−)-sector can be obtained by twisting the gauge of the LG operator using the non-trivial element of the Z 2 invariant gauge group [29].Then, one can construct the minimally entangled states (MES) on the infinite cylinder (say, L x → ∞, L y : finite) using the degenerate states in the different flux sector [36].The explicit definition of the MES is presented in the SM [29].It is worth noting that one can construct all the MESs in the S=1 model in contrast to the S=1/2 model where only the trivial and vortex sectors are accessible [35].Using the boundary theory of TN [37], we have calculated the entanglement entropy (EE) and then extracted the topological entanglement entropy (TEE) to identify the nature of each anyon sector.The result is shown in Fig. 3 (a), where the EE scales as S vN = (log 2)×L y − γ, and TEE is γ=0 for |(+, +) and γ=log 2 for all other MESs.It indicates that all topological excitations are Abelian so that the ground state possesses the Z 2 topological order.Figure 3 (b) presents the φ-dependence of the TEE of |ψ SG1 .When the SG state is close to the LG state (near φ = 0), the TEE is almost perfectly log 2. As φ increases, the TEE deviates from log 2 which might come from the stronger finite size effect (L y ).
Antiferromagnetic model.-Remarkably, the relation in Eq. ( 1) allows us to construct the ansatz for the antiferromagnetic model (K<0) without much effort.First, we prepare a TN operator, say Q LG , by inserting σ xmatrix into every bond in the Q LG -operator as illustrated in Fig. 4 (a).Then, using Eq. ( 1), one can show that the Q LG -operator also guarantees the vortex-free condition, i.e., W p Q LG = Q LG , and is related to the LG where V is a unitary transformation flipping the overall sign of the Hamiltonian, i.e., V † HV = −H.Detailed derivation and proofs are provided in the SM [29].
It implies that the deformed LG state gives exactly the same variational energy as that of the original LG state for the antiferromagnetic model: Besides, since the DG operator commutes with both Q LG and Q LG , one can also deform the SG state to have the same energy for the antiferromagnetic model.That is, the variational energy of Therefore, once the variational parameters are optimized for the ferromagnetic model, one can use the same ones for the antiferromagnetic model with insertion of the σ xmatrix in the TN state.As will be shown below, this simple transformation comes in handy when we consider a weak perturbation, e.g., the magnetic field, in the antiferromagnetic model.It is also worth noting that such TN deformation for the S=1/2 model exists but requires a larger unit-cell in the TN state [29].
Magnetic field.-Here, we discuss how the S=1 KSL is affected by the (111)-direction magnetic field, To this end, we perform the ITE from the LG and the deformed LG states for the ferromagnetic and antiferromagnetic models, respectively, in the presence of the magnetic field.We have found that the response to the field resembles the ones of the S=1/2 model [38,39].That is, the ferromagnetic KSL is easily driven to the polarized phase, in which the spins are aligned in the (111)-direction, at a weak field h D=12 c ≈ 0.013 as shown in Fig. 5 (a), where the D=12 TN is optimized to determine the critical field.On the other hand, Figure 5 (b) indicates that the antiferromagnetic KSL is quite robust against the field where there is no transition at least below h = 0.3.According to the second derivative of the energy and the first derivative of the magnetization and flux expectation value [29], the transition occurring in the ferromagnetic model is expected to be continuous.It is also noteworthy that the D-dependence of the ansatz is stronger than that of the S=1/2 model [40].The D-dependence of the results is discussed in the SM [29].
Conclusion. -We have studied the ground state properties of the S=1 Kitaev model employing the LG, SG ansatz and the ITE optimization.The ground state is found to be a gapped Z 2 spin liquid which differs from the S=1/2 model hosting the gapless KSL ground state [1].The gapped nature has been shown by mapping the LG state into the partition function of the classical O (1) LG model in the gapped phase and by directly computing the largest correlation length of the optimized ansatz.By constructing the MESs, we have identified the topological nature of the excitations, i.e., the Abelian anyons.The TN representation allows us to understand the qualitative difference between the S=1 and S=1/2 LG and SG states including the topological property.In addition, we have found a simple transformation between the tensor network wavefunctions of the ferromagnetic and antiferromagnetic models, leading to exactly the same variational energy.Therefore, once one optimizes the variational parameters with the ferromagnetic model, the ansatz for the antiferromagnetic model is obtained without additional optimization (or vice versa).The magnetic field effect on the KSL is similar to that of the S=1/2 model.That is, the antiferromagnetic KSL is much more robust than the ferromagnetic KSL against the field.We believe that our construction is applicable to the general spin-S Kitaev model.

Symmetries
At the isotropic point (K x = K y = K z ), the Hamiltonian is invariant under the (C 6 V C6 )-, (σV σ )-and timereversal (TR) transformations, where C 6 and σ denote the C 6 lattice counterclockwise rotation and inversion along a line aligned in the z-bond, respectively.And, the unitary operators V C6 = exp i 4π 3 √ 3 (S x + S y + S z ) and V σ = exp i π 2 S z transform the spin operators as follows: Due to the bond-dependent interaction in the Kitaev model, the transformation V C6 rotates the bond clockwise so that the combined transformation with C 6 , or (C 6 V C6 ), leaves the Hamiltonian invariant.Similarly, it is symmetric under the (σV σ )-transformation.
One can show that the loop gas and dimer gas operators are also invariant under those transformations.Note that V C6 and V σ transform the U γ -operator in a similar fashion: Therefore, the Q-tensor is transformed as follows: Note that the C 6 -rotation and σ-inversion permute the virtual indices of the tensor in the following way: Therefore, the combined transformations (C 6 V C6 ) and (σV σ ) leave the Q-tensor invariant, and thus the resulting loop gas operator, too.Since the TR-transformation acts trivially on U γ , i.e., T −1 U γ T = U γ , it leaves the loop gas operator intact.In a similar way, one can show that the dimer gas operator is symmetric under those transformations.

TOPOLOGICAL PROPERTY
Global flux on the loop gas (LG) and string gas (SG) states On a compact geometry, the Hamiltonian commutes with the global flux operator W P ≡ i∈P U γi where P is a non-contractible path winding the system.The eigenvalue ±1 of such global flux operator can be used to characterize the topological sector.On the torus geometry, we are able to define two global flux operators, say W h = i∈h U z and W v = i∈v U y wrapping the torus around the horizontal and vertical directions, respectively.Then, one can show that the LG state (and thus SG state) is the eigenstate of W h and W v with eigenvalue +1.To be more precise, let us consider the (2 × 2) system with the periodic boundary condition and then apply both operators on the LG state, which can be illustrated as , (10) where the z (y) on the circle stands for the multiplication of the U z (y) operator on the corresponding site.In the first equality of each equation, the equation 2 is used, i.e. the multiplication of U z (y) is replaced by the multiplication of σ x (green square) on the virtual bond.Due to (σ x ) 2 = 1, all the green squares in the second expression in each equation are cancelled, and thus W h |LG = +1|LG = W v |LG .It indicates that our LG state on the torus geometry is already the eigenstate of both global flux operators in contrary to the S = 1/2 LG state [35].Due to this difference in the S = 1 and S = 1/2 LG states, the construction of the minimally entangled state (MES) also differs from each other.

Access to other global flux sector
As shown above, the LG state obtained by the LG operator is in the (W h , W v ) = (+1, +1) sector.Note that one can access other sectors by inserting the row and column of the non-trivial element (g = σ z ) of IGG around the system, say G h and G v .Then, it changes the eigenvalue of the global flux operator.For instance, in the (2 × 2)-system, let us insert it along the vertical direction (G v ) and then apply the W h operator as depicted below: . (11) Similar to Eq. ( 10), the multiplication of U z can be replaced by multiplication of σ x (green square) on two virtual legs, and two successive green squares are cancelled.However, note that σ x σ z σ x = −σ z appears when the σ z is inserted.Therefore, inserting G v flips the horizontal global flux number W h .In a similar fashion, one can construct the LG states living in each sector (W h , W v ) = (±, ±): (12) which is very different from the S = 1/2 case where, for instance, the (+, +)-sector is obtained by summing over all of the states above [35]: (13) In the case of the n-th order SG state, the non-trivial element of IGG becomes g = σ z ⊗ n i=1 I 2 , and the construction of each sector is exactly the same as above.

Minimally entangled states on the infinite cylinder
In the presence of the Z 2 gauge structure, one can set the MES on the infinite cylinder (L x → ∞, L y : finite) as follows [36]: Since any flux sector can be accessed simply by inserting G h and G v in the original tensor network, we can construct all the minimally entangled states.

TENSOR NETWORK DEFORMATION FOR ANTI-FERROMAGNETIC MODEL
As pointed out in the main text, the deformed LG state can be used for the antiferromagnetic model.First, using Eqs.( 1) and ( 2), one can show that the deformed LG operator depicted in Fig. 1 (a) guarantees the vortex-freeness as follows.In an analogous manner to Eq. ( 4), applying the flux operator to the deformed LG operator, one obtains a similar equation

Imaginary time evolution flow of physical quantities
The imaginary time evolution (ITE) starting from the loop gas state is stable and effective in the presence of weak magnetic field.As illustrative examples, we present the ITE flow of the variational energy, magnetization and flux expectation value as a function of the ITE step in Fig. 4 at h = 0 and h = 0.02.As one can see, the magnetization and flux expectation value remain unity and zero, respectively, throughout the ITE flow at zero field.On the other hand, at h = 0.02, the ITE flow exhibits transition-like behavior at a certain ITE step, at which the ansatz seems becoming the polarized state.

Different flux sector
The LG representation allows us to access the other flux sector.For instance, an ansatz in the vortex-full sector, in which all plaquetts are occupied by the vortex W p = −1, can be obtained by inserting the non-trivial element of the invariant gauge group as illustrated in Fig. 5 (a).Then, one can perform the ITE in a given flux sector as presented in Fig. 5 (b) and (c) for the vortex-free and vortex-full sectors.The ansatz in the vortex-free sector has significantly lower energy than that of the vortex-full state.We have also checked the vortex-crystal sectors and found that all those sectors are unstable compared to the vortex-free one.As demonstrated in the main text, the ferromagnetic model exhibits a phase transition in the presence of the (111)-direction magnetic field.Figure 6 presents the bond dimension D-dependence of the energy, magnetization and flux expectation value as a function of the field strength h.Interestingly, when D is small (D < 10), the transition is discontinuous, i.e., the first order transition, while it turns out to be the continuous one with D ≥ 10.However, the critical field is still not converged up to D = 12, so that we expect the true critical field is slightly larger than h D=12 c ≈ 0.13.We also present the energy and its second derivative, and the flux expectation value and its first derivative in Fig. 7.

FIG. 1 .
FIG. 1.(a) Graphical representation of the tensor network state on the honeycomb lattice where the x, y and z bonds of the model in Eq. (1) are denoted by red, blue and yellow colors, respectively.(b) The local states depending on the local loop configuration appearing in the LG state [Eq.(3)].Here, |0 stands for the (111)-direction magnetic state, i.e. 0|S γ |0 =1/ √ 3, and U γ = e iπS γ is the SO(3) rotation operator.
This is, the equal weight superposition of all possible loop configurations where the local state depends on the direction of the local loop as depicted in Fig. 1 (b).Using the identity 0|U γ |0 =1/3, one can map the norm of |ψ LG into the partition function of the classical O(1) where the nonzero elements of ζ-tensor are simply ζ φ 000 = cos φ and ζ φ 100 = ζ φ 010 = ζ φ 001 = sin φ.Applying the DG operators onto the LG state results in the string gas (SG) state with variational parameters {φ α }, i.e., |ψ SGn ({φ α }) ≡ n α=1 R DG (φ α )|ψ

D = 4 :
FIG. 2. Variational energies of (a) the first order SG state |ψSG 1 (φ) and (b) the ITE from the initial LG state at a given bond dimension D. (c) D-interpolation of the variational energies obtained from both SG ansatz and ITE method.(d) The largest correlation length ξ of the D=10 ITE ansatz as a function of the CTM dimension χ.

1 FIG. 3 .
FIG. 3. (a) The entanglement entropy SvN of the LG state as a function of the circumference Ly in the (+, +)-flux and the MES sectors.(b) The extracted topological entanglement entropy of the first order SG state |ψSG 1 (φ) .The definition and construction of four MESs, |I , |e , |m and | are presented in the SM [29].

FIG. 4 .
FIG. 4. Schematic figure of the TN state for the antiferromagnetic model, where the green square denotes the Pauli σ x -matrix, and the original TN state (without green squares) is either the LG state or the SG state.

FIG. 5 .
FIG. 5.The magnetizations of (a) the ferromagnetic model and (b) the antiferromagnetic model as a function of the field strength h.The results are obtained by the ITE with D = 12.

FIG. 4 .
FIG. 4. Plots of the imaginary time evolution flow of (a), (c) the variational energies and (b), (d) the magnetization and flux expectation value at h = 0 and h = 0.02, respectively.The initial state is the loop gas state, and the bond dimension is D = 8 and the imaginary time step is τ = 0.01.

FIG. 5 .FIG. 6 .
FIG. 5. (a) Schematic figure of the vortex-full TN state where the yellow square stands for the non-trivial element of Z2 invariant gauge group [see Sec.].(b) The ITE flow of the variational energies in the vortex-free and vortex-full sectors with D = 8 at zero field.

FIG. 7 .
FIG. 7. Plots of (a) the variational energy and its second derivative and (b) the flux expectation value and its first derivative as a function of the field strength h.Here, the D = 10 ansatze are used.