Path to stable quantum spin liquids in spin-orbit coupled correlated materials

The spin liquid phase is one of the prominent strongly interacting topological phases of matter whose unambiguous confirmation is yet to be reached despite intensive experimental efforts on numerous candidate materials. Recently, a new family of correlated honeycomb materials, in which strong spin-orbit coupling allows for various bond-dependent spin interactions, have been promising candidates to realize the Kitaev spin liquid. Here we study a model with bond-dependent spin interactions and show numerical evidence for the existence of an extended quantum spin liquid region, which is possibly connected to the Kitaev spin liquid state. These results are used to provide an explanation of the scattering continuum seen in neutron scattering on $\alpha$-RuCl$_3$.


Introduction
The role of strong interaction between electrons in the emergence of topological phases of matter is currently a topic of intensive research. The archetypal example of a topological phase with strong electron interaction is the quantum spin liquid 1 , in which the elementary excitations are charge-neutral fractionalized particles. While a lot of progress has been made on the theoretical understanding of the quantum spin liquid phase, its direct experimental confirmation has remained elusive despite various studies on a number of candidate materials 2-6 . Significant progress, however, has recently been made due to the availabilty of a new class of correlated materials, where strong spin-orbit coupling leads to various bonddependent spin interactions [7][8][9] , thus resulting in magnetic frustation. These materials are Mott insulators with 4d and 5d transition metal elements, which include iridates and ruthenates [10][11][12][13] and come in two-or three-dimensional honeycomb variants. They have been particularly exciting because they intrinsically have a strong Kitaev interaction and therefore could potentially realize the Kitaev spin liquid (KSL) phase -an example of a Z 2 quantum spin liquid where the electron's spin- 1 2 fractionalizes into two degrees of freedom: itinerant Majorana fermions and Z 2 fluxes.
While the Kitaev interaction (K) in these materials is large, it competes with symmetry allowed nearest-neighbor (n.n.) symmetric off-diagonal (Γ) and Heisenberg (J) spin interactions 14 . For example, in α-RuCl 3 (RuCl 3 ), an actively studied KSL candidate, comprehensive ab initio computations and recent dynamical studies 15,16 suggest that ferromagnetic K and antiferromagnetic Γ interactions are dominant and comparable in magnitude, while J is negligible [17][18][19] . The balance of these and additional small further neighbor interactions causes RuCl 3 and other KSL candidates to magnetically order at low temperature; however, it is still unclear whether or not the often-large Γ interaction prefers magnetic ordering. Meanwhile, the community has attempted to revive the possibility of a KSL in RuCl 3 by applying a small magnetic field, with the effect of entering a potential spin liquid phase with no magnetization 19 .
Since a weak magnetic field takes RuCl 3 out of the ordered phase, it lends credence to the idea that the zig- zag phase is stabilized by small interactions at comparable energy scale to the magnetic field, such as a 3rd n.n. Heisenberg J 3 17-19 term or terms coming from slight trigonal distortion 20 . This calls into question the role of the Γ interaction. Interestingly, a recent analysis of the Γ model revealed a macroscopically degenerate classical ground state 21 .
In this work we will thus investigate if a model with K and Γ hosts an extended quantum spin liquid phase. A previous exact diagonalization study on a 24 site honeycomb cluster hints that the ferromagnetic KSL is unstable after perturbing with a small Γ, but the resulting phase is not orderered 14 . On the other hand, it is known that the KSL is stable upon introducing bond anisotropy, which is present in real materials as depicted in Fig. 1a. We indeed find that such anisotropy can extend the KSL phase between the −K and Γ limits, as shown in Fig. 1b.
We consider the following nearest-neighbor (n.n.) model and H x,y are defined similarly with corresponding K x,y and Γ x,y . Each H γ represents the n.n spin interactions along one of the three bond directions, γ = x, y, z. The model is parameterized by K z = −(1 + 2a K ) cos φ, K x,y = −(1−a K ) cos φ, Γ x,y,z = sin φ, with a K characterizing bond anisotropy. When φ = 0, π (i.e. Γ γ = 0), this model reduces to the exactly solvable Kitaev model with the KSL ground state.
We have studied this model using a combination of three different, corroborating, numerical methods: exact diagonalization (ED) on a 24-site honeycomb cluster, the method of thermal pure quantum states 4-8 , and infinite time-evolution block decimation (iTEBD). We have first reproduced the earlier work in the isotropic a K = 0 limit, showing a strong first-order transition between −K and Γ limits (see Supplementary Materials (SM)). We present the following results using our numerical techniques: 1) When a K ≥ 0.06, we find that the −K γ (0 ≤ φ ≤ π/2) and Γ γ (φ/π = 0.5) limits are adiabatically connected as shown in Fig. 1b. Thus we find evidence for an extended quantum spin liquid phase in the presence of anisotropy, a K . 2) An intervening magnetically ordered phase separates the spin liquid phase near the pure Γ γ limit and the antiferromagnetic KSL at φ/π = 1.
3) The specific heat C(T ) and entropy S(T ) at finite temperatures suggest a smooth crossover from the ferromagnetic Kitaev limit to the pure Γ γ limit, consistent with our ED results. 4) Zig-zag spin correlations become dominant upon perturbing the quantum spin liquid phase in 0 < φ < π/2 by J 3 , indicating the enhancement of zig-zag order by J 3 .

Results
Extended spin liquid state in global phase diagram -The ground state energy per site E 0 /N of Eq. (1) was computed for φ/π ∈ [0, 1], and for different anisotropy parameters by ED on a 24-site cluster using periodic boundary conditions (see SM). Discontinuities in 1 N ∂E0 ∂φ were used to identify possible phase transitions. Remarkably, when φ/π ∈ [0, 0.5] and 0 ≤ a K < 0.06, there is a line of first order phase transitions that terminate at a K = 0.06. Above a K = 0.06, the first derivative of the energy presents no sharp features suggesting that the ground state of the Γ-limit (φ/π = 0.5) is adiabatically connected to the ferromagnetic Kitaev spin liquid (φ/π = 0) as depicted in Fig. 1b for a K = 0.1.
In the antiferromagnetic region of phase space, there are two large discontinuities in 1 N ∂E0 ∂φ that encompass a large region of phase space separating the Γ-limit and the exactly solvable antiferromagnetic Kitaev limit at φ/π = 1. These peaks coincide with kinks in E 0 /N (solid yellow) shown in Fig. 2a. Two smaller discontinuities can also be seen near φ/π = 0.75, however these are not present when a K = 0, while the larger jumps near φ/π = 0.5 and 1 appear consistently for different a K . The small discontinuities can thus be considered spurious and a consequence of the finite cluster size. Similar finite size effects were also found for φ/π ∈ [0, 0.5] when a K = 0, as discussed in the SM.
Magnetic order and perturbations -The ground state wavefunction of Eq. (1) computed by ED is used to evaluate real-space spin-spin correlation functions S i · S j , where i and j are site indices on the honeycomb lattice. By Fourier transform, we obtain the static structure factor (SSF) given by S q = 1 N i,j e i(ri−rj )·q S i · S j where q is a vector in the reciprocal lattice. The SSF at various points in the Brillouin zone (BZ) is plotted over the phase space in the bottom panel of Fig. 2a.
The discontinuities in the SSF can be directly matched with those in 1 N ∂E0 ∂φ . Visualizations of the SSF over the BZ for representative φ in the phase diagram are presented in Fig. 2b. The SSF in Fig. 2b is obtained by computing the average of S i · S j over all n.n. bonds, 2nd n.n., etc. This calculation reflects the presence of different domains in the crystal, in which either of x, y, z bond interactions can be stronger and thus, over the whole crystal, these domains result in an isotropic S q despite the inherent bond anistropy in Eq. (1). The SSF varies adiabatically when a K = 0.1 for φ ∈ [0, π/2] and the spin correlations at the Γ-and Mpoints are comparable in intensity when Γ K γ , leading to a "star"-shaped structure in the SSF as seen in Fig. 2b (e.g, φ/π = 0.2). The extended phase separating φ/π = 0.5 and 1 is characterized by dominating spin correlations at the Kand Γ -points in the reciprocal lattice (φ/π = 0.75 in Fig.  2b). Contained within this phase is the exactly solvable point with hidden SU (2) symmetry at φ/π = 0.75 which features K-point correlations 28 consistent with the results presented here. Thus the extended spin liquid phase for ferromagnetic K is separated from the antiferromagnetic KSL at φ/π = 1 by a magnetically ordered phase. These results can be connected to real materials, particularly RuCl 3 in which a zig-zag magnetic ordering has been observed [29][30][31] . Previous studies have shown that in addition to the n.n. ferromagnetic Kitaev and antiferromagnetic Γ interactions, a 3rd n.n. antiferromagnetic Heisenberg interaction J 3 i,j S i · S j is non-vanishing and plays a role in determining the magnetic ordering in RuCl 3 17,18 . Fig. 3 shows that the effect of perturbing Eq. (1) by J 3 is to enhance (suppress) the M-point (Γ-point) spin correlations, consistent with a zig-zag magnetically ordered state observed in experiments. This result indicates that by tuning J 3 in the real material, an alternate path to achieve a spin liquid phase may be realized.
Specific heat and thermal entropy -Previous study on the finite temperature properties of the Kitaev model has shown that Kitaev spin liquids feature two peaks in the heat capacity C(T ) and a 1 2 -plateau in the entropy S(T ), which is attributed to the thermal fractionalization of spin degrees of freedom 32 . Here we go beyond the Kitaev limit and investigate the heat capacity and entropy at finite temperature in the presence of Γ, which is expected to compete with K γ in RuCl 3 , using the method of thermal pure quantum states (see SM).
The dependence of C(T ) on φ when a K = 0.1 is plotted in the top panel of Fig. 4. The expected two peak structure in C(T ) is observed when φ/π = 0, and is seen to be maintained continuously as φ/π approaches 0.5 so that the Γlimit shows a qualitatively similar behaviour in C(T ) to the Kitaev spin liquid. Evidence for a phase transition can be seen when φ 0.7 on account of the abrupt change in C(T ), resembling that of the heat capacity in trivially ordered phases 33 . This finding is consistent with our ED results and the 120-order at φ/π = 0.75 seen in Refs. 14 and 28. The dependence of S(T ) on φ is plotted in the bottom panel of Fig. 4 with a clear 1 2 -plateau observed when φ/π = 0, consistent with the expected Kitaev spin liquid behaviour. In addition, a plateau of about 1/5 the total entropy is observed when φ/π = 0.5. Another plateau is observed in the magnetically ordered phase around φ/π = 0.75; however, this feature can be attributed to finite-size effects as follows. The (N + 1)-fold ground state degeneracy at φ/π = 0.75 due to the hidden SU (2) symmetry 28 is only slightly broken away from this point, inducing a plateau in S(T ) with height given by ln(N + 1)/N ln 2 0.1935 ∼ 1/5 when N = 24. By contrast, the height of the plateau around φ/π = 0.5 is independent of N 34 .
The physical origin of the two peak structure in C(T ) and the plateau in S(T ) can be traced to the energy scales of the thermal fluctuations of the underlying quasiparticles in the spin liquid 32 . In the Kitaev spin liquid at zero temperature, the low-lying quasiparticle excitations are characterized by itinerant Majorana fermions which disperse in a background of zero flux 35,36 . It has been shown that as temperature increases, the flux degrees of freedom begin to fluctuate and lead to the low temperature peak seen in C(T ), resulting in the plateau seen in S(T ). Furthermore, the high temperature peak in C(T ) is attributed to the development of short range spin correlations 32 . Our results show that the two-peak structure in C(T ) is qualitatively maintained and further suggests that no phase transition has taken place.
Similarities on the infinite tree -We further studied Eq. (1) on an infinite Cayley tree with z = 3 connectivity, using the infinite time-evolving block decimation algorithm 9 (iTEBD; see SM). Classically, the ground state in the Γlimit on the infinite tree is macroscopically degenerate because a different state with the same energy can be constructed by flipping the sign of one spin component on an infinite string of neighboring spins. The Γ-limit on the two-dimensional (2D) honeycomb and three-dimensional (3D) hyper-honeycomb 21 lattices also feature similar classical degeneracy. The similarity at the classical level of the Γ-limit on the infinite tree to the 2D and 3D lattices prompts us to study the quantum model on the infinite tree for further insight. Figure 5 shows results of the eight-site iTEBD calculation with bond dimension χ = 10, and anisotropy a K = 0.1. In this calculation, we have also introduced an anisotropy to Γ γ such that Γ x = Γ y = (1 − a Γ ) sin φ and Γ z = (1 + 2a Γ ) sin φ in order to apply the iTEBD method (see SM). No transition is found when φ/π ∈ [0, 0.5] and the obtained state is a highly entangled paramagnet, with S E ∼ 0.8 for strong (z) bonds, while for weak (x, y) bonds, S E ∼ 0.4. Deep in the gapped phase of the Kitaev model, with large anisotropy a K , one finds S E ∼ log 2 ∼ 0.693 for the strong bonds and much smaller values of S E for the weak bonds. Both, however, increase as the anisotropy is reduced -perhaps due to a finite contributions from the Majorana fermions 13 . An increase is expected in the entanglement entropy as one approaches a phase transition, however no such peaks are seen for 0 < φ < π/2. Similarly, there are no sharp features in the ground state energy E 0 as a function of φ, which indicates that this phase is adiabatically connected to the Kitaev spin liquid at φ = 0.
There is an apparent first order transition around φ = 0.6π into a Néel state with spins ordered in the [111] direction, accompanied by a dramatic lowering of S E on both strong and weak bonds into this region. This Néel state becomes a simple product state when φ/π = 0.75, as seen by the vanishing of S E . A final transition into a paramagnetic state is seen before the antiferromagnetic Kitaev limit.

Discussion
The highlight of our numerical results is that, in the presence of bond anisotropy a K , there exists an extended quantum spin liquid region which is adiabatically connected to the ferromagnetic KSL. The model we have studied is motivated by experiments on RuCl 3 and earlier ab initio computations [17][18][19] . In a recent inelastic neutron scattering experiment on RuCl 3 , it is found that the continuum of finite energy excitations exists both below and above the magnetic transition temperature despite that the low temperature ground state is the zig-zag long-range ordered state 39 . The inelastic neutron scattering data for the continuum show the star-shape intensity that extends from the zone center towards the M points of the Brillouin zone.
Recall that the static structure factor in our ED study shows enhanced (decreased) short-range spin correlations at the M point (zone center) of the Brillouin zone as one moves from the ferromagnetic Kitaev limit to the pure Γ γ limit. When the strength of the ferromagnetic Kitaev interaction and the Γ γ interaction become comparable, both of the short-range spin correlations at the M and the zone center would show significant intensity, which leads to the star-shape structure in momentum space. This behavior may be favorably compared to the finite-energy short-range spin correlations seen in RuCl 3 . Given that the ab initio computations suggest comparable magnitudes of the ferromagnetic Kitaev and Γ γ interactions in RuCl 3 17 , it is conceivable that RuCl 3 may be very close to the quantum spin liquid phase found in our model and, as shown in our work, the introduction of small J 3 would favor the zig-zag magnetically ordered phase observed in RuCl 3 .
Finally, more analytical understanding of the connection between the pure Kitaev limit and the quantum spin liquid phases identified in our numerical work would be extremely valuable for future applications on real materials. Note that a possible incommensurate magnetic ordering cannot be ruled out due to finite cluster size. However, based on the results of a recent iDMRG study 22 , no evidence of incommensurate ordering allowed by their momentum cuts was found. It is also interesting to note that quantum fluctuations do not lift the infinite ground state degeneracy of the classical model for positive Γ, while they may lead to incommensurate ordering for negative Γ 21 . Thus it is likely that the positive Γ regime studied here possesses a spin liquid ground state, and the precise nature of the spin liquid is an excellent topic for future study.

Methods
Our results were obtained using the combination of the three independent numerical techniques listed below.
Exact diagonalization -ED was performed on a 24 site cluster with periodic boundary conditions. This cluster allows all the symmetries present in the infinite honeycomb lattice and has been used reliably in previous related classical and quantum studies 14,40 . The Hamiltonian given by Eq. (1) in the main text does not have the U (1) symmetry associated with S z conservation, making it impossible to block diagonalize into magnetization sectors. Therefore, the translational symmetry of the 24 site cluster was used to block diagonalize into different momentum sectors to gain more information about its energy spectrum. The lowest energies and corresponding wavefunctions of each block were then numerically obtained using the Lanczos method. Further details and calculations can be found in the SM.
Thermal pure quantum states -We used the method of thermal pure quantum states 7,8 in our specific heat and thermal entropy calculations. Details about the construction of thermal pure quantum states and the subsequent calculation of specific heat and entropy can be found in the SM.
Infinite time-evolving block decimation algorithm -The Hamiltonian given by Eq. (1) was studied on an infinite Cayley tree with z = 3 connectivity using the infinite time-evolving block decimation algorithm 9 . Details about the method and the construction of the ground state can be found in the SM.

Acknowledgements
This work was supported by the NSERC of Canada and the Center for Quantum Materials at the University of Toronto. Y. Y. was supported by JSPS KAKENHI (Grant Nos. 15K17702 and 16H06345) and was supported by PRESTO, JST. Y. Y. was also supported in part by MEXT as a social and scientific priority issue (Creation of new functional devices and high-performance materials to support next-generation industries) to be tackled by using post-K computer. Computations were mainly performed on the GPC supercomputer at the SciNet HPC Consortium. SciNet is funded by: the Canada Foundation for Innovation under the auspices of Compute Canada; the Government of Ontario; Ontario Research Fund -Research Excellence; and the University of Toronto. A part of the TPQ results were checked by a program package, HΦ 41 . We thank helpful discussions with Frank Pollmann, Matthias Gohlke, Shunsuke Furukawa, and Subhro Bhattacharjee. We particularly thank Natalia Perkins and Ioannis Rousochatzakis for informing us of their unpublished ED results on related models.

Competing Interests
The authors declare no competing financial or nonfinancial interests.

Data Availability
All relevant data is available from the corresponding author.

I. EXACT DIAGONALIZATION
Ground state energy and Sq in the aK = 0 limit When the Kitaev exchange is isotropic a K = 0, three discontinuities in 1 N ∂E0 ∂φ , which are absent with slight anisotropy a K = 0.1, can be seen when φ/π ∈ [0, 0.5]. These discontinuities are reflected in S q as can be seen in Fig. 6. The general qualitative picture, however, is the same as with slight anistropy: comparable Γ-and Mpoint spin correlations when φ/π ∈ [0, 0.5] lead to a "star"shaped pattern in S q in the Brillouin zone (BZ). As when a K = 0.1, two large discontinuities in 1 N ∂E0 ∂φ and a dramatic change in S q when K γ > 0 are evidence of an extended magnetically ordered phase with dominating K-and Γ -point correlations. The inconsistent appearance of the other discontinuities in 1 N ∂E0 ∂φ (and associated small discontinuities in S q ) with changing anisotropy parameter are likely finite-size effect manifestations. Indeed, the location of these small discontinuities depends on system size (N = 18, 24). We investigated the phase diagram between φ/π = 0 and φ/π = 0.5 for different a K as shown in Fig. 8 calculated by ED. The discontinuities seen in ∂E0 ∂φ when a K = 0 disappear for a K ≥ 0.06. The smooth connection between φ/π = 0 and φ/π = 0.5 limits is found to be robust for a wide range of a K . Black vertical bars label the location of small spurious discontinuities in ∂E 0 ∂φ which appear specifically when N = 24. These transitions do not appear when N = 18 nor in recent iDMRG calculations 1 . We conclude that these transitions are likely due to finite-size effects, however, they do not affect our main conclusion, which is that the KSL and Γ limits are adiabatically connected for small aK .

Finite size spectrum
The Hamiltonian was block diagonalized in terms of eigenvalues of the translation operators T 1,2 acting on the spin configuration on the 24 site cluster. In this notation, any spin configuration in the S z basis, |s ≡ |S z 1 . . . S z 24 , on the 24 site cluster can be translated as T n 1 T m 2 |s where n = 0, ...5 and m = 0, 1. Periodic boundary conditions are imposed such that T 6 1 = 1 and T 2 2 = 1. Thus there are 12 unique eigenvalues of T n 1 T m 2 (there are two atoms in the unit cell) of the form e 2πi(kx/N1+ky/N2) , with N 1 = 6 and N 2 = 2 on the 24-site cluster. The momentum eigenvalues are labeled by the integer values k x and k y . For all φ, the ground state is in the q ≡ (k x , k y ) = 0 momentum sector. The energy spectra corresponding to φ/π = 0.5 and 0 are shown in Fig. 9a. In the Γ-only limit, the ground state is accompanied close in energy by 3 excited states in the q = 0 sector, the first excited state being doubly degenerate. These states are separated by a gap from the rest of the spectrum. The Kitaev limit hosts a doubly degenerate ground state within ED.
It is found that when varying φ and keeping the Kitaev term isotropic, the four states in the q = 0 sector experience level crossing, shown in Fig. 9b, and account for the firstorder phase transitions found when φ ∈ [0, π/2] seen in Fig. 6. The location of these level crossings when φ ∈ [0, π/2] depends on system size (N = 18, 24) and are thus attributed to finite size effect. By contrast, there is no level crossing when introducing small anisotropy a K = 0.1 to the Kitaev exchange. The energy spectrum surrounding the exactly solvable point φ/π = 0.75 tends to be degenerate and is qualitatively different from the spectrum in the rest of the phase space. This result serves as a way to contrast the energy spectrum of the magnetically ordered phase with the rest of the phase diagram.

J3 perturbation
Under the effect of a 3rd n.n. J 3 > 0, the intensity in S q at the M-points in the BZ was shown in the main text to be enhanced. This result suggests that a phase transition into the experimentally observed zig-zag magnetically ordered state could occur for sufficiently large J 3 > 0. Fig. 10 shows signatures of a possible phase transition when investigating the second derivative of the energy (dashed lines) and fidelity susceptibility 2 (solid lines) χ F when φ/π = 0.5, 0.1 and 0.2. The broad peaks in both χ F and − 1 N ∂ 2 E0 ∂φ 2 could be a result of finite size effect, and makes it difficult to estimate the necessary J 3 that would induce a transition into the zig-zag ordered phase.

II. THERMAL PURE QUANTUM STATES
To capture thermal excitations of Eq. (1) in the main text, and examine how the Kitaev spin liquid is deformed in the presence of finite Γ, we systematically calculate heat capacity by using typical pure states 3 . Several pioneering studies 4-6 have shown that averages over a few tens of pure states can replace a canonical ensemble. Recently, a statistical mechanics proof for the replacement of a canonical ensemble by a single pure state was given in Refs.7 and 8, where such a pure state that replaces a canonical ensemble is called a thermal pure quantum (TPQ) state 7,8 .
We briefly summarize the construction of TPQ states following Ref. 7. A TPQ state of N quantum S = 1/2 spins at infinite temperatures is simply given by a random vector, where |i is represented by the real-space S = 1/2 basis and specified by a binary representation of a decimal number and {c i } is a set of random complex numbers with the normalization condition 2 N −1 i=0 |c i | 2 = 1. Then, by multiplying a wave vector by the hamiltonianĤ, we construct the TPQ states at lower temperatures as follows: Starting with an initial (or the 0-th step) vector |Φ 0 = |φ +∞ , the k-th step vector |Φ k (k ≥ 1) is constructed as where Λ is a constant larger than the largest eigenvalue of H. The above k-th step vector is a TPQ state at a finite temperature T . The corresponding inverse temperature β = (k B T ) −1 is determined through the following formula 7 , where k B is the Boltzmann constant. In other word, a TPQ state at temperatures T is given as, The heat capacity and entropy ofĤ are then estimated by using a TPQ state |φ T . In thermodynamic limit, a single TPQ state |φ T indeed replaces a canonical ensemble and gives us exact heat capacity and entropy. However, for finite N , contribution from atypical states is not negligible. Therefore, in the present study, we estimate errors that originate from the atypical states by using standard deviation of the results obtained from several initial random vectors. Here, we calculate the heat capacity C by using fluctuations of internal energy. Ideally, the following equation holds: The entropy S is estimated by integrating C/T from high temperatures as

III. iTEBD METHOD
We also studied Eq. (1) in the main text on an infinite Cayley tree with z = 3 connectivity, using the infinite time-evolving block decimation algorithm 9 (iTEBD). Within this approach, the ground state is approximated by a tensor product state (TPS) defined on an infinite tree lattice. A Cayley tree TPS |ψ may be defined by site-tensors T (r),sr irjrkr , and diagonal bond-matrices λ rr iri r . Here r, r denote sites, s r =↑, ↓ is a site's spin index, and i r , j r , k r are its bond indices, for the x, y, z bonds, respectively, each running up to the bond dimension χ. The TPS is formally given by |ψ = r irjrkrsr T (r),sr irjrkr × rr ∈x λ rr iri r rr ∈y λ rr jrj r rr ∈z λ rr krk r |s r (9) In the TPS canonical form, the diagonal values of λ rr iri r determine the state's Schmidt decomposition at the corresponding bond. Therefore, the entanglement entropy for partitioning the system at that bond is simply given by S E = − i λ 2 ii log λ 2 ii . Assuming the ground sate is periodic, the Hilbert space can be greatly reduced by considering only a subset of TPSs defined by a periodic arrangement of a small number of independent tensors. In our calculation we consider an eight-site periodicity which is compatible with the simplest magnetically ordered states expected in this system. Under such an assumption, there are eight independent T site-tensors and twelve independent λ bond-matrices.
The ground state is obtained by projection, iteratively evolving the TPS in imaginary time until convergence is reached. Technically, the time evolution operator U = e −∆τ H , where H is the Hamiltonian, and ∆τ is a small time step, is applied using a fourth-order Suzuki-Trotter decomposition 10 , separately evolving each (independent) bond at a time. Each time the evolution operator is applied, the bond dimension (or Schmidt rank) χ generally increases, making the algorithm unmanageable. To avoid this, a fixed χ is maintained by applying a truncated singular value decomposition (SVD) to the evolved TPS. Thus, larger values of χ yield a more accurate ground state. Additional pure SVD steps, applied between time evolution operations, are needed in order to bring the TPS closer to the canonical form 11 . We begin by iteratively evolving a random TPS with a time step of ∆τ = 0.1 until the entanglement entropies S E on all bonds have converged. We then reduce ∆τ by a half and continue evolving, again until convergence is achieved. This process is repeated until ∆τ < 10 −5 . Since this method is defined directly in the thermodynamic limit, there are no finite size effects, and one may arrive at ground states which spontaneously break underlying symmetries.
Generally, TPSs are well suited to describe gapped phases, but not gapless ones. In the Kitaev limit, it is possible to have a gapped phase by introducing bond anisotropy 12 . When applying the iTEBD method to Eq.
(1), we similarly found that the calculation was unstable in the isotropic limit, and a small bond anisotropy was required to stabilize it. By calculating the zero energy density of states for Majorana fermions hopping on the anisotropic infinite tree, Kimchi et al. 13 were able to determine that in the K limit there is a gapped phase for a K 0.05.