Universal Thermodynamics in the Kitaev Fractional Liquid

In the Kitaev honeycomb model, the quantum spin fractionalizes into itinerant Majorana and gauge flux spontaneously upon cooling, leading to rich experimental ramifications at finite temperature and an upsurge of research interest. In this work, we employ the exponential tensor renormalization group approach to explore the Kitaev model under various perturbations, including the external fields, Heisenberg, and the off-diagonal couplings that are common in the Kitaev materials. Through large-scale manybody calculations, we find a Kitaev fractional liquid at intermediate temperature that is robust against perturbations. The fractional liquid exhibits universal thermodynamic behaviors, including the fractional thermal entropy, metallic specific heat, and an intermediate-temperature Curie law of magnetic susceptibility. The emergent universal susceptibility behavior, with a modified Curie constant, can be ascribed to the strongly fluctuating $\mathbb{Z}_2$ fluxes as well as the extremely short-ranged and bond-directional spin correlations. With this insight, we revisit the susceptibility measurements of Na$_2$IrO$_3$ and $\alpha$-RuCl$_3$, and find evident signatures of finite-temperature fractionalization and ferromagnetic Kitaev couplings. Moreover, the peculiar spin correlation in the fractional liquid corresponds to a stripy structure factor which rotates in the extended Brillouin zone as the spin component changes. Therefore, our findings encourage future experimental exploration of fractional liquid in the Kitaev materials by thermodynamic measurements and spin-resolved structure factor probes.


I. INTRODUCTION
Quantum spin liquids realize a novel class of quantum states of matter [1][2][3] which represent the nonmagnetic states with emergent long-range entanglement and fractionalized excitations in frustrated magnets [4,5]. After extensive search for decades, the exactly soluble Kitaev model on the honeycomb lattice [6], which has three types of bond-directional Ising couplings in the trivalent, was found as a prototypical system that can realize both gapped and gapless spin liquid states with spins fractionalized into itinerant and localized Majorana fermions [6,7]. In particular, in the gapless spin liquid phase, time-reversal-symmetry breaking perturbations such as a magnetic field can immediately induce a gapped non-Abelian chiral spin liquid, which has potential application for topological quantum computation [8].
To pursue the materialization of the intriguing Kitaev spin liquid (KSL) states, a race has been triggered for searching crystalline Kitaev-like materials [9][10][11]. Fortunately, there are Mott insulators with strong spin-orbit coupling and effective pseudospin J eff = 1/2 on the honeycomb lattice that can realize the highly anisotropic couplings and have been intensively studied in recent years, including the iridate family such as A 2 IrO 3 (A = Li and Na) [12][13][14][15][16][17][18] and H 3 LiIr 2 O 6 [19], as well as the ruthenate family such as α-RuCl 3 [20][21][22][23][24]. Nonetheless, due to the couplings beyond the Kitaev model such as the Heisenberg, off-diagonal Gamma, and further-neighbor inter- Fig. 1. The cylindrical honeycomb lattice of size YC4 × 6 × 2 is illustrated with the snake-like quasi-1D mapping path (underlying grey lines), adopted in the XTRG calculations. Three types of Kitaev bonds K x , K y and K z are marked in blue, green, and red colors, respectively. The Z 2 flux W P = σ z 1 σ x 2 σ y 3 σ z 4 σ x 5 σ y 6 = +1 (hexagon in white) or = −1 (purple), with σ γ i = 2S γ i , can be pair-flipped by applying a magnetic field h // [0 0 1] in spin-space cubic axis. Typical nearest and next-nearest bonds, d x , d z , d x+y , and d x+z , are indicated by the dashed lines. actions in real materials, semi-classical magnetic orders appear at low temperature (T ) in most of these KSL candidate compounds. Surprisingly, signatures of fractionalization in quantum spins can nevertheless be detected at intermediate temperature in some Kitaev materials [17,[22][23][24], even when the low-T states are classically ordered. These exciting observations and the spin liquid phase found in the presence of magnetic fields strongly suggest that these materials should be in proximity to the KSL phase, and imminently call for systematic understanding on the thermodynamic properties of these Kitaev materials.
For the pure Kitaev model, the finite-T fractionalization of spins has been well-established via the Majorana-based quantum Monte Carlo (QMC) simulations [25][26][27], which found two temperature scales and an intermediate regime with unconventional properties [25][26][27][28][29]. This may shed light on the understanding of the exotic equilibrium [24] and transport measurements [30] in the Kitaev materials. Nonetheless, the gap between experiments and theoretical understanding is still apparent. Since the non-Kitaev interactions and magnetic field [31,32] spoil the exact solution and restrict QMC simulations within relatively high temperature [33], little is known on the fractionalization of spins in the perturbed Kitaev models.
In this work, we employ the exponential tensor renormalization group (XTRG) [34,35] at finite temperature to explore the Kitaev model with various additional perturbations. We find that the well-recognized characteristic features at intermediate temperature, including the fractional thermal-entropy plateau and the linear-T specific heat scaling clearly persist. More interestingly, we find an emergent universal Curie law of magnetic susceptibility below a high-T scale T H (and above the low-T scale T L ), due to the fluctuating gauge fluxes and the predominant bond-directional nearest-neighboring (NN) spin correlations. Through exploring these robust finite-T features, we find a generic scenario in the perturbed Kitaev model: A fractionalization of the quantum spin takes place spontaneously as temperature decreases, and half of the spin degrees of freedom, i.e. the itinerant Majorana fermions, release large part of the entropy below T H , accompanied by the saturated short-range spin correlation; meanwhile the other half, the gauge fluxes, fluctuate strongly until near the much lower temperature T L . In between the two temperature scales, there reside a unique finite-T quantum state dubbed the Kitaev fractional liquid.
Based on our systematic characterization, we reveal that the fractional liquid regime does exist at intermediate temperature, for an extended range of external fields or non-Kitaev terms including the Heisenberg or Gamma couplings, even when the KSL ground state is altered. These results therefore pave the way towards filling the gap between theoretical and experimental studies of the Kitaev physics. In particular, with insight from the universal susceptibility scaling, we revisit the magnetic susceptibility data of Na 2 IrO 3 and α-RuCl 3 and find evident signature of Curie susceptibility scaling in the two prominent Kitaev materials, suggesting the existence of the fractional liquid therein. We also propose the spincorrelation diagnostics of the fractional quantum liquid, which is available through spin-resolved neutron or X-ray diffusion scatterings.
The rest part of the paper is organized as follows. In Sec. II, we introduce the Kitaev model with perturbative couplings and the XTRG method for finite-T simulations. Secs. III, IV, and V present our main results on the universal thermodynamic properties and spin-correlation characterization of the intermediate-T Kitaev fractional liquid that is robust against additional interactions. In Sec. VI, we discuss the relevance of universal susceptibility scaling in two Kitaev materials, Na 2 IrO 3 and α-RuCl 3 , and propose the spin structure diagno-sis of the fractional liquid. Finally, we summarize in Sec. VII our main conclusions and discuss their implications in future studies of the Kitaev model and materials.

II. MODEL AND METHODS
We study the celebrated Kitaev honeycomb model with various perturbative terms as follows where S i = {S x i , S y i , S z i } represent the spin-1/2 vector operator at site i. The three cubic coordinates {α, β, γ} represent {x, y, z} up to a cyclic permutation, and i, j γ denotes the NN γ-bond [see Fig. 1] with an Ising coupling K γ . The non-Kitaev perturbations include the NN isotropic Heisenberg coupling J and an off-diagonal Gamma (Γ) interaction. Moreover, we apply an external field h coupled to the S z component, which is along [0 0 1] spin-space direction and tilts an angle of 45 • from the honeycomb ab plane (see Fig. 1) in typical Kitaev materials.
For the pure Kitaev model, spin can be represented as localized (b γ i ) and itinerant Majorana fermions (c i ) through the parton construction S γ i = i 2 b γ i c i [8]. The itinerant Majorana can be related to the building up of spin correlation S z i S z j = u z i j c i c j , and u γ i j = b γ i b γ j = ±i represents a Z 2 gauge connection. The product of u γ i j on the six edges of a hexagon gives rise to the flux W P = ±1 as illustrated in Fig. 1, which conserves in the pure Kitaev case.
For the typical Kitaev materials Na 2 IrO 3 and α-RuCl 3 , the predominant Kitaev interactions have been widely accepted to be ferromagnetic (FM) [10,24,33,[36][37][38]. Thus, we choose K γ < 0 henceforth and set K z = −1 as the fixed energy scale. The Gamma term is suggested to be of next-to-leading order and with Γ > 0 in α-RuCl 3 [39,40]. On the other hand, in Na 2 IrO 3 the Heisenberg coupling is believed to play an important role and is proposed to be antiferromagnetic (AF) [36,41]. Thus, in our simulations below we introduce additional AF Heisenberg or Gamma term with J, Γ > 0, and focus on the finite-T quantum states of the perturbed Kitaev model.
In the XTRG simulations [34,35], we follow the snake-like path and map the cylinder into a quasi-one-dimensional (1D) system as illustrated by the grey lines in the shadow of Fig. 1. We work on the cylindrical lattice YC W × L × 2 (see Fig. 1), with width W (up to 4) and length L (up to 6), thus total site N = 2WL, which reproduces results in excellent agreement with benchmark exact diagonalization (ED) on small system size as well as large-scale Majorana QMC simulations (see Appendix A). From comprehensive benchmark results, we see finite-size effects are in most cases modest, at least for the intermediate-T regime of main focus in the current study. This is in contrast from the ground-state calculations and can be ascribed to the extremely short-range correlations at intermediate temperature, even under perturbations (see Appendix B).
Moreover, as there are two characteristic temperature scales in the Kitaev model that are typically separated by more than one orders of magnitude, our XTRG approach that cools down the system exponentially in temperature scale turns out to be very efficient to address this problem. Therefore, we can exploit the XTRG method that is free of sign problem in simulating the Kitaev model under various perturbations down to extremely low temperature, T 0.005. In practical calculations we retain bond dimension up to D = 600 that is capable to describe the thermal quantum states accurately and produce well converged thermodynamic results (Appendix C). We have also employed the density matrix renormalization group (DMRG) method to simulate the ground state, confirming the XTRG results at very low temperature.

III. KITAEV FRACTIONAL LIQUID
A. Ferromagnetic Kitaev model at finite temperature It has been well-established that the itinerant and localized Majorana fermions can give rise to the two temperature scales, T H and T L , respectively. At intermediate temperature T L T T H , an exotic fractional quantum liquid emerges. To start with, we simulate finite-T properties of the isotropic FM Kitaev model (K x,y,z = −1) by XTRG on a YC4 × 6 × 2 cylinder. In Fig. 2(a), we plot the spin correlation S z i S z j measured along the NN bonds d z , d x (d y always leads to the same  results as d x , thus not shown) and next-nearest ones d z+y and d x+y (as depicted in Fig. 1). From the zero-field XTRG results, we find only the nearest d z -bond correlations is established below T H , and the rest remain small down to the lowest temperature. Thus the spin correlations S γ i S γ j are nonzero only on the NN γ-bond, i.e., they are extremely short-range and bondoriented, in full agreement with the analytical results [42,43].
Matter of fact, the itinerant Majorana fermions release most of their entropy and the NN spin correlation is established at around T H [ Fig. 2(a)]. On the other hand, the rest "half" of spins, Z 2 fluxes, strongly fluctuate until the low-T scale T L , as evidenced by small W P above T L in Fig. 2(c). Correspondingly, in Fig. 2(d) the thermal entropy exhibits a quasi-plateau with S 1 2 ln 2, which constitutes a thermodynamic signature of the intermediate-T Kitaev fractional liquid.

B. Kitaev fractional liquid at finite fields
Recently, the magnetic field effects on the Kitaev materials have raised great interest [30][31][32][33]44]. For example, thermal quantum Hall signature has been detected in the field-induced spin liquid phase of α-RuCl 3 [30]. Given the fractional liquid scenario in pure Kitaev model, it is also interesting to further study how it evolves under external fields [33]. In the following, we mainly focus on the field direction h // [0 0 1], and the effects of external fields along other directions are discussed in Appendix D, where qualitatively the same results are seen.
The finite-field entropy results are contour plotted in Fig. 3(a), with the S vs. T curves shown explicitly in Fig. 2(d).
The two isentropic lines with S / ln 2 = 0.75 and 0.25 in Fig. 3(a) correspond to the crossover temperatures T H and T L (or T L ) very well. The intermediate-T regime, which corresponds to the 1 2 ln 2 entropy plateau and signifies the fractional liquid, extends to a rather wide field range regardless if the low-T phase is an asymptotic KSL or not. The two-step release of entropy gives rise to the double-peaked specific heat curves as illustrated in Figs. 3(b) and 4.
While T H remains stable as the magnetic field h increases, the low-T crossover temperature T L drops at h = h c 0.02, and rises up again as h further increases [T L in Fig. 3(b)]. As a matter of fact, T L represents a different low-T scale since for h > h c the ground state is no longer KSL but a topologically trivial field-induced FM phase. This can be seen evidently in Fig. 3 [45,46], although [111] fields are usually considered.

C. Spin correlation and structure factor
In robust fractional liquid regime, we again observe extremely short-ranged and bond-directional spin correlations that characterize this exotic finite-T quantum state even under external fields. In Fig. 2(b), we plot the spin correlations S z i S z j − M 2 with the field-induced uniform background subtracted. Although the longer-range correlations rise up at lower temperature and spoil the extremely short correlation length "inherited" from the pure Kitaev model, we nevertheless can observe that the d z correlation is predominate over a very wide range of intermediate temperatures.
Correspondingly, we compute the static spin structure factor from the equal-time correlation function where r 0 is the position vector of fixed reference site (i 0 ) and j (and thus vector r j ) runs over the lattice sites. We observe in Fig. 3(c) a stripy structure factor at intermediate temperature for both h = 0 and h 0, resembling that in the KSL ground state [17,47,48]. For h > 0, there exists an additional bright spot at the Γ point, which represents the field-induced FM background in the system and can be observed clearly at intermediate temperature. Besides the two-point correlation, from Figs. 2(c) and 3 we observe that the Z 2 flux operator W P remains small at intermediate temperature, suggesting that the gauge fluxes can still flip "freely" under various fields in the fractional liquid regime.

D. Linear-T specific heat in the itinerant Majorana metal
In the Kitaev fractional liquid, itinerant Majorana couples to the thermally fluctuating fluxes and leads to a finite density of states emerging at the Fermi surface, which can be characterized by a linear T -dependence behavior of the specific heat C V [26]. It is therefore of great interest to check if this intermediate-T Majorana metal state could persist at finite external fields, by examining the specific heat behaviors.
In Fig. 4, we show the XTRG results of C V that indeed exhibit linear behaviors in the lower part of the intermediate-T regime, in excellent agreement to previous QMC result at zero field [26], and thus extends the conclusion there to a finite range of fields. As temperature further decreases, for h h c the C V curve deflects the linear-T behavior slightly above T L (or T L ), and when the fields h h c , the low-T scale T L rises up as h increases, resulting in the reducing of the temperature range showing metallic specific heat. We reveal this behavior more clearly in the inset of Fig. 4, where a plateau appears in C V /T for h up to h c 0.02, and it then turns into a shoulderlike structure for larger fields up to h = 0.045, a reminiscent of the itinerant Majorana metal.

IV. EMERGENT CURIE LAW IN THE FRACTIONAL LIQUID
The linear-T specific heat scaling constitutes an exotic thermodynamics of the Kitaev fractional liquid. However, both itinerant Majorana and Z 2 flux have contributions in the specific heat and thus the Majorana metallic behavior is masked by the entropy release of flux freezing at low temperature near T L . Moreover in experiments, there are contributions other than the Majorana fermions, e.g., magnon excitations on top of the long-range magnetic order and the excess phonon contribution not fully deducted, that may enter the intermediate-T specific heat data and thus mask the expected universal linear-  T behavior. Therefore, it calls for in-depth analysis of the other readily accessible thermodynamic measurement in experiments that reflects more of the intrinsic magnetic properties of the system, i.e., the susceptibility χ. From the analysis, we reveal that there emerges an universal Curie-law behavior in the intermediate-T susceptibility that is robust against perturbations and origins from the fractionalization of spins.
A. Magnetic susceptibility in the isotropic and anisotropic Kitaev models As a direct characterization of unconventional paramagnetism of the fractional liquid, in Fig. 5(a), we show the susceptibility of the isotropic FM Kitaev model under various magnetic fields. For small fields h h c , the susceptibility data collapse into a universal Curie form χ = C T in the fractional liquid regime, further consolidated by plotting 1/χ vs. T in the inset of Fig. 5(a). This unexpected universal Curie behavior extends across the fractional liquid regime, spanning a very wide temperature range of about one order of magnitude (T ∼ 0.03-0.3), constitutes an emergent phenomena distinct from conventional paramagnetism at high temperature. The latter corresponds to a Curie constant C = S (S + 1)/3 = 1/4 relating to the intrinsic moment of spin S = 1/2, while the emergent Curie constant C ≈ 1/3 [see Fig. 5(a)] in the fractional liquid regime, suggesting the presence of spin correlations that renormalizes the effective magnetic moments.
Besides the isotropic case, we further turn to the anisotropic Kitaev model and also find emergent Curie behaviors in Fig. 5(b,c) that can extend down to even lower temperatures. The Curie-Weiss fittings give C ≈ 0.50 for K x,y = −0.1 and −0.25, whose susceptibility is already very close to the ensemble of isolated FM spin pairs [red lines in Fig. 5(b, c)], and the Curie constant is close to twice a single spin-1/2, i.e., C = 1/2. The emergent Curie susceptibility there can extend down to temperatures significantly lower than T L (of the isotropic case), due to the reduction of flux gap. Therefore, starting from the trivial K x,y = 0 limit, as we turn on the FM coupling K x,y in the system, the Curie behavior remains, although the temperature range gets narrowed and the Curie constant becomes modified. For example, with K x,y = −0.5, there also exists a Curie-law regime with C reduced to 0.42, as indicated by the blue dashed lines in Fig. 5(b, c). This universal Curie scaling extends all the way to the isotropic point, implying spins can be pair-flipped nearly freely in the intermediate-T regime.
Matter of fact, the emergent Curie behaviors can be observed as long as the spin correlations S z i S z j are mostly restricted within the d z -bond [see Fig. 2(a,b)] and in a quasiplateau vs. T at intermediate temperature shown in Fig. 5(d). We note that, according to the perturbation theory, the Z 2 fluxes remain approximately conserved, and spin correlations beyond the d z -bonds only appear at the h 2 c order. At the same time, Z 2 fluxes should be flippable in a nearly free fashion, which can be observed as W P 1 in Fig. 2(c). These conditions are satisfied even when perturbative external fields are applied, that gives rise to the intermediate-T Curie law and will be elaborated in Sec. IV B below.

B. Fluctuating fluxes and diverging susceptibility
The diverging susceptibility reflects the strongly fluctuating fluxes in the fractional liquid below T H , and deviation occurs when the flux starts to freeze out at about T L . To see this, we exploit the Kubo formula of susceptibility and assume M = 0 in the fractional liquid for the sake of discussion (Appendix E), where S z i 0 (τ) = e τH S z i 0 e −τH , and j runs over NN sites of i 0 by d z bond (as well as i 0 itself) in the fractional liquid regime due to the extremely short-range correlations. As illustrated in Fig. 1, the spin operator S z j flips a pair of fluxes on two neighboring hexagons (color in purple) and changes the flux configuration of energy eigenstate. After evolving along the imaginary time of τ, S z i 0 can flip the Z 2 fluxes back and thus restore the flux configurations, which results in a nonzero correlation.
We denote the energy eigenstates by |E n {W P } , where {W P } represents the flux configurations, and n labels the individual states within the {W P } sector. By inserting the orthonormal basis, we arrive at the Lehmann spectral representation where |E n {W P } represents a state n in the flux-flipped sector {W P }. Since the itinerant Majorana fermions are weakly coupled to the flux at intermediate temperature [26] and thus the excitation energy is close to the flux gap, ∆ n,{W P };n , Consequently, S z i 0 (τ)S z j β is virtually τ-independent at intermediate temperature and can be approximated as S z i 0 S z j β in Eq. (3), thus recovering the conventional susceptibility for- which holds only for T T L . As observed in Fig. 2(a) and Fig. 5(d), the correlation S z i 0 S z j β exhibits a quasi-plateau vs. temperatures between T L and T H , resulting in a 1/T -divergent susceptibility that accounts for the emergent Curie law in the fractional liquid regime. On the other hand, in the KSL regime at T T L , albeit the correlation remains extremely short there, the (flux) excitation energy ∆ n,{W P };n ,{W P } ∼ T L T leads to an exponentially suppressive factor in τ, and thus a non-divergent susceptibility χ below T L in Fig. 5.
Above we have discussed the emergent susceptibility scaling in the FM Kitaev model, and below we briefly discuss the case of AF model, and refer the details to Appendix F. For the AF Kitaev model, the spin correlations are also extremely short-ranged, whose absolute values are identical to those of the FM case, but with opposite sign (AF correlation). Therefore, the emergent Curie law is absent in the uniform susceptibility χ of the AF Kitaev model, while it appears instead in the stagger susceptibility. This can be understood intuitively as follows: one can flip simultaneously the spins on one of the two sublattices on the honeycomb lattice, tuning the AF Kitaev model into an FM one, and the staggered susceptibility χ s corresponds to the uniform χ exhibiting universal scaling at intermediate temperature. Such an interesting correspondence between the FM Kitaev model under uniform field and the AF model under staggered field has been noticed in other studies under different circumstances [49]. Our results indicate a sensitive way to determine the sign of Kitaev coupling via analyzing the susceptibility measurements.

V. FRACTIONAL LIQUID UNDER NON-KITAEV PERTURBATIONS
The Heisenberg coupling is proposed to be relevant in the Kitaev materials Na 2 IrO 3 [13,15], and the Kitaev-Heisenberg model has thus been extensively studied for both the groundstate [13,16,45,47,50] and finite-T properties [51,52]. On the other hand, in α-RuCl 3 there exist a symmetric off-diagonal exchange -the so called Gamma term [39,46,48,53] -which turns out to be important non-Kitaev interactions. In this section, we explore the Kitaev-Heisenberg (on YC4 × 6 × 2 cylinders) and Kitaev-Gamma (YC4 × 4 × 2) models [See Eq. (1)] at finite temperature, from which we see that although the low-temerapture fate of the Kitaev model may be drastically affected by the Heisenberg and Gamma terms, the behavior at high and intermediate temperatures, especially in the Kitaev fractional liquid regime, remain qualitatively unchanged.

A. Fractional thermal entropy and specific heat hehaviors
We start with the thermal entropy and specific heat calculations of the Kitaev-Heisenberg and Kitave-Gamma models, and show the results in Figs. 6 and 7. Between the two entropy-release steps [ Fig. 7(c,d)], there exists an extended regime with thermal entropy S 1 2 ln 2 [ Fig. 6(a,b)] for both models, which indicate the existence of Kitaev fractional liquid at intermediate temperature. The 1/2-entropy quasi-plateau remains to be evident till Heisenberg coupling J 0.16 and for Gamma perturbation of Γ 0.06. Beyond this Γ or J range, a shoulder-like feature with slow entropy release can still be recognized in both models. In Fig. 6(a,b), the crossover temperature scales, T H and T L (T L ), can be determined from the isentropic lines of S / ln 2 = 0.75 and 0.25, which coincide with the higher and lower crossover temperatures in the double-peaked specific heat curves [see Fig. 7(a,b) and also Appendix G].
In Fig. 7(a,b), we show that the linear-T behavior of specific heat also survives in the presence of a finite J or Γ perturbation, which can be also recognized from the C V /T plots in the insets. This universal specific heat behavior suggests the emergence of Majorana-metal states in the fractional liquid regime. Nevertheless, similar to what was observed in Fig. 4, the J and Γ perturbations can shift the flux excitations as well as other low-energy spin fluctuation modes towards higher temperatures and thus may mask the linear-T C V , leaving eventually only a shoulder-like reminiscent in the insets of Fig. 7(a,b).
The finite-T properties can sensitively reflect the quantum phase transition in the ground state. In Fig. 6(a), T L (T L ) is non-monotonic vs. J: For J J c 0.12, T L represents a crossover from the KSL to the fractional liquid; while for J > J c , the stripe anti-ferromagnetic (SAF) order sets in below T L , where correlations along d x , d x+y , etc, also increase rapidly as shown in Fig. 6(c). Therefore, we infer from the finite-T results that the Kitaev-Heisenberg model undergoes a transition from the KSL to SAF phase at around J c 0.12, consistent with previous ground-state studies [13,50,53].

B. Bond-directional spin correlation and stripy structure factor
From the spin-correlation data shown in Fig. 6(c,d), we emphasize that the two models share common features in the Ki-  taev fractional liquid regime. First, the spin correlation remains extremely short-ranged, and mostly present on the d z bonds. This is consistent with a perturbation-theory analysis: except on the d z bonds, spin correlations are suppressed because Z 2 fluxes are approximately conserved. Furthermore, such correlations appear at the first order in J but only at the second order in Γ, consistent with the fact that J has a larger effect on spin correlations than Γ, as shown in Fig. 6(c,d). This prominent spin-correlation feature consequently leads to the stripy structure factor in the insets of Fig. 6(a,b).
Under the Heisenberg perturbations, we observe in the inset of Fig. 6(a) a clear structure peak at the X point, which indicates the enhanced short-range SAF correlation, with the classical spin pattern shown rightside. Besides the SAF peak, we notice that there exists a stripy structure background which is related to the predominate bond-directional NN correlations at intermediate temperature. On the other hand in Fig. 6(b), for the Kitaev-Gamma model, the low temperature scale T L only slightly changes. The KSL-like stripy spin structure can be observed with no clear signature of spin ordering, which, together with the nonzero fluxes, W P 0, at low temperature in Fig. 6(f), suggest that the low-T quantum state remains in the asymptotic KSL phase at least for Γ 0.1, in agreement with recent ground-state studies of the Kitaev-Gamma model [46,53].

C. Universal susceptibility scaling
Now we switch to the analysis of the magnetic susceptibility scalings in both models. At intermediate temperature, we again observe bond-directional spin correlations that are mostly restricted within the NN sites [see Fig. 6(c,d) and discussions in Sec. V B], as well as the strongly fluctuating flux [ Fig. 6(e,f)] under an extended range of J or Γ perturbations. Following the line of discussions in the pure Kitaev model, we expect an emergent Curie law of susceptibility in the fractional liquid, and this is indeed observed in Fig. 6(g,h). Although the temperature window of the emergent Curie behaviors gets narrowed as J or Γ increases, the universal susceptibility scaling can be very clearly identified across the intermediate-T regime, regardless of whether the ground state is a KSL or semiclassically ordered.
Nevertheless, there are still some difference in the susceptibility scalings between the Heisenberg and Gamma perturbations in the susceptibility scalings. At the temperature above T L , the susceptibility curves in the Kitaev-Gamma model collapse into the same universal behavior, with C ≈ 1/3 (thus C /C ≈ 1.33, the same to the pure Kitaev model) and with the same small θ . This is in contrast to the case of the Kitaev-Heisenberg model, where although C 1/3 remains the value of the pure Kitaev model, the Curie-Weiss temperature θ grows as J increases, which suggests the fluxes are no longer "free" but coupled with each other since the Heisenberg perturbation J introduces longer-range correlations to the system [ Fig. 6(c)] . Accordingly in Fig. 6(e), the gauge flux for J J c behaves very similarly to that of the pure Kitaev model, i.e., W P freeze at around T L ; while it gets suppressed at low temperature for J > J c .
Overall, we observe that the emergent Curie law of susceptibility scaling appears to be robust against the perturbations at intermediate temperature, which is found to be suppressed only when the finite flux gap or the longer-range correlation induced by the non-Kitaev terms takes effects at the low temperature scale. Therefore, such universal susceptibility scaling can serve as a useful tool probing the fractional liquid in the Kitaev materials.

A. Kitaev paramagnetism in the Kitaev materials
Thermodynamic features, including the double-peaked specific heat curve, linear-T behavior, and the fractional thermal entropy S integrated from C V /T have been exploited in detecting finite-T fractionalization in the Kitaev materials [18,24,[54][55][56]. However, the magnetic response to external fields -the susceptibility χ -has been much less analyzed, except for the high-T behaviors in the Kitaev materials. With insight from the emergent intermediate-T Curie susceptibility from the above model studies, we revisit the experimental χ data of two prominent Kitaev materials, i.e., Na 2 IrO 3 [14] and α-RuCl 3 [54], which are believed to realize the Kitaev model under the Heisenberg and/or Gamma perturbations.
In Fig. 8, we analyze the in-plane and out-of-plane susceptibility data of both materials, where the emergent Curie behaviors can be observed right below certain high temperature scales, i.e., T Ir  higher crossover temperatures in the magnetic specific heat curves of two materials. For example, in Na 2 IrO 3 the higher specific heat peak locates at around 110-120 K [17,18] which is in very good agreement with the T Ir H in Fig. 8(a). Similarly, the high-T entropy rapidly releases at around 100-110 K in α-RuCl 3 [24,54,56], also coinciding with T Ru H in Fig. 8(b). Through Curie-Weiss fittings of the intermediate-T susceptibility results in Fig. 8, we find the the ratios between intermediate-and high-T Curie constants C /C ∼ 1.3 in both materials (and either within the ab-plane or parallel to the c axis), suggesting the existence of the Kitaev fractional liquid in these two Kitaev materials. The emergent diverging susceptibility at intermediate temperature is very sensitive to the sign of Kitaev interaction K, and for the AF case (K > 0) there exhibits a convex uniform susceptibility curve without any emergent Curie behavior at intermediate temperature (Appendix F). Therefore, the analysis in Fig. 8 suggests the FM Kitaev interactions in both compounds Na 2 IrO 3 and α-RuCl 3 , in agreement with previous theoretical [10,36,57,58] and experimental [38] estimates.
Meanwhile, in Fig. 8 we observe the intermediate Curie-Weiss behaviors are deviated at the lower temperature T Ir L 50-55 K and T Ru L 60-65 K. Therefore, in these two materials there are fractional liquid regimes stretching appreciable temperature ranges as wide as 65-70 K (Na 2 IrO 3 ) and 45-50 K (α-RuCl 3 ), respectively. In addition, by comparing Fig. 8(a) to Fig. 6(g), we find a relatively wider intermediate-T range for the emergent Curie-Weiss in Na 2 IrO 3 . This resembles that of the Kitaev-Heisenberg model, where universal scaling is deviated at around the susceptibility peak temperature. This observation suggests that indeed the Heisenberg-type perturbation can play an important role in this material [13]. On the other hand, in Fig. 6(g) the Gamma term deflect the Curie-Weiss behavior way before the susceptibility peaks, and this bears a resemblance to the case of α-RuCl 3 shown in Fig. 8(b), supporting thus the proposal that the Γ term (together with other perturbations) is particularly relevant in this material [40,59].
Overall, through analysis of the magnetic susceptibility data of Na 2 IrO 3 and α-RuCl 3 , we find a universal susceptibility scaling that may be ubiquitous in other members of the intriguing family of Kitaev materials and can serve as a sensitive thermodynamic signature of the Kitaev fractional liquid. Besides, the susceptibility scaling indicates the existence of ferromagnetic Kitaev interactions, together with other non-Kitaev couplings, in the two prominent Kitaev materials.

B. Spin-correlation diagnostics of the Kitaev fractional liquid
The short-range correlated Kitaev fractional liquid leaves unique features in the spin correlations and thus structure factors. At intermediate temperature, the nonzero correlation S γ i S γ j is restricted within the NN γ-bond, with the rest correlations negligibly small, even under magnetic fields and non-Kitaev perturbations, as shown in Figs. 2(a,b) and 6(c,d). Correspondingly, in Figs. 3(c) and 6(a,b) we see this peculiar bond directional correlation pattern results in a stripy-like background in the z-direction spin structure factors S z (q).
Matter of fact, this stripy spin structure also appears in the structure factors S γ (q) with γ = x, y, and the direction is "locked" with the associated spin components γ. As shown in Fig. 9(a,e), the S x (q) stripe extends from the left-top to the right-bottom, which rotates counterclockwise by 120 • when switched to S y (q) [ Fig. 9(b,f)], followed by another 120 • rotation to S z (q) in Fig. 9(c,g). Note that the spin structure factors results in Fig. 9 are evaluated at an intermediate temperature T 0.15 for the pure Kitaev model (the upper row) and Kitaev-Heisenberg model (lower). In the latter Kitaev-Heisenberg case with J = 0.14, there emerges a structure factor peak representing the semicalssical SAF correlation, which also rotate as the spin component switches. Such a rotation of short-range correlation peak has previously been observed in the spin-resolved magnetic X-ray measurements right above the ordering transition temperature (T 17 K) [17], which provides a direct evidence of the dominant bonddirectional interactions in the Kitaev material Na 2 IrO 3 . Therefore, our findings in Fig. 9(e-g) confirm that this experimental observation is indeed related to the presence of Kitaev interactions in the material, although the SAF order involved here is different from the zigzag type in Na 2 IrO 3 .
Moreover, our results in Fig. 9(d,h) also shed light on the renowned star-like spin structure factor observed in inelastic neutron scattering measurement of α-RuCl 3 [23,24]. From the full structure factor S (q) = γ S γ (q), we observe in the pure Kitaev model that the superposition of the three stripes leads to an approximately rotationally symmetric shape in Fig. 9(d). By adding bright spots representing short-range SAF correlations on top of that, we instead see a six-fold symmetric S (q) in Fig. 9(h). This clearly suggests that the star-shape structure factor (with six-fold symmetry) observed in inelastic neutron scatterings on α-RuCl 3 can be ascribed to the coexistence of the Kitaev bond-directional correlation and the short-range zigzag order. As the temperature/energy scale gradually increases, it approaches asymptotically the rotationally symmetric spin structure in Fig. 9(d), which has also been observed in experiments.
In general, we propose that the bond-directional spin correlation and related intermediate-T stripy feature provide a useful tool for diagnosing the Kitaev materials. Specifically, diffuse X-ray (for Na 2 IrO 3 ) or spin-resolved neutron (for α-RuCl 3 ) measurement at intermediate temperature, say, T = 50-100 K, can be exploited to observe the stripy spin structures S γ (q) that serves as a direct evidence of the Kitaev fractional liquid.

VII. DISCUSSIONS AND OUTLOOK
In conventional spin models and magnetic materials, a "crystallization" of spins occurs spontaneously, from the high-T gas-like paramagnetic to the low-T phase with long-range order. Interestingly, in the frustrated Kitaev model and related Kitaev materials there exists an alternative scenario, where a fractionalization of spins takes place upon cooling. This gives rise to an intermediate-T Kitaev fractional liquid and eventually to the low-T quantum spin liquid phase with long-range entanglement.
Exploiting the XTRG method in the finite-T simulations, aided by the ground-state DMRG calculations, we reveal unambiguously a Kitaev fractional liquid in the Kitaev model under perturbations, including the magnetic field, Heisenberg and Gamma terms. In the exotic finite-T quantum states, there exist fractional thermal entropy, linear-T specific heat, and, in particular, an emergent Curie law in magnetic susceptibility. These thermodynamic features are robust against perturbations, regardless if the ground state is still in a spin liquid phase or not.
Remarkably, this universal Curie susceptibility that reflects the thermally free flux degrees of freedom, together with the previously recognized metallic specific heat behaviors of itinerant Majorana, each revealing one aspect of the spin fractionalization, constitute a comprehensive thermodynamics characterization of the exotic Kitaev fractional liquid. In light of this, and from analyzing the experimental susceptibility data, we find signatures of the Kitaev fractionalization liquid in two Kitaev materials Na 2 IrO 3 and α-RuCl 3 , which support ferromagnetic Kitaev interactions in the two materials.
Matter of factor, the diverging susceptibility is a robust thermodynamic signature also related to the extremely shortranged and bond-directional spin correlations in the fractional liquid regime, which corresponds to a stripy spin structure factor. Therefore, besides with the magnetic susceptibility measurement, our findings also encourage further material exploration of the Kitaev fractional liquid using spin-resolving neutron and resonant X-ray scatterings. Lastly, our work paves the way for the precise determination of interaction parameters in the Kitaev materials by fitting the measured thermodynamic data, which we leave in a future study.   Comparisons between the large-scale XTRG and QMC simulations [27] are presented for (c) the thermal entropy S / ln 2 (orange) and flux W P (blue), (d) specific heat C V , and (e) spin correlations S z i S z j . We also compare the XTRG susceptibility results to the QMC data on the L = 20 torus [29], where the high-and intermediate-T fittings are also performed to the latter (dashed lines). We show in (c-f) both the YC4 × 6 × 2 results (solid curves) and the YC4 data (dotted curves) obtained by subtracting length L = 6 by L = 4 results. XTRG results of susceptibility are computed under [0 0 1] field, which agree excellently with the QMC results in (f) above the low temperature T L , although the latter is along the [1 1 1] direction.
Through this expansion, H n can also be expressed as an MPO using tensor compression techniques [61], followed by a sumand-compression procedure that finally leads to the compact MPO form of ρ(τ).
Following the XTRG idea of exponential cooling procedure, we square the density matrix and obtain at the n-th step ρ(2 n+1 τ) = ρ(2 n τ) · ρ(2 n τ), with which we can compute the free energy at β ≡ 2 n+2 τ via the thermo-field double formalism by a bilayer tensor trace [62]. Given the free energy, other thermodynamic quantities including the specific heat, thermal entropy, and magnetic susceptibility, etc, can also be computed.
We first benchmark the XTRG results with exact diagonalization (ED) on a 12-site cluster (with the same width and boundary condition as YC4 geometry but has only three sites along the zigzag edge). In Fig. A1(a), the relative errors of free energy δ f / f reduces as the retained bond dimension D increases, reaching a very high accuracy with D = 600. On top of that, we show the XTRG results of specific heat computed with D = 600, of pure Kitaev, Kitaev-Heisenberg, and Kitaev-Gamma models, and show them in Fig. A1(b), where we also find excellent agreement with ED results.
As shown in Fig. A1(c-f), we compute the thermal entropy, flux expectation, specific heat, spin correlations and magnetic susceptibility, and compare our XTRG results with the QMC data taken from Refs. 27 and 29. Despite the very different system sizes, i.e., YC4 cylinder vs. L = 20 torus, excellent agreement between two methods can be observed at high and intermediate temperatures, covering the fractional liquid regime. The modest deviation of XTRG from QMC data near the lower temperature scale T L can be ascribed to the cylindrical boundary condition and the finite width, which nevertheless does not lead to any qualitative difference for our discussion.
Moreover, we have also analyzed the large-scale QMC susceptibility data in Fig. A1(f), and from the fitting we find again an emergent Curie behavior clearly at intermediate temperature. The Curie constant C 1/3 is also very close to the YC4 XTRG results in Fig. 5(a) of the main text.

Appendix C: Matrix Product Operator Entanglement
In the MPO form of density matrix ρ(β ≡ 1/T ), we compute the bipartite entanglement entropy S E (β) = −Tr(ρ A lnρ A ), whereρ A is the reduced matrix of subsystem A (chosen as half the system) in the superstate |ρ(β) purified from the thermal density matrix ρ(β). In XTRG, S E (β) can be used to measure the required computational resource, since the bond dimension scales as D ∼ e S E (β) for an accurate simulation down to temperature T ≡ 1/β. For a gapped system, S E typically saturated as T is below the excitation gap ∆, and thus one can in principle simulate the low-T properties even down to the ground state with a sufficiently large, but finite, bond dimension D. The logarithmic entanglement scaling S E (β) ∼ η ln β requires D ∝ β η with η an algebraic exponent. The logarithmic entanglement scaling at low temperature usually corresponds to gapless excitations. Some typical examples include the low-T thermal states near the 1+1D conformal critical point and in 2D Heisenberg model with Goldstone modes [34,35]. In the former case, the exponent has been found to be universal, η = c/3, with c the central charge of the conformal field theory [34].
In Fig. A3(a), we show the bipartite entanglement entropies S E across the center of the system under various magnetic fields h, which exhibits a logarithmic scaling as S E ∼ 2 3 ln β in the intermediate-T regime, as indicated by the black dashed line. In Fig. A3(b,c), we have also computed the MPO entanglement for the Kitaev-Heisenberg and Kitaev-Gamma models, where similar logarithmic behaviors are also observed in the fractional liquid regime. As the perturbative strength increases, for all three cases in Fig. A3(a-c), the entanglement S E deviates the logarithmic behavior, exhibits a shoulder/peak and eventually converges to a constant at sufficiently low temperature.
The exponent η ≈ 2/3 suggests that the required computational resource for simulating the Kitaev models at intermediate temperature is similar to a conformal critical point with central charge c = 2. Besides, in Fig. A3(a) the maximal values of S E appears at the lowest temperature, and is less than 3.5, which reveal the accessibility of the YC4 × 6 × 2 systems by retaining not too many bond states. Matter of fact, we have plotted in Fig. A3(a)   In Fig. A4, we show the χ data along various field directions, as well as the case h // [0 0 1] considered in the main text. We label the field direction lĥ x + mĥ y + nĥ z by [l, m, n], whereĥ x ,ĥ y ,ĥ z are three cubic axes in the spin space, therefore the Zeeman term reads Firstly, we apply an in-plane field h // [-1 0 1] whose direction is shown in Fig. A4(a), and compute the susceptibility. The results are shown in Fig. A4(b,c), where an intermediate Curie behavior can be clearly identified. For the out-of-plane fields h perpendicular to the ab plane, i.e., h // [1 1 1], there also emerges a Curie susceptibility [see Fig. A1(f)], as already discussed in Appendix A. As for more fields directions, we perform ED calculations on systems of very limited size, N = 12 cluster, and show the results in Fig. A4(d). A diverging susceptibility behavior near the high-T scale T H is also evident, confirming that the emergent Curie-law susceptibility is robust and ubiquitous in the Kitaev model. This universal susceptibility scaling can be observed under various fields of different directions, on different lattice sizes/geometries, and by different methods including ED, QMC, and XTRG approaches.

Appendix E: Kubo Formula of Magnetic Susceptibility
Here we derive the Kubo formula of the magnetic susceptibility, in terms of imaginary-time dynamical correlations. Suppose the magnetization of the system is uniform, the magnetic moment per spin, under a perturbative magnetic field h along z-axis, is the expectation value of S z on arbitrary site i 0 where H(K, h) = H 0 (K) − hS z tot , H 0 (K) is the pure Kitaev Hamiltonian, and S z tot = j S z j is the total magnetization operator.
The magnetic susceptibility χ = ∂M/∂h can be obtained by taking derivative of Eq. (E1). Note that since [H, S z tot ] 0, the derivative of the exponent −β∂H/∂h = −βS z tot can not be taken down straightforwardly from the exponential, i.e. (E5) where on the second line we have made the variable substitution τ ≡ λβ, and in the last second line we perform a cyclic permutation within the trace. For the sake of simplicity (and without loss of generality), we assume M = S z i 0 β = 0 in the discussion of fractional liquid, and thus arrive at the Kubo formula Eq. (3) in the main text.
If the total spin S z tot commutes with the Hamiltonian H, e.g. in the Heisenberg model, the two factors e τH and e −τH cancel each other and the Kubo formula reduces to the conventional susceptibility expression in terms of equal-time spin correlators. Here in the Kitaev model, although S z tot is not a good quantum number, e τH and e −τH nevertheless approximately cancel in the fractional liquid regime (see discussions in Sec. IV), and Eq. (5) of the main text also holds, which in turn gives rise to a Curie behavior with a modified Curie constant C > 1/4. However, as temperature further decreases to the KSL regime, "dynamical" correlations decays exponentially as the gauge fluxes freeze at low temperature below the finite flux excitation gap, which can suppress the universal Curie-law susceptibility at intermediate temperature.

Appendix F: Uniform and Staggered Magnetic Susceptibilities in the Antiferromagnetic Kitaev Model
For the isotropic Kitaev model with the AF Kitaev interaction K x,y,z = 1, we also compute the uniform magnetic susceptibility χ by applying small fields that coupled to S z component, i.e., along the [0 0 1] direction. In sharp contrast to the FM Kitaev model results in Fig. 5 of the main text, the uniform susceptibilities χ in Fig. A5(a) do not exhibit a universal Curie behavior in the intermediate-T regime. The susceptibility curves there are non-diverging convex curves and can not be fitted by the Curie-Weiss formula. Our cylindrical results of AF Kitaev susceptibility in Fig. A5(a) are again in excellent agreement with previous QMC results [29]. The underlying reason for the absence of emergent Curie behaviors under uniform fields is that the nearest spin pairs in the AF Kitaev model are also correlated along d z bond, while in an anti-parallel manner. Therefore, the change of total spin that relates to uniform susceptibility, corresponds to flipping one of the two spins, which has to overcome a finite amount of AF binding energy along the d z bond and does not reflect the fluctuating fluxes. Matter of fact, as shown in Fig. A5(b), the staggered susceptibility χ S under a staggered magnetic field h S (which has alternating sign on A and B sublattices of the Kitaev honeycomb) exhibits a universal Curie behavior at intermediate temperature.
In this section, we provide additional thermodynamic results of the Kitaev-Heisenberg and Kitaev-Gamma models. In Fig. A6, we plot the inverse susceptibility data 1/χ vs. T , under the Heisenberg [ Fig. A6(a)] and Gamma [ Fig. A6(b)] perturbations. The emergent Curie-Weiss susceptibility scaling is evident at intermediate temperature in Fig. A6, complementing the χ vs. T fittings in Fig. 6(g,h) of the main text.
In Fig. A7, we show contour plots of the specific heat C V and gauge flux W P results of both models. From Fig. A7(a,b), we see that the C V curves are always double-peaked, which defines two temperature scales T H and T L that bound the fractional liquid regime. The corresponding flux expectations W P are plotted in Fig. A7(c,d), where the asymptotic KSL regimes with nonzero W P can be clearly identified.
For the Kitaev-Heisenberg model, the low-T scale T L determined from the lower peak position of C V exhibits a clear dip at J c 0.12 in Fig. A7(a), and for J J c W P quickly establishes itself at around T L as shown in Fig. A7(c). The scenario is altered for J J c , where different low-T scale T L emerges, which represents the temperature scale between the SAF phase and Kitaev fractional liquid.
On the contrary, for the Kitaev-Gamma model in Fig. A7(b), the low-T scale T L remains stable for various Γ values, suggesting that the low-T KSL regime is very robust against Gamma perturbations. This is also supported by the flux results in Fig. A7(d), where W P remains nonzero at low temperature at least up to Γ = 0.1.