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

TL;DR

This paper investigates a model with bond-dependent spin interactions in honeycomb materials, providing numerical evidence for an extended quantum spin liquid phase that could explain experimental neutron scattering results in $
$-RuCl$_3$.

Contribution

It offers new numerical evidence for an extended quantum spin liquid region in a model relevant to spin-orbit coupled materials, connecting theory with experimental observations.

Findings

Numerical evidence for an extended quantum spin liquid phase.

Possible connection to the Kitaev spin liquid state.

Explanation of neutron scattering continuum in $
$-RuCl$_3$.

Abstract

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$.

Figures

Equations

$$H=\sum _{\gamma \in x,y,z}{H}^{\gamma },$$

$$H^z=\sum _{⟨ij⟩\in z-bond}[K_z{S}_i^z{S}_j^z+{\Gamma }_z(S_i^x{S}_j^y+S_i^y{S}_j^x)]$$

$$|{\varphi }_{+\infty }⟩=\sum _{i=0}^{2^N-1}{c}_i|i⟩,$$

$$|{\Phi }_k⟩=\frac{(\Lambda -\hat{H})|{\Phi }_{k-1}⟩}{\sqrt{⟨{\Phi }_{k-1}|(\Lambda -\hat{H}{)}^2|{\Phi }_{k-1}⟩}},$$

$$\beta =\frac{2k_{B}k}{\Lambda -⟨{\Phi }_k|\hat{H}|{\Phi }_k⟩}+O(1/N),$$

$$|{\varphi }_T⟩=|{\Phi }_k⟩.$$

$$C(T)=k_{\mathrm{B}}\frac{⟨{\varphi }_T|{\hat{H}}^2|{\varphi }_T⟩-⟨{\varphi }_T|\hat{H}|{\varphi }_T{⟩}^2}{(k_{\mathrm{B}}T{)}^2}.$$

$$S=N{k}_B\ln 2-{\int }_T^{+\infty }d{T}^{\prime }\frac{C}{T^{\prime }}.$$

$$|\psi ⟩$$

Full text

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 liquidBalents (2010), 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 materialsShimizu et al. (2003); Helton et al. (2007); Okamoto et al. (2007); Yamashita et al. (2010); Han et al. (2012). 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 bond-dependent spin interactionsJackeli and Khaliullin (2009); Witczak-Krempa et al. (2013); Rau et al. (2016), thus resulting in magnetic frustation. These materials are Mott insulators with 4$d$ and 5$d$ transition metal elements, which include iridates and ruthenatesSingh et al. (2012); Plumb et al. (2014); Modic et al. (2014); Takayama et al. (2015) 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 $\mathbb{Z}_{2}$ quantum spin liquid where the electron’s spin-$\frac{1}{2}$ fractionalizes into two degrees of freedom: itinerant Majorana fermions and $\mathbb{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 ($\Gamma$) and Heisenberg ($J$) spin interactionsRau et al. (2014). For example, in $\alpha$-RuCl3 (RuCl3), an actively studied KSL candidate, comprehensive ab initio computations and recent dynamical studiesRan et al. (2017); Wang et al. (2017) suggest that ferromagnetic $K$ and antiferromagnetic $\Gamma$ interactions are dominant and comparable in magnitude, while $J$ is negligibleKim and Kee (2016); Winter et al. (2016); Yadav et al. (2016). The balance of these and additional small further neighbor interactions causes RuCl3 and other KSL candidates to magnetically order at low temperature; however, it is still unclear whether or not the often-large $\Gamma$ interaction prefers magnetic ordering. Meanwhile, the community has attempted to revive the possibility of a KSL in RuCl3 by applying a small magnetic field, with the effect of entering a potential spin liquid phase with no magnetizationYadav et al. (2016).

Since a weak magnetic field takes RuCl3 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}$Kim and Kee (2016); Winter et al. (2016); Yadav et al. (2016) term or terms coming from slight trigonal distortionRau and Kee (2014). This calls into question the role of the $\Gamma$ interaction. Interestingly, a recent analysis of the $\Gamma$ model revealed a macroscopically degenerate classical ground stateRousochatzakis and Perkins (2017).

In this work we will thus investigate if a model with $K$ and $\Gamma$ 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 $\Gamma$, but the resulting phase is not ordereredRau et al. (2014). 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 $\Gamma$ limits, as shown in Fig. 1b.

We consider the following nearest-neighbor (n.n.) model on a two-dimensional honeycomb lattice:

[TABLE]

where

[TABLE]

and $H^{x,y}$ are defined similarly with corresponding $K_{x,y}$ and $\Gamma_{x,y}$. Each $H^{\gamma}$ represents the n.n spin interactions along one of the three bond directions, $\gamma=x,y,z$. The model is parameterized by $K_{z}=-(1+2a_{K})\cos\phi$, $K_{x,y}=-(1-a_{K})\cos\phi$, $\Gamma_{x,y,z}=\sin\phi$, with $a_{K}$ characterizing bond anisotropy. When $\phi=0,\pi$ (i.e. $\Gamma_{\gamma}=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 statesImada and Takahashi (1986); Jaklic and Prelovsek (1994); Hams and De Raedt (2000); Sugiura and Shimizu (2012, 2013), 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 $\Gamma$ limits (see Supplementary Materials (SM)). We present the following results using our numerical techniques:

1. When $a_{K}\geq 0.06$, we find that the $-K_{\gamma}$ ($0\leq\phi\leq\pi/2$) and $\Gamma_{\gamma}$ ($\phi/\pi=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 $\Gamma_{\gamma}$ limit and the antiferromagnetic KSL at $\phi/\pi=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 $\Gamma_{\gamma}$ limit, consistent with our ED results.

4. Zig-zag spin correlations become dominant upon perturbing the quantum spin liquid phase in $0<\phi<\pi/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 $\phi/\pi\in[0,1]$, and for different anisotropy parameters by ED on a 24-site cluster using periodic boundary conditions (see SM). Discontinuities in $\frac{1}{N}\frac{\partial E_{0}}{\partial\phi}$ were used to identify possible phase transitions. Remarkably, when $\phi/\pi\in[0,0.5]$ and $0\leq 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 $\Gamma$-limit ($\phi/\pi=0.5$) is adiabatically connected to the ferromagnetic Kitaev spin liquid ($\phi/\pi=0$) as depicted in Fig. 1b for $a_{K}=0.1$.

In the antiferromagnetic region of phase space, there are two large discontinuities in $\frac{1}{N}\frac{\partial E_{0}}{\partial\phi}$ that encompass a large region of phase space separating the $\Gamma$-limit and the exactly solvable antiferromagnetic Kitaev limit at $\phi/\pi=1$. These peaks coincide with kinks in $E_{0}/N$ (solid yellow) shown in Fig. 2a. Two smaller discontinuities can also be seen near $\phi/\pi=0.75$, however these are not present when $a_{K}=0$, while the larger jumps near $\phi/\pi=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 $\phi/\pi\in[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 $\langle{\bf S}_{i}\cdot{\bf S}_{j}\rangle$, where $i$ and $j$ are site indices on the honeycomb lattice. By Fourier transform, we obtain the static structure factor (SSF) given by $S_{\textbf{q}}=\frac{1}{N}\sum_{i,j}e^{i(\textbf{r}_{i}-\textbf{r}_{j})\cdot\textbf{q}}\langle{\bf S}_{i}\cdot{\bf S}_{j}\rangle$ 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 $\frac{1}{N}\frac{\partial E_{0}}{\partial\phi}$. Visualizations of the SSF over the BZ for representative $\phi$ in the phase diagram are presented in Fig. 2b. The SSF in Fig. 2b is obtained by computing the average of $\langle S_{i}\cdot S_{j}\rangle$ 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 $\phi\in[0,\pi/2]$ and the spin correlations at the $\Gamma$- and M-points are comparable in intensity when $\Gamma\simeq K_{\gamma}$, leading to a “star”-shaped structure in the SSF as seen in Fig. 2b (e.g, $\phi/\pi=0.2$). The extended phase separating $\phi/\pi=0.5$ and $1$ is characterized by dominating spin correlations at the K- and $\Gamma^{\prime}$-points in the reciprocal lattice ($\phi/\pi=0.75$ in Fig. 2b). Contained within this phase is the exactly solvable point with hidden $SU(2)$ symmetry at $\phi/\pi=0.75$ which features K-point correlationsChaloupka and Khaliullin (2015) consistent with the results presented here. Thus the extended spin liquid phase for ferromagnetic K is separated from the antiferromagnetic KSL at $\phi/\pi=1$ by a magnetically ordered phase.

These results can be connected to real materials, particularly RuCl3 in which a zig-zag magnetic ordering has been observedSears et al. (2015); Johnson et al. (2015); Cao et al. (2016). Previous studies have shown that in addition to the n.n. ferromagnetic Kitaev and antiferromagnetic $\Gamma$ interactions, a 3rd n.n. antiferromagnetic Heisenberg interaction $J_{3}\sum_{\langle\langle\langle i,j\rangle\rangle\rangle}{\bf S}_{i}\cdot{\bf S}_{j}$ is non-vanishing and plays a role in determining the magnetic ordering in RuCl3Kim and Kee (2016); Winter et al. (2016). Fig. 3 shows that the effect of perturbing Eq. (1) by $J_{3}$ is to enhance (suppress) the M-point ($\Gamma$-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 $\frac{1}{2}$-plateau in the entropy $S(T)$, which is attributed to the thermal fractionalization of spin degrees of freedomNasu et al. (2015). Here we go beyond the Kitaev limit and investigate the heat capacity and entropy at finite temperature in the presence of $\Gamma$, which is expected to compete with $K_{\gamma}$ in RuCl3, using the method of thermal pure quantum states (see SM).

The dependence of $C(T)$ on $\phi$ 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 $\phi/\pi=0$, and is seen to be maintained continuously as $\phi/\pi$ approaches $0.5$ so that the $\Gamma$-limit shows a qualitatively similar behaviour in $C(T)$ to the Kitaev spin liquid. Evidence for a phase transition can be seen when $\phi\gtrsim 0.7$ on account of the abrupt change in $C(T)$, resembling that of the heat capacity in trivially ordered phases Yamaji et al. (2016). This finding is consistent with our ED results and the 120-order at $\phi/\pi=0.75$ seen in Refs. Rau et al., 2014 and Chaloupka and Khaliullin, 2015. The dependence of $S(T)$ on $\phi$ is plotted in the bottom panel of Fig. 4 with a clear $\frac{1}{2}$-plateau observed when $\phi/\pi=0$, consistent with the expected Kitaev spin liquid behaviour. In addition, a plateau of about 1/5 the total entropy is observed when $\phi/\pi=0.5$. Another plateau is observed in the magnetically ordered phase around $\phi/\pi=0.75$; however, this feature can be attributed to finite-size effects as follows. The $(N+1)$-fold ground state degeneracy at $\phi/\pi=0.75$ due to the hidden $SU(2)$ symmetryChaloupka and Khaliullin (2015) is only slightly broken away from this point, inducing a plateau in $S(T)$ with height given by $\ln(N+1)/N\ln 2\simeq 0.1935\sim 1/5$ when $N=24$. By contrast, the height of the plateau around $\phi/\pi=0.5$ is independent of $N$Not .

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 Nasu et al. (2015). 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 fluxKitaev (2006); Knolle et al. (2014). 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 correlationsNasu et al. (2015). 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 algorithmVidal (2007) (iTEBD; see SM). Classically, the ground state in the $\Gamma$-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 $\Gamma$-limit on the two-dimensional (2D) honeycomb and three-dimensional (3D) hyper-honeycombRousochatzakis and Perkins (2017) lattices also feature similar classical degeneracy. The similarity at the classical level of the $\Gamma$-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 $\chi=10$, and anisotropy $a_{K}=0.1$. In this calculation, we have also introduced an anisotropy to $\Gamma_{\gamma}$ such that $\Gamma_{x}=\Gamma_{y}=(1-a_{\Gamma})\sin\phi$ and $\Gamma_{z}=(1+2a_{\Gamma})\sin\phi$ in order to apply the iTEBD method (see SM). No transition is found when $\phi/\pi\in[0,0.5]$ and the obtained state is a highly entangled paramagnet, with $S_{E}\sim 0.8$ for strong ($z$) bonds, while for weak ($x,y$) bonds, $S_{E}\sim 0.4$. Deep in the gapped phase of the Kitaev model, with large anisotropy $a_{K}$, one finds $S_{E}\sim\log 2\sim 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 fermionsKimchi et al. (2014). An increase is expected in the entanglement entropy as one approaches a phase transition, however no such peaks are seen for $0<\phi<\pi/2$. Similarly, there are no sharp features in the ground state energy $E_{0}$ as a function of $\phi$, which indicates that this phase is adiabatically connected to the Kitaev spin liquid at $\phi=0$.

There is an apparent first order transition around $\phi=0.6\pi$ 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 $\phi/\pi=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 RuCl3 and earlier ab initio computationsKim and Kee (2016); Winter et al. (2016); Yadav et al. (2016). In a recent inelastic neutron scattering experiment on RuCl3, 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 stateBanerjee et al. (2016). 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 $\Gamma_{\gamma}$ limit. When the strength of the ferromagnetic Kitaev interaction and the $\Gamma_{\gamma}$ 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 RuCl3. Given that the ab initio computations suggest comparable magnitudes of the ferromagnetic Kitaev and $\Gamma_{\gamma}$ interactions in RuCl3Kim and Kee (2016), it is conceivable that RuCl3 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 RuCl3.

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 studyGohlke (2017), 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 $\Gamma$, while they may lead to incommensurate ordering for negative $\Gamma$Rousochatzakis and Perkins (2017). Thus it is likely that the positive $\Gamma$ 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 Chaloupka et al. (2010); Rau et al. (2014). 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 statesSugiura and Shimizu (2012, 2013) 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 algorithmVidal (2007). 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$\Phi$Kawamura et al. (2017). 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.

 

Contributions

 

A.C. and Y.Y. performed the exact diagonalization calculations. Y.Y. performed the thermal pure quantum states calculations. G.W. performed the iTEBD calculations. H.-Y.K. and Y.B.K. supervised the study. All authors contributed to the writing of the manuscript.

 

Competing Interests

 

The authors declare no competing financial or non-financial interests.

 

Data Availability

 

All relevant data is available from the corresponding author.

Supplementary Materials

I Exact Diagonalization

Ground state energy and $S_{q}$ in the $a_{K}=0$ limit

When the Kitaev exchange is isotropic $a_{K}=0$, three discontinuities in $\frac{1}{N}\frac{\partial E_{0}}{\partial\phi}$, which are absent with slight anisotropy $a_{K}=0.1$, can be seen when $\phi/\pi\in[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 $\Gamma$- and M-point spin correlations when $\phi/\pi\in[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 $\frac{1}{N}\frac{\partial E_{0}}{\partial\phi}$ and a dramatic change in $S_{q}$ when $K^{\gamma}>0$ are evidence of an extended magnetically ordered phase with dominating K- and $\Gamma^{\prime}$-point correlations. The inconsistent appearance of the other discontinuities in $\frac{1}{N}\frac{\partial E_{0}}{\partial\phi}$ (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$).

Phase diagram for various $a_{K}\neq 0$

We investigated the phase diagram between $\phi/\pi=0$ and $\phi/\pi=0.5$ for different $a_{K}$ as shown in Fig. 8 calculated by ED. The discontinuities seen in $\frac{\partial E_{0}}{\partial\phi}$ when $a_{K}=0$ disappear for $a_{K}\geq 0.06$. The smooth connection between $\phi/\pi=0$ and $\phi/\pi=0.5$ limits is found to be robust for a wide range of $a_{K}$.

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\rangle\equiv|S_{1}^{z}\dots S_{24}^{z}\rangle$, on the 24 site cluster can be translated as $T_{1}^{n}T_{2}^{m}|s\rangle$ where $n=0,...5$ and $m=0,1$. Periodic boundary conditions are imposed such that $T_{1}^{6}=1$ and $T_{2}^{2}=1$. Thus there are 12 unique eigenvalues of $T_{1}^{n}T_{2}^{m}$ (there are two atoms in the unit cell) of the form $e^{2\pi i(k_{x}/N_{1}+k_{y}/N_{2})}$, 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 $\phi$, the ground state is in the $q\equiv(k_{x},k_{y})=0$ momentum sector. The energy spectra corresponding to $\phi/\pi=0.5$ and [math] are shown in Fig. 9a. In the $\Gamma$-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 $\phi$ 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 first-order phase transitions found when $\phi\in[0,\pi/2]$ seen in Fig. 6. The location of these level crossings when $\phi\in[0,\pi/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 $\phi/\pi=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.

$J_{3}$ 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 susceptibilityAlbuquerque et al. (2010) (solid lines) $\chi_{F}$ when $\phi/\pi=0.5,0.1$ and $0.2$. The broad peaks in both $\chi_{F}$ and $-\frac{1}{N}\frac{\partial^{2}E_{0}}{\partial\phi^{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 $\Gamma$, we systematically calculate heat capacity by using typical pure statesReimann (2007). Several pioneering studies Imada and Takahashi (1986); Jaklic and Prelovsek (1994); Hams and De Raedt (2000) 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.Sugiura and Shimizu, 2012, 2013, where such a pure state that replaces a canonical ensemble is called a thermal pure quantum (TPQ) stateSugiura and Shimizu (2012, 2013).

We briefly summarize the construction of TPQ states following Ref. Sugiura and Shimizu, 2012. A TPQ state of $N$ quantum $S=1/2$ spins at infinite temperatures is simply given by a random vector,

[TABLE]

where $|i\rangle$ 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 $\sum_{i=0}^{2^{N}-1}|c_{i}|^{2}=1$. Then, by multiplying a wave vector by the hamiltonian $\hat{H}$, we construct the TPQ states at lower temperatures as follows: Starting with an initial (or the 0-th step) vector $|\Phi_{0}\rangle=|\phi_{+\infty}\rangle$, the $k$-th step vector $|\Phi_{k}\rangle$ $(k\geq 1)$ is constructed as

[TABLE]

where $\Lambda$ is a constant larger than the largest eigenvalue of $\hat{H}$. The above $k$-th step vector is a TPQ state at a finite temperature $T$. The corresponding inverse temperature ${\beta}=(k_{B}T)^{-1}$ is determined through the following formula Sugiura and Shimizu (2012),

[TABLE]

where $k_{B}$ is the Boltzmann constant. In other word, a TPQ state at temperatures $T$ is given as,

[TABLE]

The heat capacity and entropy of $\hat{H}$ are then estimated by using a TPQ state $|\phi_{T}\rangle$. In thermodynamic limit, a single TPQ state $|\phi_{T}\rangle$ 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:

[TABLE]

The entropy $S$ is estimated by integrating $C/T$ from high temperatures as

[TABLE]

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 algorithmVidal (2007) (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 $|\psi\rangle$ may be defined by site-tensors $T_{i_{r}j_{r}k_{r}}^{(r),s_{r}}$, and diagonal bond-matrices $\lambda_{i_{r}i_{r^{\prime}}}^{\langle rr^{\prime}\rangle}$. Here $r,r^{\prime}$ denote sites, $s_{r}=\uparrow,\downarrow$ 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 $\chi$. The TPS is formally given by

[TABLE]

In the TPS canonical form, the diagonal values of $\lambda_{i_{r}i_{r^{\prime}}}^{\langle rr^{\prime}\rangle}$ 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}=-\sum_{i}\lambda_{ii}^{2}\log\lambda_{ii}^{2}$. 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 $\lambda$ 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^{-\Delta\tau H}$, where $H$ is the Hamiltonian, and $\Delta\tau$ is a small time step, is applied using a fourth-order Suzuki-Trotter decompositionSornborger and Stewart (1998), separately evolving each (independent) bond at a time. Each time the evolution operator is applied, the bond dimension (or Schmidt rank) $\chi$ generally increases, making the algorithm unmanageable. To avoid this, a fixed $\chi$ is maintained by applying a truncated singular value decomposition (SVD) to the evolved TPS. Thus, larger values of $\chi$ 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 formOrús and Vidal (2008). We begin by iteratively evolving a random TPS with a time step of $\Delta\tau=0.1$ until the entanglement entropies $S_{E}$ on all bonds have converged. We then reduce $\Delta\tau$ by a half and continue evolving, again until convergence is achieved. This process is repeated until $\Delta\tau<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 anisotropyKitaev (2006). 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.Kimchi et al. (2014) were able to determine that in the $K$ limit there is a gapped phase for $a_{K}\gtrsim 0.05$.