Thermodynamics of 3D Kitaev quantum spin liquids via tensor networks

We study the 3D Kitaev and Kitaev-Heisenberg models respectively on the hyperhoneycomb and hyperoctagon lattices, both at zero and finite-temperature, in the thermodynamic limit. Our analysis relies on advanced tensor network (TN) simulations based on graph Projected Entangled-Pair States (gPEPS). We map out the TN phase diagrams of the models and characterize their underlying gapped and gapless phases both at zero and finite temperature. In particular, we demonstrate how cooling down the hyperhoneycomb system from high-temperature leads to fractionalization of spins to itinerant Majorana fermions and gauge fields that occurs in two separate temperature regimes, leaving their fingerprint on specific heat as a double-peak feature as well as on other quantities such as the thermal entropy, spin-spin correlations and bond entropy. Using the Majorana representation of the Kitaev model, we further show that the low-temperature thermal transition to the Kitaev quantum spin liquid (QSL) phase is associated with the non-trivial Majorana band topology and the presence of Weyl nodes, which manifests itself via non-vanishing Chern number and finite thermal Hall conductivity. Beyond the pure Kitaev limit, we study the 3D Kitaev-Heisenberg (KH) model on the hyperoctagon lattice and extract the full phase diagram for different Heisenberg couplings. We further explore the thermodynamic properties of the magnetically-ordered regions in the KH model and show that, in contrast to the QSL phase, here the thermal phase transition follows the standard Landau symmetry-breaking theory.

The two-dimensional (2D) Kitaev model on the honeycomb lattice [3] is one of the famous examples of a QSL which has played a major role in a deeper understanding of the physics of quantum phases of matter, both theoretically and experimentally. The Kitaev model is a quantum spin system with anisotropic Ising-like bond directional interactions that naturally arises as an interplay of crystal-field and strong spin-orbit coupling in a variety of 4d and 5d materials [20][21][22]. Due to the bonddirectional interactions, the Kitaev model on trivalent lattices is highly frustrated with dominant quantum fluctuations. At zero temperature, the ground state of these Mott insulators forms a highly entangled QSL in which the original spin degrees of freedom are fractionalized into itinerant (non-interacting) Majorana fermions and an emergent static Z 2 gauge field. The resulting QSL is a gapless state which becomes gapped by breaking timereversal symmetry (TRS) [3]. This so-called Kitaev spin liquid is a topologically ordered state that is known to host non-abelian Ising anyons, gapped flux excitations (visions) [3,5], and a chiral gapless Majorana edge mode * saeed.jahromi@dipc.org which gives rise to a quantized thermal quantum Hall effect (QHE) at low-temperature regime [23].
The exact solvability of the Kitaev model remains valid for 3D trivalent lattices [24][25][26]. This has largely motivated the study of the Kitaev model on lattices such as the hyperhoneycomb [27,28] and hyperoctagon [29] (see Fig. 1 and also Ref. [24] for a full list of relevant 3D Kitaev lattices) which are extensions of the honeycomb [3] and square-octagon [30] lattices to 3D. Most importantly, the recent discovery of a 3D material with strong bond directional spin-orbit interactions in β − Li 2 IrO 3 and γ − Li 2 IrO 3 compounds [28,31] has attracted considerable attention to the theoretical and experimental study of 3D Kitaev materials. Similar to the 2D Kitaev honeycomb model, the low-energy physics of 3D Kitaev QSLs is also described as a gapless Majorana metal in the background of Z 2 gauge fields [24]. However, depending on the underlying lattice geometries, the Fermi surfaces of the Majorana metals are described as topological Majorana Fermi surfaces, nodal lines, or Weyl points [24,27,29].
Away from the exactly solvable point or at finite temperature, the analytical tractability of the Kitaev model becomes highly non-trivial or impossible. In such situations, ground-state properties of the model can only be studied by advanced numerical techniques. For example, studying the phase diagram of the Kitaev model in the presence of the Heisenberg interaction, which naturally arises as the next-leading interaction in the spin-orbit Mott insulators [32,33], is highly challenging with analytical techniques. Besides, studying the thermodynamic properties of the Kitaev QSL on different 3D lattices at finite temperatures is numerically challenging. While the 2D Kitaev model on different lattices has been studied largely by state-of-the-art numerical methods such as exact diagonalization (ED) [34][35][36][37], quantum Monte Carlo (QMC) [38], and tensor network (TN) algorithms [39], the study of a generic Kitaev model on different 3D lattices is only limited to mean-field treatment [40], series expansion [41], and more recently a QMC which remains sign-free as long as the gauge fields are static and the Majorana representation remains valid [38,[42][43][44][45].
Although TN methods have been shown to be one of the most promising techniques for accurate simulation of the 2D strongly correlated systems both at zero- [10,[46][47][48][49][50][51][52][53][54] and finite-temperature [55][56][57][58][59][60][61][62][63][64], their application to 3D lattices, and in particular at finite temperature, has largely been left behind mostly due to technical challenges. It is therefore crucial to develop new efficient tools to simulate generic 3D quantum many-body systems that are not directly tractable by, say, QMC methods. In this paper, we use our recent graph-based infinite projected entangled-pair state algorithm (gPEPS) [54] to study the ground-state properties of the 3D Kitaev model in the thermodynamic limit. More specifically, we simulate the spin-1/2 Kitaev model, i.e., Eq.(1) on the hyperhoneycomb lattice, computing their zero temperature phase diagram. Next, we use a variant of our TN algorithm for calculating the thermal density matrix of infinite-size quantum systems, the so-called thermal gPEPS (TgPEPS) [53] and study the thermodynamic properties of the systems at finite temperature. We show that the TgPEPS can faithfully capture the intermediateto-high temperature regimes in the thermodynamic limit, which are the most relevant ones in experimental probes of Kitaev materials [42].
We particularly demonstrate how fractionalization of the original spin degrees of freedom to Majorana fermions and gauge fields leaves its fingerprint on the local observables such as nearest-neighbor correlations and bond entanglement. In order to crosscheck and supplement our TN simulations, we use the Majorana representation of the Kitaev model [3] and extract the thermodynamic properties of the system for the hyperhoneycomb lattice particularly at very low-temperatures, which is the challenging regime for TN algorithms. We further calculate the Chern number and thermal Hall conductivity of the Kitaev model and capture the thermal phase transition beyond which the gauge degrees of freedom are stabilized in the background and the ground state ends up being a highly entangled QSL. Away from the exactly-solvable point, we study the 3D Kitaev-Heisenberg (KH) model, i.e., Hamiltonian (11), in the hyperoctagon lattice, and extract the full phase diagram of the KH model in different regimes of the Heisenberg couplings. We show how different phases and phase boundaries can be identified by measuring different quantities such as magnetization, entanglement entropy, and spin-spin correlations. Our study is further complemented by exploring the thermodynamic properties of different magnetically ordered regions in the phase diagram of the KH model. The paper is organized as follows: In Sec. II we introduce the Kitaev Hamiltonian and review the ground state properties of the model on different 3D lattice structures. In Sec. III we discuss the details of our TN algorithm for simulating both ground state and the thermal density matrix of local Hamiltonians on any arbitrary lattice structure. Next in Sec. IV A we discuss the TN phase diagram of the pure 3D Kitaev hyperhoneycomb model at zero-temperature. The thermodynamic properties of the Kitaev QSL at finite-temperature are discussed in Sec. IV B, and the 3D Kitaev-Heisenberg hyperoctagon model is studied in Sec. V. Finally, Sec. VI is devoted to the discussion and conclusions.

II. MODEL
The Kitaev model was first introduced by Alexei Kitaev on the 2D honeycomb lattice in the context of topological quantum computation [3]. The intriguing properties of the model soon attracted huge attention from the condensed matter perspective since it was suggested that the ground state of the system is a QSL [1]. The Hamiltonian of the Kitaev model is given by where the sum runs over three subclasses of bonds labeled by γ = x, y, z, denoting the three Ising-like interactions on the corresponding link. While the original Kitaev model was introduced using spin-1/2 Pauli operators, our definition here is based on generic spin operators.
The Kitaev model is exactly solvable on trivalent lattices, i.e., lattices with three links connected to each site. Using the local transformation S γ i = i 2 b γ i c i , the original spin operators can be represented by four Majorana fermions. The original interacting spin model takes then the bilinear form, whereû γ ij = −û γ ji = ib γ i b γ j are bond operators obtained from grouping Majoranas b γ along nearest-neighbor links i, j . The bond operators commute with each other and with the Hamiltonian (2). They are, therefore, conserved integrals of motion which square to identity and their eigenvalues are given by u ij = ±1. Within this parton construction, the new degrees of freedom are now the non-interacting Majorana fermions c i coupled to a static Z 2 gauge field u ij . Given a fixed gauge configuration {u ij }, the spectrum of the system can be readily extracted by diagonalizing the skew-symmetry matrix iA. Its eigenvalues come in pairs ± µ as a result of the inherent particle-hole symmetry of the (single-particle) Majorana Hamiltonian (2). The ground state of the system is thus given by a configuration of gauge fields which minimizes the energy. According to Lieb's theorem, depending on the lattice structure, its symmetry, and dimensionality, the ground state is generically given by 0-or π-flux configuration of gauge fields.
While Lieb's theorem is valid for 2D trivalent lattices, its extension to 3D trivalent structures is hampered since the symmetry requirement for the theorem is not fulfilled in general, except for hyperhoneycomb lattice. One may therefore use numerically exact QMC simulations to unambiguously determine the ground-state gauge pattern. A careful analysis in this regard shows that all bipartite trivalent 3D lattices still follow the general conclusions of Lieb's theorem (in spite of not being applicable) in the absence of additional geometrical constraints leading to "gauge frustration" (see Refs. [42] for details).
The zero-temperature phase diagram of the Kitaev model on the 2D honeycomb lattice is known to host two QSL phases. The first one is gapless and emerges near the isotropic point of Hamiltonian (1), i.e., K ≡ K x = K y = K z . The second one is actually three equivalent gapped phases that arise whenever one of the couplings dominates the other two. The low-energy effective theory of the three gapped phases is given by an abelian topological phase known as the toric code [5]. The gapless phase is also known to become gapped by applying time-reversal symmetry (TRS) breaking perturbations to the Hamiltonian (1) such as a uniform magnetic field, − i B. S i along the [111] direction. When the field strength B is small compared to the vison gap ∆, one can derive a low-energy effective model using third-order perturbation theory in the ground-state flux sector, yielding where the three-spin coupling constant κ ∼ B x B y B z /∆ 2 encodes the strength of the effective magnetic field, and triples i, j, k indicate three neighbouring sites with strictly different bond-type α, β and γ.
In the effective Majorana representation (3), the threespin interaction turns into a next-nearest-neighbour hopping term between sites i and k connected by site j. In the 2D Kitaev models, such a hopping term opens a gap at the Dirac points of the bulk spectrum and yields a chiral edge mode with a non-vanishing Chern number ν (other third-order terms are irrelevant to the Diracgap opening in the renormalization group sense). The resulting gapped state is Z 2 topologically ordered with non-abelian anyon excitations.
In the low-temperature regime, the thermal Hall conductance of the 2D Kitaev model shows a quantized value with respect to the field strength κ, i.e., I xy = π 12 νT , which is exactly half of the two-dimensional thermal Hall conductance in the integer QHE. Such a half-integer quantization is a signature of a topologically-protected chiral edge mode of charge-neutral Majorana fermions, which have half the degrees of freedom of conventional fermions.
The Kitaev model on 3D lattices still shares some of the properties of the 2D honeycomb version, including a rich variety of gapless Z 2 QSL phases. The main difference with the 2D version though is that, according to the symmetry classification of free-fermion systems [65], breaking the TRS will not give rise to topologically protected gapped phases for 3D Kitaev models. Despite not being topological in a strict sense, these models exhibit however a finite but non-quantized thermal quantum Hall effect in the presence of a magnetic field (see Sec. IV B). Moreover, the nodal structure of the Majorana fermions on different 3D Kitaev lattices have distinct properties such as nodal lines with a flat surface band, topologically protected Weyl points with the so-called surface Fermi arcs, and a Majorana semi-metal with a Majorana Fermi surface characterized by a finite Majorana density of states at lowest-lying energy levels [24].
In this paper we focus on two specific examples of 3D Kitaev structures, namely the hyperhoneycomb and hyperoctagon lattices, shown in Fig. 1. Both lattices are among the most studied examples of 3D Kitaev materials. While the hyperoctagon lattice is the 3D extension of the square-octagon lattice, the hyperhoneycomb is the natural extension of the honeycomb lattice to 3D. The latter has a experimental realization in β − Li 2 IrO 3 compounds [28,31]. Besides, both hyperhoneycomb and hyperoctagon lattices contain two independent types of loop operators of length 10 around their plaquettes which guarantee a ground state with vanishing total flux [24,42].
While the phase diagram of both lattices is composed of a gapless region around the isotropic point surrounded by three gapped phases at the corners (see Fig. 4), they have different Majorana Fermi surfaces. The gapless region of the hyperoctagon phase diagram harbours two distinct Majorana Fermi surface regions: a trivial one and another region with topological protected Fermi surface enclosing a Weyl node at finite energy. In contrast, the gapless modes in the bulk spectrum of the hyperhoneycomb lattice form a closed line of Dirac nodes which is protected by TRS. Breaking the TRS will gap out the Majorana Fermi line, leaving pairs of Weyl points in the bulk, exactly at zero energy, and gapless Fermi arcs on the surface [27].

III. METHODS
Let us now briefly sketch how TNs can be used to obtain the phase diagram and thermodynamic properties of the Kitaev model on the aforementioned 3D lattices at zero and finite-temperature.
A. Tensor network algorithms for 3D lattices TN methods [66][67][68][69][70][71] have played a major role in the discovery of many exotic phases of matter in recent years. The basic idea at work in these methods is to write the wave function of local quantum many-body (QMB) Hamiltonians in a very efficient way based on their entanglement structure [66]. Matrix product states (MPS), which provide a natural TN representation of 1D QMB systems, are probably the most famous example of TN states, since they are at the core of well-known algorithms such as density matrix renormalization group (DMRG) [72,73]. Another well-known example is that of projected entangled-pair state (PEPS) [66,67,71,74] and its variant for infinite lattices, i.e, infinite-PEPS (iPEPS).
The TN representation for the ground state of a local lattice Hamiltonian can generally be obtained by variational methods [75] or approaches based on evolution in imaginary-time [76][77][78]. In spite of their great promise and success, these techniques are restricted due to the implementation challenges. Examples of this are the lattice geometry as well as the efficient approximation of environments [76,79]. During past years tremendous progress has been put forward to simulate different 2D lattice structures with PEPS. However, there have been not so many examples of 3D simulations with TNs. Thanks to our recent graph-based iPEPS algorithm (which we call gPEPS) [54], we have managed to simulate local quantum lattice models in arbitrary dimensions and lattice geometries.
The gPEPS algorithm relies on the so-called structurematrix (SM) construction, which resolves the geometrical implementation challenges [54,64]. The SM provides a compact way for storing the connectivity information of the underlying TN, i.e., each column of the SM corresponds to one of the links of the lattice and contains all the details about the neighbouring tensors, their interconnecting indices, and their bond dimensions. One can then fully automatize the TN update by looping over the columns of the SM in a very systematic way, without the burden of complications due to geometry (see Refs. [54] for detailed discussions and examples of SM for different lattices). The ground state of local nearestneighbour Hamiltonians is further approximated by using the simple-update (SU) technique based on imaginarytime evolution (ITE). Within the SU scheme, the environments of local tensors are given by bond matrices obtained from local singular-value decompositions (SVD), which provide a mean-field-like approximation to the environment and correlations around local tensors. While more accurate techniques such as full-update (FU) [47][48][49] have also been designed to capture the full correlation in the system, their implementation and truncation is largely limited by lattice geometry particularly for 3D lattices. Nevertheless, the mean-field approximation of the environment in the SU has been shown to be very good for higher-dimensional systems as well as thermal states, in turn making the SU a quite accurate option in these situations.
In order to study the thermodynamic properties of QMB systems at finite-temperature we target the thermal density-matrix (TDM) of the corresponding Hamiltonian H, i.e., ρ β = e −βH , β = 1/T being the inverse temperature. The TDM of the system (see Fig. 2 for the TDM of the hyperhoneycomb lattice) can then be approximated by evolving in imaginary-time for a time β/2 both the bra and ket degrees of freedom starting from the maximally-entangled infinite temperature state, i.e., ρ β = e −βH/2 · I · e −βH/2 . Similar to the ground state simulation, our thermal gPEPS algorithm (TgPEPS) [64] uses both SU and mean-field environment for local optimization of the TDM in arbitrary dimensions. The expectation values of local operators and correlators can then be calculated as Ô β = Tr(ρ βÔ )/Tr(ρ β ) where, Tr denotes the trace operation (see Ref. [64] for more details on the method).
In order to simulate the 3D Kitaev model on the hyperhoneycomb and hyperoctagon lattices, we used translationally invariant unit-cells of 32-sites. Our SU optimization technique for both gPEPS and TgPEPS was based on imaginary-time evolution accompanied by proper truncations at the level of local SVDs. Our algorithms have further been stabilized by proper choice of gauge-fixing and super-orthogonalization of local tensors [54].

IV. HYPERHONEYCOMB KITAEV SPIN-LIQUIDS
Let us start by presenting our results for the 3D Kitaev model in the hyperhoneycomb lattice at zero temperature. Results for the hyperoctagon lattice are very similar nd can be obtained analogously.
The full phase diagram of the 3D Kitaev model on both Hyperhoneycomb and Hyperoctagon lattices pretty much looks the same as Fig. 4 which is composed of a gapless phase surrounded by three gapped regions. The phase boundaries can be captured analytically by locating the location of the closure of the gap of the spectrum of the quadratic Majorana Hamiltonian (2) in the momentum space [3]. In the TN representation, we only have access to the ground state wave function of the 3D Kitaev model in the real space. The phase boundaries, however, can still be located by different physical observables and their derivatives such as the 2nd-derivative of the ground state energy per-site ε 0 , as shown in Fig. 3 for the hyperhoneycomb lattice. The two curves show the behaviour of ∂ 2 J ε 0 along scan lines J with fixed K z = 0, 0.6 in the K x + K y + K z = 2 plane (dashed lines in Fig. 4). The 2nd-derivatives of the energy show sharp discontinuities when crossing a phase boundary, being this a clear signal of second-order quantum phase transitions.

B. Kitaev QSL at finite-T
It has already been known that the specific heat C v of the Kitaev model on generic lattices exhibits a doublepeak behaviour [38,44,45] at two different temperature regimes. One peak corresponds to the order of the Majorana bandwidth at T ∼ K, and the other is correlated with the size of the vison gap, typically of the order T c ∼ K/100 (see the schematic diagram of Fig. 5). This double-peak feature is indeed a beautiful signature of spin fractionalization to Majorana fermions and a gauge field that occurs in two separate temperature regimes, leaving their fingerprint on the C v curve. In the following, we discuss the thermodynamics of the 3D Kitaev model on the hyperhoneycomb lattice and present our tensor network results. In particular, we show how the signatures of spin fractionalization can be observed in different physical quantities such as spin correlations, thermal entropy, and bond entanglement entropy. In the very low-temperature regime T ∆, where the TN simulations are challenging, we provide the results for the Chern number and thermal Hall conductivity directly from the single-particle spectrum, in order to offer a complete picture for the thermodynamics of the 3D Kitaev model.

C. High-temperature spin-ordering crossover
Let us start our discussion by focussing on the isotropic point in the gapless region of the phase diagram, i.e., Approaching from the hightemperature regime where the system is in a spin paramagnet (spin gas) phase, we first hit the second peak in the specific heat, C v = ∂ε 0 (T )/∂T . For the case of the hyperhoneycomb lattice this is located at the crossover temperature T ≈ 0.256 (see Fig. 6-(b)).
In order to understand the nature of this peak, let us remind that in the parton configuration, the Kitaev model is characterized by Majorana fermions, which are local objects, coupled to a Z 2 gauge field (see also Sec. II). In this language, the kinetic energy of the Majorana fermions, i c i c j , is precisely equivalent to the nearestneighbor (NN) spin correlation, S γγ (T ) ≡ S γ i S γ j T . We, therefore, calculated the NN spin correlation of the Kitaev model on the hyperhoneycomb lattice as shown in Fig. 6-(c). Interestingly, the high-temperature crossover coincides with the onset of NN spin correlations at T ≈ 0.256, in the same location as that of the second peak of the specific heat, which is the temperature at which the spin degrees of freedom begin to be locally fractionalized (see also Fig. 5). Our TN simulation further shows that the system is a conventional paramagnet in the high-temperature region T > T and the NN spin correlation obeys the Curie-Weiss behavior, S γγ ∝ K γ /T . In the opposite extreme limit, however, the NN spin correlation reaches its T = 0 saturation value, i.e., S zz = 0.5248/4, characteristic of the gapless Kitaev QSL [80].
This NN spin ordering is a purely local phenomenon, and its features and its location do not strictly depend on things like the lattice size or geometry [38,42,45]. It is therefore a thermal crossover which can be revealed by the second peak of the specific heat. In our TN simulations, this thermal crossover at the second peak can be better identified using entanglement scaling. Fig. 8-(a) demonstrates the scaling of the location of the hightemperature peak T versus inverse bond dimension 1/D. One can clearly see that the location of T is invariant with respect to the change of bond dimension, implying the crossover nature of T , instead of being a true thermal transition.
The thermal entropy of the system is also given by where S ∞ = ln 2 is the maximum entropy density corresponding to an equally weighted, infinite-temperature Gibbs state. As expected, the thermal entropy saturates at S ∞ in the deep paramagnetic phase. However, the system releases exactly half of its entropy at the crossover temperature T and reaches the plateau ln 2/2. This entropy is related to the contribution of the Majorana fermions to the specific heat of the Kitaev model which is released by moving from hightemperature to the low-temperature regime across the crossover point. The plateau of thermal entropy for the Kitaev model on the hyperhoneycomb lattice is shown in Fig. 6-(d).
In the corners of the 3D Kitaev phase diagram, in the toric code limit, the system is fully gapped. It is, therefore, reasonable to expect that the crossover to high-temperature paramagnetic spin gas occurs at a relatively larger temperature. The lower panels (e)-(h) of Fig. 6 illustrate the energy, specific heat, NN spin correlation, and thermal entropy of the Kitaev model for K x = 0.8, K z = K y = 0.1 on the hyperhoneycomb lattice. The crossover temperature for this gapped region is located at T = 0.338 as shown in the plot of specific heat ( Fig. 6-(f)). Similar to the isotropic point, here the NN spin correlation saturates as well to S zz ≈ 0.248 for T < T down to zero temperature which is accompanied by release of half of the entropy, implying again the fractionalization of original spins.
Last but not least, we have benchmarked our TN simulations against the corresponding quantities extracted from the effective tight-binding Hamiltonian (2). By assuming static fixed-flux configurations, thermal fluctuations of the gauge fields can be ignored and only Majorana fermions are retained. It is then straightforward to see that in the Majorana basis, the thermal average over   (2) ln (2) spin correlation and energy reads where ψ µ (i) is the i th component of the normalized complex eigenvector of the Kitaev Hamiltonian, and the sum on µ runs over half of the single-particle spectrum with non-negative µ > 0. The solid lines in all panels of Fig. 6 correspond to the results for the ground-state gauge ansatz, i.e., the free-flux sector, as well as those averaged over random flux configurations. We find remarkable agreements between the disorder-averaged data and the TN results, in particular above the gauge-fielddisordering transition. Besides, the results of free-flux and random gauge configurations become indistinguishable above the spin-disordering crossover in the thermal paramagnet phase.

D. Low-temperature gauge-ordering transition
Next, we discuss the thermodynamics of the Kitaev model on the hyperhoneycomb lattice below the crossover temperature, i.e., T < T . As already pointed out previously, the spin degrees of freedom are fractionalized to the itinerant Majorana fermions and a Z 2 gauge field at the crossover temperature T . While the formation of Majorana fermions can be captured perfectly by NN spin correlation (since the fermions are local), there is however no direct local measure to reveal the thermodynamics of gauge fields. The gauge fields in the 3D Kitaev model are defined by non local strings which go around the plaquettes of the hyperhoneycomb lattice. Associating a nonlocal flux operator to each plaquette p of the lattice as the gauge structure can be captured by the ±1 eigenvalues of the W p operators, which are integrals of motion of the Kitaev Hamiltonian (1). The QSL ground state of the Kitaev model at T = 0 is therefore distinguished by eigenvalues W p = 1 = e i0 (−1 = e iπ ), indicating the presence of a Z 2 gauge field on the plaquettes with uniform zero (π) flux. Below the low-temperature phase transition, i.e., T < T c , the system is therefore expected to be in a state with an ordered gauge structure. By increasing the temperature above the T > T c , thermal fluctuations break the patterns of gauge loops and allow the formation of loops with different shapes and sizes, hence a disordered gauge background with W p = 0 is stabilized in the system. The intermediate disordered gauge region continues to persist for T c ≤ T ≤ T until it is totally thermalized above the crossover temperature to the trivial paramagnet phase.
The low-temperature gauge ordering transition is revealed by the first peak in the specific heat of the Kitaev model as sketched in Fig.5. Previous Monte Carlo studies [38,[42][43][44][45] have revealed that the thermal gauge ordering transition is located at very low temperature of the order 10 −2 K − 10 −3 K depending on the lattice size, geometry and spatial dimension. In contrast to the crossover temperature, the low-temperature gauge ordering is an actual thermal phase transition.
The local simple-update and the mean-field environ- ment that we used in the TgPEPS algorithm does not allow us to accurately simulate very low temperatures of the order T T c ∼ K/100 due to poor convergence of the SU in this temperature regime. We where therefore unable to calculate the expectation value of the W p operators as a direct probe for capturing the gauge ordering. However, as an indirect signature of the entangled gauge loops around the plaquettes, we calculated the bond entropy, The bond entropy is zero in the high-temperature regime in the spin paramagnet phase, as expected from the infinite-temperature limit. As the system is cooled down and the spin fractionalization occurs at the crossover temperature, the bond entropy starts to grow until it reaches the thermal transition point at T c below which the system is in the Kitaev QSL phase. While our TN specific heat ( Fig.6-(b)) is insensitive to the first peak at T c (due to local SU update and finite bond dimension), it is remarkable that the location of the gauge ordering transition can be captured from the derivative of S b (7-(b)). In contrast to the crossover temperature, physical observables, and their derivatives are sensitive to the bond dimension, and are suitable for entanglement scaling (see (a)), T c depends clearly on the bond dimension, which is a typical signature of a diverging correlation length and hence a true quantum phase transition. Our results for the largest considered bond dimension, D = 7, yield T c = 0.006, which is slightly away from the previous QMC estimate T QM C c = 0.0024 [42]. However, our simulation also shows a tendency to reach the QMC value by increasing the bond dimension D.
Let us further note that the remaining thermal entropy of the system will be further released at the thermal transition point due to the full ordering of the gauge structure. The entropy plot will then show another drop from the ln 2/2 plateau to S T = 0 [42]. Unfortunately, due to the limited convergence of our TN results for T T c ∼ K/100, we are not able to capture the first entropy release, caused by the thermal fluctuation of vi-sons.
Between the two thermal transitions, there is an intermediate temperature regime that spans over two orders of magnitude, T c < T < T ∼ K and might, in fact, be the most relevant temperature regime in experimental probes of Kitaev materials. In this regime, one expects to observe the first signatures of fractionalization with the original spins already broken apart into Majorana fermions and a Z 2 gauge field. The latter, however, is still highly disordered in the intermediate regime which prevents the formation of a clean Majorana band structure (as it is the case at strictly zero temperature or, more precisely, the low-temperature transition). Instead one expects to see a disordered Majorana metal (thermal metal), which has already been observed numerically in certain 2D settings [38,[43][44][45].
While the crossover mechanism at the hightemperature regime is linked to the spin ordering, the low-temperature thermal transition is associated with the fluctuations of gauge fields. As previously pointed out, the ground state of the Kitaev model (at the isotropic point) on different 3D lattice geometries is given in the zero-or π-flux sector, which are sectors with loop-like objects constructing the boundaries of closed volumes. At the critical point T c , these loops break apart and span the whole lattice. The transition between these two regimes has been argued to be of the second-order type, and in the 3D Ising universality class. Moreover, recent QMC simulations have shown that the thermal phase transition in the gapped regions of the phase diagram of the Kitaev model, e.g., K ≡ K x = K y 1, K z → 1, occurs at even lower temperature, of the order T c = 1.925K eff K/100 with K eff ∝ 7K 6 256K 5 z [38], which is beyond the reach of both our TN and exact calculations.
Here we have calculated the bond entanglement entropy of the gapped phase at [0.8, 0.1, 0.1] (see Fig. 7-(c)). As seen in the figure, this shows a similar behaviour to the isotropic case, indicating an increase of correlations from the paramagnet regime to the gauge disordered region. We expect that at the thermal gauge ordering transition, S b will show another increase due to the QSL ground state at T = 0.

E. Thermal Hall effect
Let us now show how the gauge-ordering thermal transition at the isotropic point can be captured precisely by thermal Hall conductivity and Chern number. As we pointed out in Sec. II, Breaking TRS by a field term opens a gap at the Dirac points of the bulk spectrum of the model and and yields a chiral edge mode with a non-vanishing Chern number, giving rise to a finite thermal Hall conductance when subjected to a thermal gradient. In order to shed light on the relationship between the Majorana band topology and thermal Hall response, we evaluate the Chern number of the effective Majorana Hamiltonian (3) by slicing the 3D Brillouin zone (BZ) into 2D planes passing through three points k = (0, 0, k z ), k + q 2 /2 and k + q 3 /2, where q 2 and q 3 are reciprocal lattice vectors of hyperhoneycomb lattice. This allows to define the total Chern number for any fixed value of k z , where Ω nk = i ∇ k u nk | × |∇ k u nk is the non-Abelian Berry curvature for the n th band with energy dispersion nk in the free-flux sector, and |u nk denotes the corresponding eigenvector. Here the integral is taken over a single slice of the BZ parameterized by the momentum k z , and the sum runs over all occupied bands with negative energy. The definition of the Chern number can further be extended to finite temperature as ν(k z , T ) by replacing Ω z nk → f ( nk , T )Ω z nk in Eq. (7), i.e., Ω z nk is weighted by the Fermi-Dirac distribution, f ( , T ) = 1/(1 + exp( /T )). Fig. 9(a) shows the k z -dependence of the Chern number as a function of temperature at a fixed magnetic field strength κ/K = 0.001. Such a small non-zero value of κ allows to capture the low-temperature characteristic features of the system in the absence of magnetic field, for which QMC results are already at hand. In the high temperature region, the system is topologically trivial with ν(k z ) = 0 due to Berry curvature cancellation of almost equally populated bands. As the system is cooled down to the gauge-ordering temperature T c , ν(k z , T ) begins to show characteristic features that are expected to appear at exactly zero temperature: once a plane passes through a gapless Weyl node, located at (0, 0, ±k 0 ), the Chern number jumps discontinuously by an amount given by the charge of the Weyl point (see the inset of Fig. 9(a)).
Next, we discuss the thermal Hall effect, in order to investigate how the changes in topological nature of the system and the presence of Weyl nodes can manifest themselves via nontrivial transport features. The expression for the thermal Hall conductivity of a general noninteracting Majorana Hamiltonian is given by [81][82][83][84] (we set K B = = 1 for simplicity), where V is the volume of the system, f ( , T ) is the Fermi-Dirac distribution, and, is the zero-temperature anomalous Hall coefficient for a system with the chemical potential . According to Eqs. (8) and (9), a finite Berry curvature is the essential ingredient for generating the Hall conductance. It is also instructive to look at the momentum-resolved representation of the thermal Hall conductance (MRHC) on each 2D plane of the sliced BZ, to identify the contribution of each plane separately. The MRHC can be defined as κ xy (T ) = dk z F (k z , T ). Figs. 9(b,c) show the k zvariation of the integrand F (k z , T ) and the corresponding κ xy (T ), respectively. One can clearly see that the detailed feature of F (k z , T ) is strongly tied to the behavior of the Chern number ν(k z , T ). Remarkably, the F (k z , T ) exhibits an abrupt change in the vicinity of the transition temperature T c , and reaches a finite value in the momentum interval with nontrivial topological features, i.e., |k z | ≤ k 0 . Lastly, in Fig. 9(c), we investigate the behaviour of κ xy in the low-temperature limit, indicating a linear T dependence with a finite but non-quantized coefficient. To leading order in T → 0, the expression of thermal Hall conductivity in Eq. (8) is reduced to lim T →0 κ xy = πT 12 σ xy (0).
The thermal Hall coefficient of such Weyl spin liquids is proportional to the anomalous Hall coefficient at the Fermi energy, σ xy (0) = dkz 2π ν(k z ), which is proportional to the distance between the Weyl nodes (or equivalently the length of surface Fermi arcs), and hence is not generally an integer topological invariant. The appearance of such non-quantized plateau is reminiscent of Weyl superconductors [85], e.g., in engineered heterostructures with alternating conventional (s-wave) superconductor and topological insulator layers. This behavior can also occur spontaneously (with no applied magnetic field) in some variants of 3D chiral superconductors [86], and 2D field-driven U(1) spin liquids with Dzyaloshinskii-Moriya interactions [87].
It is worth noting that both the precise location and the total number of Weyl point depend on the strength of the magnetic field, which in turn can leave its finger-print in the magnitude of low-temperature thermal Hall plateau. This can be best perceived from the inset of Fig. 9(c), where the value of thermal Hall coefficient (in units of π/12) is shown as a function of magnetic field κ. Upon increasing κ < 1, the value of thermal Hall plateau remains fairly constant up to κ c = 1 2 3 5 [27], at which it drastically drops due to the change in the total number of Weyl points from two to six (each pair of Weyl points splits into three, conserving zero net chirality).

V. HYPEROCTAGON KITAEV-HEISENBERG MODEL
It has already been argued that the interplay between the crystal-field and spin-orbit coupling in iridates not only results in the emergence of bond-anisotropic Kitaev interactions but also induces isotropic exchange interactions of Heisenberg type [21,32,33]. It is therefore of particular interest to study the interplay between the Kitaev interaction and Heisenberg exchange coupling. Let us further note that the applicability of QMC simulations is restricted to those 3D Kitaev models whose parton description does not exhibit a sign problem, being this a feature that generally gets lost due to the Heisenberg interaction. This motivates us to analyze the ground state and finite-temprature properties of the Kitaev-Heisenberg (KH) model on 3D lattices in the thermodynamic limit. The Hamiltonian of the KH model is given by   where K = sin θ and J = cos θ are the Kitaev and Heisenberg exchange couplings, respectively.
In the following we consider the hyperoctagon lattice. We first elaborate on the zero-temperature phase diagram of the model, and discuss different magnetically ordered phases hosted in the vicinity of the QSL region in the full parameter space of Hamiltonian (11), i.e, θ ∈ [0, 2π]. Additionally, we provide further insight into the finite temperature phase transition in the magnetically ordered phase.
A. T = 0 Phase Diagram Lets us start by investigating the limiting cases of the KH Hamiltonian (11). In the extreme regime where the Heisenberg interaction is switched off, i.e., θ = π 2 (θ = 3π 2 ), the Hamiltonian (11) reduces to the antiferromagnetic (AFM) (ferromagnetic (FM)) pure Kitaev model and the ground state of the system in these regimes is given by a Z 2 QSL phase. In contrast, the opposite limit where θ = 0 (θ = π) the system host a trivial antiferromagnetic (ferromagnetic) phase. Additionally, it has been shown that there is a duality (so-called Klein duality [32,88]) between specific points (angles) of the phase diagram of the Kitaev-Heisenberg models: for specific angles θ there is a transformation to another anglẽ θ such that tanθ = − tan θ − 1. The mapping immediately reveals that by transforming the exact FM angle at θ = π, another hidden SU(2)-symmetric point is revealed at θ = − π 4 for which the Kitaev term vanishes. Solving the resulting Heisenberg Hamiltonian at θ = − π 4 yields a FM ground state which, after transforming back to the original spins, maps to a stripy state. Thanks to this duality, we find that the AFM point (θ = 0) is also isomorphic to θ = 3π 4 , which after the same transformation and rotations leads to the zigzag state.
A sketch of the magnetically ordered phases on the hyperoctagon lattice is shown in Fig. 10. One can empirically check the Klein duality by considering a 4-sublattice transformation between the FM and stripy phases as well as the AFM and the zigzag states. Let us further note that this transformation holds true for both 2D and 3D bipartite lattices, with the exception of those systems in which the AFM phase is frustrated.
Away from these points, the KH model is no longer exactly solvable. We therefore resort to TN simulations in such regimes. Let us further stress that the KH model is not tractable by QMC, except at the pure Kitaev limit θ = π 2 , 3π 2 . Even in this limit, the exact analytical solution is not available due to the inapplicability of Lieb's theorem on the hyperoctagon lattice [24,29,43]. However, the 32-site unit cell used in the thermodynamiclimit gPEPS/TgPEPS methods allows for the Klein duality and is thus expected to capture all the symmetrybreaking patterns in both zero and finite-T phase diagrams. on a number of different signatures such as the ground state energy and its derivatives, magnetization and bond entropy, as shown in Fig. 11-(a-c) ] regions, respectively in the phase diagram of the KH model. The phase boundaries of the transition between the QSL phases and the ordered phase as well as the transition between the magnetically ordered phases seems to be of first-order, since it shows a discontinuity in the first-order derivative of the ground state energy (also visible with sharp discontinuous peaks in the second derivative of energy as depicted in the inset of Fig. 11-(a)). These findings are in agreement with previous mean-field studies on other 3D Kitaev structures such as the hyperhoneycomb lattice [40]. Besides, the sharp discontinuity at other quantities such as total magnetization ( Fig. 11-(b)) and S b (Fig. 11-(c)) is typical of first-order transitions.
The two spin liquid phases in the phase diagram are best identified by a vanishing total magnetization as shown in Fig. 11-( Fig. 11-(b), we examined the QSL region with our TN simulations with more resolution (tiny steps between the couplings) and confirmed the location of the QSL phase boundaries with the gPEPS technique (not shown here). See a similar simulation for the hyperhoneycomb lattice in Ref. [54].

B. Magnetically ordered phases at finite-T
In Sec. IV A and IV B, we investigated the QSL phase of the pure Kitaev model for the AFM exchange couplings. In the absence of external field, the FM QSL ground state has also been shown to have the same ther-modynamic properties as that of the AFM case [35]. Here we further elaborate on the thermodynamics of the ordered phases in the phase diagram of the KH model on the hyperoctagon lattice. More specifically, we will focus on the thermal phase transition in the FM and AFM ordered regions of the phase diagram in the thermodynamic limit, keeping in mind that the stripy and zigzag phases respectively are related to the FM and AFM states according to the Klein duality.
In contrast to the QSL phases, in which the doublepeak nature of the specific heat is associated with the two species of elementary excitations (namely, itinerant Majoranas and visons) with different energy scales, the magnetically ordered phases of the KH model are trivial, showing no signature of fractionalization. We, therefore, expect to observe only one peak in the specific heat of the KH model in the finite-temperature regime of the ordered regions.
In what follows, we fix θ = 34 • and θ = 217 • deep in the AFM and FM phases, respectively, and use the temperature T as tuning parameter. Figure. 12-(a-c) shows the energy, specific heat, and the bond entropy in both the FM and AFM regions of the KH model on the hyperoctagon lattice obtained with the TgPEPS method. As expected, we observe only a single transition between the zero-temperature spin-ordered and high-temperature spin-disordered phases. This transition happens at T AF c = 0.51 and T FM c = 0.53 for the AFM and FM cases, respectively. Cooling down the system from the high-temperature regime spontaneously breaks the spin inversion symmetry down at T c for both FM and AFM phases, indicating a thermal phase transition between the high-and low-temperature regimes. Besides, the T c in both FM and AFM phases is of the same order of magnitude, indicating a similar robustness in the presence of thermal fluctuations.

VI. CONCLUSIONS AND DISCUSSION
The Kitaev model is one of the first examples of quantum spin liquids with fascinating properties such as topo-logical order and fractional excitations. While the 2D version on the honeycomb lattice has been shown to host a non-abelian topologically ground state at zero temperature [3], the 3D version is known to have distinct properties in the nodal manifold such as topologically protected Weyl nodes [24]. Most importantly, it has been suggested that iridate compounds such as β − Li 2 IrO 3 and γ − Li 2 IrO 3 are relevant platforms for experimental realizations of 3D Kitaev materials [28,31]. On top of that, the Kitaev model is a fascinating playground for a deeper understanding of the fractionalization mechanism which occurs in QSLs when the system is cooled down to zero temperature. It is the perfect theory lab to study quantum phases of matter.
While the 2D Kitaev model has been widely studied both analytically and numerically, its 3D version is less explored due to the lack of efficient numerical techniques capable of simulating 3D structures. Current state of the art QMC techniques can only simulate the pure Kitaev model in specific gauge sectors, where the system remains sign-free, and are inapplicable in the presence of relevant perturbations such as Heisenberg interactions and magnetic fields. Other techniques such as exact diagonalization or mean-field approximation are also not free from limitations, such as large finite-size effects and/or deficiency in capturing thermal properties of experimental relevance.
In this paper we proposed an advanced and efficient tensor network algorithm based on graph-projected entangled-pair state for simulating both the ground state and thermodynamic properties of 3D Kiatev qantum spin liquids. In order to demonstrate the accuracy and efficiency of our algorithms, we simulated the 3D Kitaev model on the hyperhoneycomb lattice and mapped out the full phase diagram of the system at zero temperature. We further simulated its thermal density matrix and investigated the thermodynamic properties of both gapless and gapped regions in the phase diagram. We showed how the original spin degrees of freedom are fractionalized to itinerant-Majorana fermions and a static gauge field, leaving their fingerprints on quantities such as specific heat, spin correlations, thermal entropy, and bond entropy. In particular, we found that by cooling down the system from the high-temperature regime, where the system is in a spin paramagnet phase, we first hit a crossover temperature at T ∼ K (K being the Kitaev exchange coupling), below which the spins are fractionalized to a spin-ordered phase equivalent to itinerant-Majorana fermions and a disordered gauge field. While the spins remain ordered down to T = 0, the gauge field remains disordered in an intermediate regime T c < T < T ∼ K due to thermal fluctuations, which prevent the gauge fields (loop structures) from becoming ordered. Eventually, for T T c ∼ K/100 the gauge structure becomes ordered as well, and the system ends up being in a loop soup, i.e, the Kitaev QSL.
Let us note that the spin-ordering thermal crossover at T and the thermal gauge ordering transition at T c can be identified with two separate peaks in the specific heat and each is accompanied by releasing half of the thermal entropy of the system. While the crossover temperature T has been shown to be captured readily both numerically and experimentally, the thermal transition point T c is beyond the reach of both experiments and our TN simulation. However, we showed convincing proof that our TN simulations work perfectly in the intermediate to high-temperature regimes which is also the most relevant regime for the experiment, as well. The poor convergence of our thermal TN algorithm in the low-temperature regime is typical of all thermal TN algorithm based on purification and is related to the large entanglement growth during real-time evolution close to the critical point and in the ground state of QSL phases [39,89,90].
In order to provide a complete picture about the thermodynamic properties of the Kitaev spin liquid particularly in the low-temperature regime, we calculated the Chern number and the thermal Hall conductivity. We showed that the low-temperature gauge-ordering transition at T c is associated with the non-trivial Majorana band topology, yielding a chiral edge mode that has a non-vanishing Chern number for T < T c . We further saw that the changes in the topological nature of the system below T c is related to the presence of Weyl nodes, which is manifested via nontrivial transport features, i.e., a finite thermal Hall conductivity that is non-zero in the QSL phase (T < T c ) and vanishes above the low-temperature gauge-ordering thermal transition.
Away from the pure Kitaev point, we studied the phase diagram of the 3D Kitaev-Heisenberg model on the hyperoctagon lattice at zero and finite-temperatures. and captured the magnetically ordered phases in the vicinity of the QSL regions. Our simulations also confirme that, in contrast to the QSL phase which has a doublepeak feature in the specific heat, the single thermal phase transition in magnetically ordered phases occurs in the context of Landau symmetry-breaking.
The tensor network techniques introduced in this study can be used as efficient tools for a deeper understanding of the thermodynamics of strongly correlated systems on any dimension and lattice geometry. We, therefore, believe that they are of potential interest for benchmarking both numerical and experimental studies, and will become essential for the discovery of new phase of matter.