Hunting Majorana Fermions in Kitaev Magnets

A Majorana fermion is a fermionic particle that is its own antiparticle. Since the theoretical discovery by Ettore Majorana in 1937, the exotic particle has long been searched in particle physics. In the last few decades, however, it has attracted renewed interest in condensed matter physics, where it can be realized as an elementary excitation (quasiparticle) in quantum states of matter, such as the fractional quantum Hall states and topological superconductors. In this review, we discuss another platform for Majorana fermions, the quantum spin liquid. The quantum spin liquid is a bizarre quantum phase of insulating magnets, firstly proposed by Philip Anderson in 1973, in which interacting magnetic moments remain disordered down to the lowest temperature under strong quantum fluctuations. They are characterized by topological entanglement and fractional excitations, whose possible application to topological quantum computation is recently discussed intensively. As a prime candidate for such exotic states, we here focus on the Kitaev magnets, a subgroup of the spin-orbit Mott insulators, which have been a subject of intense research initiated by the seminal works by Alexei Kitaev in 2006 and by George Jackeli and Giniyat Khaliullin in 2009. After a brief overview of the Kitaev model and the fractionalization of spins in the exact ground state, we review recent explosive development in this rapidly growing field, with a focus on numerical solutions of the Kitaev model at finite temperatures and the comparison with experiments. The key concept is thermal fractionalization — two types of fractional excitations manifest themselves at largely different temperatures. This leads to distinct thermodynamics and spin dynamics in a variety of experimentally measurable quantities. We discuss such peculiar behaviors as the signatures of fractional quasiparticles, in careful comparison with the available experimental data for the candidate materials of the Kitaev magnets. Our review gives an overview of the current status of the identification of Majorana fermions in the Kitaev magnets, which would serve as a basis for further experimental and theoretical studies toward the manipulation of the exotic particles for topological quantum computation.


Introduction
Majorana fermions are charge-neutral spin-1=2 particles that are their own antiparticles. They were theoretically discovered by Ettore Majorana in 1937 in a real solution for the Dirac equation. 1) The Majorana fermions are distinguished from the ordinary fermions in the complex solution, called the Dirac fermions. The Dirac fermions are not their own antiparticles, and can be described by the annihilation and creation operators, f and f y , respectively. Two Majorana operators are defined by using f and f y as The definitions immediately yield that their creation and annihilation are equivalent: and they satisfy the anticommutation relation where ij is the Kronecker delta (i; j ¼ 1; 2). Equation (1) indicates that the occupied and unoccupied states of the Dirac fermion can be described by a pair of Majorana fermions. This means that one Majorana fermion carries half degrees of freedom of one Dirac fermion. Since the intriguing proposal by Ettore Majorana, the physical example of the exotic particles has long been sought in particle physics. Within the standard model, all the fermionic particles are the Dirac fermions, except for the neutrino. Thus, the neutrino has long been studied as a prime candidate for the Majorana fermion, but its nature is not settled yet. [2][3][4] Another candidates have been discussed for superpartners in the supersymmetry model, but no evidence was established to date.
In the last few decades, the Majorana fermions have attracted renewed interest by their possible realization in condensed matter physics. 5) In this case, they appear not as elementary particles but as elementary excitations (quasiparticles) in quantum states of matter. In general, quantum many-body states under electron correlations can host emergent quasiparticles, which have distinct nature from the constituent electrons. In some cases, the elementary excitations are described by more than one types of quasiparticles, which looks like the electrons are fractionalized into several particles. This is called fractionalization. For instance, in the two-dimensional (2D) fractional quantum Hall states, the elementary charge −e is fractionalized into fractional charges, e.g., e=3, and as a result, the elementary excitations of the system are described by emergent quasiparticles called anyons that do not obey either Dirac-Fermi or Bose-Einstein statistics.
In the context of the fractionalization, the emergence of Majorana fermions has been discussed for several different quantum states, such as the edge modes in the ¼ 5=2 fractional quantum Hall state, [6][7][8][9][10] the zero modes in p-wave superconductors, 8,11) and the bound states in topological superconductors. [12][13][14][15] Since these Majorana fermions originate from the fractionalization of fundamental particles, i.e., electrons, they acquire topological entanglement and intrinsically nonlocal nature. Owing to the unusual properties, the emergent Majorana fermions have drawn a great attention for the possible application to topological quantum computation. 16,17) In this review, we focus on another realization of Majorana fermions in insulating magnets called Mott insulators. In these systems, electrons are spatially localized due to strong electron correlations, and hence, the charge degree of freedom is inactive. Instead, what can be fractionalized here is the spin degree of freedom. Such a possibility of spin fractionalization has been discussed to take place in the quantum spin liquid (QSL), which is a quantum disordered state in the Mott insulators, firstly proposed by Philip Anderson in 1973. 18) In the QSL, any conventional symmetry breaking is precluded by strong quantum fluctuations, and the localized spins remain disordered but quantum entangled. Several types of QSLs have been predicted depending on the symmetry of the system, and they host their own fractional quasiparticles. 19,20) For instance, in the so-called Z 2 QSLs, the spin excitations are supposed to be fractionalized into two types of elementary excitations, spinons and visons; the spinons are charge-neutral spin-1=2 particlelike excitations, while the visons are topological excitations defined by their stringlike traces. 21,22) Most of such arguments, however, lack rigorous grounds, as there are less well-defined QSLs in more than one dimension. Thus far, tremendous efforts have been made for geometrically-frustrated antiferromagnets in two and three dimensions, but there are few examples where the ground state is rigorously shown to be a QSL. [23][24][25][26] A main difficulty lies in the lack of suitable theoretical methods: Any approximate theories may miss the essential aspects of the quantum entanglement in QSLs, and numerical methods require extremely high precisions to select out the true ground state from other quasi-degenerate states under strong frustration. Thus, it has remained a big challenge to identify fractional spin excitations in QSLs.
The situation has been changed dramatically over the past decade through two breakthroughs. One is the proposal of the exactly-solvable model in the seminal paper by Alexei Kitaev in 2006, 27) which is now called the Kitaev model. The model is a spin-1=2 model defined on a 2D honeycomb structure with bond-dependent interactions. The ground state is exactly obtained to be a QSL, in which the spin excitations are fractionalized into two types of quasiparticle excitations: itinerant spinon-like excitations, which are described by the Majorana fermions, and localized ones that constitute visonlike excitations. The other breakthrough was brought by Jackeli and Khaliullin in 2009. 28,29) They pointed out that the Kitaev model can be materialized in a class of the Mott insulating magnets with the strong spin-orbit coupling. Stimulated by their argument, several materials have been nominated as the candidates for the Kitaev QSL, such as iridium oxides A 2 IrO 3 (A = Li and Na) and a ruthenium trichloride α-RuCl 3 . These two breakthroughs have driven intense research for the Kitaev QSL from both theoretical and experimental viewpoints.
In the present article, we give an overview of the recent progress in this rapidly growing field. Several review articles are already available for the Kitaev QSL and its candidates. [30][31][32][33][34][35][36] Here we particularly focus on the finite-temperature (T) aspects of the fractional excitations, which are relevant to identify them in the candidate materials. Since the exact solution of the Kitaev model is limited to the ground state, the authors and their collaborators have developed several numerical techniques to study the finite-T properties, [37][38][39][40][41][42] and calculated the experimental observables, such as the specific heat and entropy, static spin-spin correlations, 38) magnetic susceptibility, inelastic neutron scattering spectra, spin-lattice relaxation rate in the nuclear magnetic resonance (NMR), [39][40][41] Raman scattering spectra, 43) and longitudinal and transverse components of the thermal conductivity. 44) Through the detailed comparison of the theoretical results with experimental data, signatures of the fractional excitations have been accumulated for the Kitaev candidate materials. We will discuss in detail such comparisons in this review.
The structure of this article is as follows. In Sect. 2, we introduce the Kitaev model and the fractional excitations derived from the exact solution for the QSL ground state. After introducing the Hamiltonian in Sect. 2.1, we briefly discuss the origin of the peculiar bond-dependent interaction in the Kitaev model in Sect. 2.2. In Sect. 2.3, we describe a Majorana representation of the spin operators, which is different from the original one introduced by Kitaev but useful for numerical techniques developed for finite-T calculations. After an overview of the exact QSL ground state and the fractional excitations in Sects. 2.4 and 2.5, respectively, we discuss the effects of finite T, an external magnetic field, and other exchange interactions in Sects. 2.6, 2.7, and 2.8, respectively. These additional effects on the Kitaev QSL are schematically summarized in the potential phase diagrams in Sect. 2.9.
In Sect. 3, we discuss one of the distinct aspects in the thermodynamics of the Kitaev model, which we call thermal fractionalization. In Sect. 3.1, as the prototypical behaviors, two successive crossovers are discussed for the Kitaev model on the 2D honeycomb structure. Then, a peculiar phase transition with time-reversal symmetry breaking is overviewed for a 2D triangle-honeycomb structure in Sect. 3.2. In Sect. 3.3, we showcase several unconventional phase transitions found for three-dimensional (3D) extensions of the Kitaev model, which can be regarded as gas-liquid-solid transitions in terms of the spin degree of freedom. We also briefly discuss spontaneous breaking of time-reversal symmetry in the 3D cases. These crossovers and phase transitions are summarized in Sect. 3.4. In Sect. 4, we introduce several candidate materials for the Kitaev QSL. We discuss the fundamental aspects of quasi-2D iridium oxides in Sect. 4.1, a ruthenium trichloride in Sect. 4.2, and 3D iridium oxides in Sect. 4.3.
In Sect. 5, we compare theoretical results for the Kitaev model with experimental data for the candidate materials, focusing on the quasi-2D materials. We discuss the two successive crossovers in the specific heat and entropy in Sect. 5.1, the saturation of static spin correlations measured from optical probe in Sect. 5.2, and peculiar T dependence of the magnetic susceptibility in Sect. 5.3. Then, we turn to the signatures of the fractional excitations in the spin dynamics: the dynamical spin structure factor measured in inelastic neutron scattering in Sect. 5.4 and the NMR relaxation rate in Sect. 5.5. From the comparison, we discuss the dichotomy between static and dynamical spin correlations as a signature of the thermal fractionalization. More direct signatures of fermionic excitations are discussed for the thermal conductivity in Sect. 5.6 and the Raman scattering in Sect. 5.7; in the latter, the unusual fermionic nature is clearly identified in a wide-T range. Finally, in Sect. 5 When the angular momentum l p ¼ 1 is coupled with the spin angular momentum s ¼ 1=2 by the SOC, the t 2g manifold is further split into the low-energy j eff ¼ 3=2 quartet and the high-energy j eff ¼ 1=2 doublet. Thus, the low-spin d 5 state ends up with the one-hole state in the j eff ¼ 1=2 doublet, as shown in the right panel of Fig. 2(a). The j eff ¼ 1=2 doublet comprises a time-reversal Kramers pair, which is described by

Majorana representation
In the seminal paper, Kitaev showed that the ground state of the model in Eq. (4) is exactly obtained by introducing a Majorana representation of the spin operators. 27) In the exact solution, each spin-1=2 operator is represented by four Majorana fermion operators. Later, another Majorana representation was introduced, which gives the same exact solution. [57][58][59] In this case, the spin-1=2 operator is represented by two Majorana fermions. In this article, we briefly review the latter Majorana representation, as it is used in the numerical simulations in the later sections. The advantage of the latter is in the size of the Hilbert space. The former Kitaev's representation doubles the Hilbert space and requires a projection to the original subspace to obtain physical results. It is not straightforward to deal with the projection in the numerical methods. 60) On the other hand, such a projection is not necessary in the latter representation, as the size of the Hilbert space is retained.
In the Majorana representation, we first apply the Jordan-Wigner transformation to the model in Eq. (4), by regarding the system as a one-dimensional chain composed of the x and y bonds; see Fig. 3. In the Jordan-Wigner transformation, the spin operators are rewritten by spinless fermion operators as ð1 À 2n i 0 Þa y i ; where a y i and a i are the creation and annihilation operators for the spinless fermions, respectively, and n i ¼ a y i a i is the number operator; a y i and a i satisfy the anticommutation relations as fa y i ; a j g ¼ ij ; fa y i ; a y j g ¼ 0; fa i ; a j g ¼ 0; where ij is the Kronecker delta. Then, by considering that the honeycomb structure is bipartite, the Hamiltonian in Eq. (4) is transformed into X hr 0 ;w;r;bi x ða r 0 ;w À a y r 0 ;w Þða r;b þ a y r;b Þ À J y 4 X hr;b;r 0 ;wi y ða r;b þ a y r;b Þða r 0 ;w À a y r 0 ;w Þ À J z 4 X r ð2n r;b À 1Þð2n r;w À 1Þ; where the subscripts b and w label the two sublattices in the rth unit cell with one z bond (see Figs. 1 and 3); the sums of r; b and r 0 ; w in the first and second terms are taken for all NN pairs on the x and y bonds, colored by blue and green in Fig. 3, respectively. Note that the so-called boundary terms (Color online) Schematic picture of the one-dimensional chains consisting of the x and y bonds, which are shown by the thick blue and green lines, respectively. The honeycomb structure in Fig. 1 is deformed into a brick-wall structure. The dotted square represents the unit cell including the r th z bond; the two sites are denoted as r; b and r; w [see Eqs. (10) and (13)].
In the Jordan-Wigner transformation in Eq. (8), the sites are numbered from the bottom left as partly shown in the figure.
appear in the Jordan-Wigner transformation for the systems under periodic boundary conditions. The boundary terms are nonlocal and hard to treat in the numerical simulations. One way to avoid this is to consider the systems under open boundary conditions. Another is just to neglect the boundary terms; their contributions are expected to be negligible in the thermodynamic limit. Next, we replace the spinless fermion operators by Majorana fermion operators. This is done by where γ and " are the Majorana fermion operators. These are the same as Eq. (1). By using Eqs. (11) and (12), Eq. (10) is rewritten into where r in the last term is defined on the z bond as The bond variable r in Eq. (14) satisfies the following relations: This means that each r is a constant of motion and takes AE1. Thus, f r g are Z 2 conserved quantities. It is worth noting that they are related with another conserved quantities called the Z 2 fluxes denoted by W p , which were introduced in the paper by Kitaev. 27) W p is defined for each hexagonal plaquette on the honeycomb structure as where the product is taken for the six sites on the plaquette p in the clockwise manner [see Fig. 4(a)]; " is the index for the bond connected to the site j which is not included in the sides of p, and j is the μth component of the Pauli matrices (S j ¼ ħ 2 j , where ħ is the reduced Planck constant and taken to be unity hereafter). By using the algebra of the Pauli matrices and the equations above, W p is also expressed as where the product is taken for the two z bonds belonging to the hexagonal plaquette p [see Fig. 4(b)].

Quantum spin liquid ground state
The Majorana representation of the Hamiltonian in Eq. (13) shows that the original spin model in Eq. (4) is mapped to the system with itinerant Majorana fermions f j g coupled with the Z 2 conserved variables f r g or the Z 2 fluxes fW p g via Eq. (18). The situation is schematically shown in Fig. 5. The Hamiltonian is in a bilinear form of f j g, namely, there is no quantum interactions between the Majorana fermions f i g; they interact only with the Z 2 variables f r g. This means that the Hamiltonian can be written in a block diagonalized form classified by different configurations of f r g or fW p g as follows. The total Hamiltonian matrix with the dimension 2 N is decomposed into a direct sum of the sectors specified by fW p g configurations. The number of fW p g configurations is 2 N=2 . The block Hamiltonian in each sector has thus the dimension 2 N =2 N=2 ¼ 2 N=2 , and it is represented by a N Â N bilinear form of Majorana operators with hopping matrix elements including fW p g as c-numbers. This decomposition enables one to find the ground state, in principle, by comparing the energies in all the sectors, as the energy in each sector is easily obtained for the noninteracting fermion problem.
For this problem, a mathematical proof, called Lieb's theorem, offers the exact solution for the lowest energy state. 61) This theorem tells the flux configuration which gives the lowest energy state in the systems with mirror symmetry with respect to the plane not including the lattice sites. In the present model on the honeycomb structure, we can apply this  theorem to the cases when at least two of three J are equal. The exact ground state for these symmetric cases is given in the sector with all W p ¼ þ1, which is called the flux-free state. On the other hand, Lieb's theorem does not apply to the cases with generic J . Nevertheless, by comparing the energies for different configurations of fW p g, the flux-free state is deduced to be the ground state in the entire parameter space of J . 27) The flux-free state is a QSL. This was explicitly shown by calculating the spin correlations. 62) The spin correlations have nonzero values only for the μ components on the NN μ bonds as well as the same sites, namely, hS i S j i ≠ 0 only for ¼ and i; j 2 hi; ji : ð19Þ All other spin correlations vanish. This is concluded from the fact that a spin operator S i flips two neighboring W p sandwiching the μ bond including the site i; only the combinations of S i and S j satisfying the condition in Eq. (19) conserve the flux-free configuration of W p (see Fig. 6). Thus, the spin correlations are extremely shortranged in the flux-free state. This means that the flux-free ground state does not break any symmetry of the system, and hence, it is a rare realization of the exact QSL in more than one dimension.
Note that Eq. (19) holds for arbitrary flux configurations. This suggests that spin correlations other than those in Eq. (19) are always zero even at finite T where fluxes with W p ¼ À1 are thermally excited. This is indeed confirmed by numerical studies introduced in Sect. 2.6 and Appendix.

Fractional excitations
For the flux-free ground state, there are two types of excitations. One is the excitations in terms of the itinerant Majorana fermions f j g, and the other is for the Z 2 fluxes fW p g. These are quasiparticle excitations arising from the fractionalization of the spin degree of freedom.
The former excited states are constructed by exciting complex fermions ff y k g, which are comprised as linear combinations of Majorana fermions f j g with complex amplitudes (see Appendix A.1). They are noninteracting fermions traversing on the honeycomb structure with the NN hopping. The dispersion relation is given by 27) Here, are the primitive translation vectors (see Fig. 1), whose lengths are taken to be unity. The dispersion relation in Eq. (20) is depicted in Fig. 7 for several sets of the parameters J x , J y , and J z . In the isotropic case with J x ¼ J y ¼ J z , EðkÞ becomes gapless at the point nodes located at the corners of the Brillouin zone (K and KA points), as shown in the upper panel of Fig. 7(d). Near the gapless nodal points, EðkÞ has a linear dispersion, similar to the Dirac nodes in the dispersion of π electrons in graphene, as shown in the inset of Fig. 7(d). This leads to the ω-linear dependence of the density of states (DOS) in the low-energy limit, as shown in the lower panel of Fig. 7(d). The gapless nature is retained for small anisotropy in J x , J y , and J z , despite a shift of the nodal points; see Figs. 7(c) and 7(e). The two nodal points approach each other while increasing the anisotropy by increasing jJ z j, and finally merge at some point, as exemplified for jJ z j ¼ jJ x j þ jJ y j in Fig. 7(b). With a further increase of the anisotropy, EðkÞ is gapped in the entire Brillouin zone, as exemplified in Fig. 7(a). 6. (Color online) Configurations fW p g for the states (a) S z i jÉi and (b) S y j jÉi, where jÉi represents the flux-free state. S i flips two W p on both sides of the μ bond connected to the site i. As the states with different fW p g are orthogonal to each other, hÉjS z i S y j jÉi ¼ 0. The magnitude of the excitation gap Á is plotted in the entire parameter space in Fig. 8(a). Á is zero in the center triangle defined by the conditions jJ x j jJ y j þ jJ z j, jJ y j jJ z j þ jJ x j, and jJ z j jJ x j þ jJ y j (dashed lines in the figure). Meanwhile, Á becomes nonzero in the other three outer triangles and increases as increasing the anisotropy in the Kitaev coupling; the contours are parallel to the gaplessgapped boundaries. The J z dependence of the gap is shown in Fig. 8(b) along the center vertical line in Fig. 8 indicating that Á increases linearly with J z in the gapped phase for J z > 1:5. Thus, there are two different phases with respect to the fermionic excitations associated with the itinerant Majorana fermions: the gapless phase including the isotropic point and the gapped one including the anisotropic limits.
On the other hand, the other types of excitations are generated by flipping W p from the flux-free ground state. 27) It turns out that they are always gapped and dispersionless reflecting the localized nature of W p . The lowest-energy excited state is given by a pair flip of neighboring two W p . The excitation gap Á f is plotted on the J x -J y -J z phase diagram in Fig. 8(c). The gap is nonzero in the entire parameter space, except for the anisotropic limits at the three corners of the phase diagram; it remains small in the gapped phases in Fig. 8(a) but becomes large rapidly in the gapless phase. As shown in Fig. 8(d), along the center vertical line in Fig. 8(c), Á f becomes maximum at the isotropic point with Thus, the two different types of the fractional excitations have distinct excitation spectra. The fermionic excitations from the itinerant Majorana fermions are dispersive and become both gapless and gapped depending on the anisotropy in the exchange coupling constants. Meanwhile, the Z 2 flux excitations are always gapped with a flat dispersion. The energy scales are also largely different for these two excitations; the bandwidth for the former is roughly set by half of the sum of three J , while the excitation gap for the latter is much smaller by more than one order of magnitude. This large difference in the energy scales affects the thermodynamics and the spin dynamics in a peculiar fashion, as discussed in the later sections.

Effect of finite temperature
The exact solution and related arguments above are limited to zero temperature (T ¼ 0). At finite T, the Z 2 flux excitations are generated by thermal fluctuations, and the exact solution is no longer available. As discussed in Sect. 2.4, however, the model in Eq. (13) is defined by noninteracting fermions coupled with thermally-fluctuating Z 2 variables f r g. As f r g are regarded as classical variables taking AE1, the situation is similar to the Falicov-Kimball model 63) and the double-exchange model with Ising spins. 64,65) This enables us to study the finite-T properties by developing numerical techniques similar to those used for such fermion models. The authors and their collaborators have developed the quantum Monte Carlo (QMC) method free from the negative sign problem 37,38,42) and the cluster dynamical mean-field theory (CDMFT). 39,40) These Majorana-based techniques are efficient to compute thermodynamic quantities, but they cannot be applied to the quantities not commuting with f r g, e.g., dynamical spin correlations. To overcome this difficulty, the authors and their collaborators have also developed the continuous-time QMC (CTQMC) method based on the Majorana representation. [39][40][41] The details of each method are presented in Appendix.
As will be described in detail in Sect. 3, an interesting finding at finite T is that the two distinct fractional excitations manifest themselves clearly in the thermodynamic behavior of the system. Specifically, in the 2D honeycomb case, the two largely different energy scales lead to two crossovers at largely different temperatures. One appears at T ¼ T H in the order of the characteristic energy scale of the itinerant Majorana fermions [more precisely, the center of mass (COM) of the DOS for the fermion band; see Sect. 3], and the other takes place at T ¼ T L in the order of the excitation gap in terms of the localized Z 2 fluxes. These two characteristic temperatures show up in many observables, not only thermodynamic quantities, but also the spin dynamics, as discussed in the later sections.

Effect of a magnetic field
Let us return to the flux-free ground state and discuss the effect of an external magnetic field at T ¼ 0. The Zeeman coupling to the magnetic field, Àh Á P i S i , spoils the exact solvability, because it makes the flux operators W p in Eq. (17) and r in Eq. (14) no longer conserved. (Note that the sign of the g factor is opposite to that for electron spins, as discussed in Sect. 2.2.) Nonetheless, Kitaev suggested an interesting possibility by using the perturbation theory with respect to the field strength. 27) In the perturbation theory, the lowest-order relevant term is in the third order of h as where h ¼ ðh x ; h y ; h z Þ and the Kitaev couplings are set to be isotropic, where the sites p 1 -p 6 and the bonds b 1 and b 2 are defined for the plaquette p as shown in Fig. 9(b). 44) Equation (23) shows that the weak magnetic field induces the complex secondneighbor hopping of the itinerant Majorana fermions coupled with the Z 2 bond variables f r g. This modulates the dispersion relation from Eq. (20) to 27) Thus, while the fermionic excitation in the isotropic case with J x ¼ J y ¼ J z has the gapless nodal points at the K and KA points [see Fig. 7(d)], the magnetic field opens a gap proportional toh / h x h y h z as Á ¼ 3 4 ffiffiffi 3 ph (see Fig. 10).
On the other hand, the flux gap Á f is almost independent ofh. These behaviors are plotted in Fig. 11. Interestingly, the model in the presence of the secondneighbor hopping in Eq.     equivalence shows that the gapped fermion band in the magnetic field is topologically nontrivial. The gapped state is a Majorana Chern insulator with the Chern number C ¼ AE1 (the sign is set by that of h x h y h z 27) ). Thus, similar to other topologically-nontrivial insulators with nonzero Chern numbers, the gapped topological state in the weak magnetic field is predicted to possess gapless chiral edge modes. 27) In contrast to the integer quantum Hall states, such chiral edge currents cannot be detected by electromagnetic measurements, as the Majorana fermions do not carry any electric charges; however, they could be observed by heat measurements (see Fig. 12). There are two distinct features in this thermal Hall effect by the Majorana fermions. One is that the thermal Hall conductivity divided by T is predicted to be quantized at half of that for the integer quantum Hall state. 27) This is because each Majorana fermion carries half degrees of freedom of an electron, as mentioned in Sect. 1. The other feature is that the half-quantized thermal Hall effect can be induced by the magnetic field in any direction, even in-plane directions. This is because the chiral Majorana edge currents are induced by the Zeeman effect enhanced on the spins near the edges (see Sect. III in Supplemental Material in Ref. 44), in contrast to the electric edge currents from skipping orbits by the Lorentz force. This interesting phenomenon specific to the Majorana fermions will be discussed in Sect. 5.8.
Beyond the perturbation theory, any rigorous argument is not available thus far. Nonetheless, many numerical studies have been performed to clarify the effect of the magnetic field at T ¼ 0. One of the earliest studies was done by the density matrix renormalization group for the Kitaev-Heisenberg model (see Sect. 2.8). 67) The Kitaev coupling was assumed to be isotropic and FM (J x ¼ J y ¼ J z ¼ J > 0), and the magnetic field was applied along the [111] direction with the strength h. The results indicate that the topologicallynontrivial QSL state predicted by the perturbation theory survives up to the critical field h c ' 0:018J, and turns into a topologically-trivial forced FM state above h c . This has been confirmed, e.g., by the exact diagonalization and other density matrix renormalization group calculations. [68][69][70][71][72] Recently, considerable attention has been drawn to the case with AFM Kitaev couplings. While the perturbation theory above is common to the FM and AFM cases, different aspects appear between the two cases when going beyond the perturbation. The most intriguing aspect is the possibility of another topological QSL in the intermediate-field region. [69][70][71][73][74][75][76] It was argued that the AFM Kitaev model undergoes successive phase transitions from the low-field QSL connected to the topological QSL in the perturbed region to another topological QSL, and to the forced FM state, while increasing the field. Although candidate materials with the AFM Kitaev couplings have not been identified thus far, this interesting possibility has attracted much interest. Note that recently there are several theoretical proposals for material realization of the AFM Kitaev couplings, for instance, by using f electrons 54) and polar asymmetry perpendicular to the honeycomb plane. 77) Finite-T calculations under a magnetic field are more difficult. For instance, we cannot apply the sign-free Majorana-based QMC method, since it assumes the conservation of fW p g and f r g. Nonetheless, one can study finite-T properties of the Hamiltonian with the effective magnetic field in Eqs. (22) and (23) derived from the perturbation, by using the sign-free Majorana-based QMC method. Such applications will be discussed in Sect. 5.8. In addition, a CTQMC technique has recently been developed and applied to the region where the negative sign problem is not severe, as discussed in Sects. 5.1, 5.4, and 5.5. 78)

Effect of other exchange interactions
As briefly mentioned in Sect. 2.2, in reality, there exist other types of the exchange couplings. A generic Hamiltonian proposed for realistic compounds is given by where ij denotes the type of ij bond, and the 3 Â 3 matrixĴ ij is parametrized, e.g., for the z bond aŝ Here, J Heis is the coupling constant for the isotropic Heisenberg interaction, and Γ and À 0 are for the symmetric  In (a), under the magnetic field h perpendicular to the sample plane, an electric field E (thermal gradient ÀrT ) causes unbalance in the edge electric (thermal) currents J high (J high Q ) and J low (J low Q ), which leads to the quantized (thermal) Hall effect. In contrast, in (b), the Majorana fermions do not carry electric charge, and hence, edge electric currents do not appear under the electric field; however, edge thermal currents can appear under a thermal gradient. In this case, the magnetic field is not necessarily perpendicular to the plane; any direction, even in the plane, leads to the thermal Hall effect. In this situation, each excited flux is associated with a Majorana zero mode as schematically shown in (b), which behaves as a nonabelian anyon (see Sect. 6).
off-diagonal interactions.Ĵ x andĴ y are obtained by C 3 rotations.
In the early stage of the research of the Kitaev QSL, the case with À ¼ À 0 ¼ 0 has been intensively studied. The model is called the Kitaev-Heisenberg model. Figure 13 displays the ground-state phase diagram obtained by the exact diagonalization. 79,80) In this case, the model exhibits at least four magnetically-ordered phases in addition to two regions of the Kitaev QSLs: FM, zigzag, Néel, and stripy phases. An important finding in this phase diagram is that the Kitaev QSL is found in narrow but finite parameter windows with nonzero J Heis in both FM and AFM Kitaev cases. The result suggests that the Kitaev QSL is not a singular property limited to the pure Kitaev model but survives against additional exchange couplings. This has encouraged material exploration for the Kitaev QSL.
Through such experimental exploration as well as computational studies of the spin Hamiltonians on the basis of first-principles calculations, it has been realized that beside the Heisenberg interaction, the symmetric off-diagonal interaction Γ plays a role. Indeed, Γ can be larger than J Heis from the perturbative argument. 81) Thus, the model including Γ has also been studied, 81,82) for which the ground-state phase diagram becomes richer. The effect of À 0 was also studied. 82) From the materials perspectives, the crucial question is how these other exchange interactions affect the QSL behavior in the exact solution for the Kitaev model. Unfortunately, in most of the candidate materials found thus far, the lowest-T state shows a magnetic order, such as the zigzag type and an incommensurate noncollinear type (see Sect. 4). An exception was recently found for H 3 LiIr 2 O 6 . 83) The absence of long-range ordering in this material was discussed on the basis of the stacking manner of the honeycomb layers, 84) the role of the hydrogens, 85,86) and the relevance of disorder in the exchange interactions, 87) but it remains unclear how to reconcile the sharp NMR lines observed down to the lowest T 83) to these scenarios. Figure 14 summarizes the arguments in the previous sections into the schematic phase diagrams. The phase diagrams are displayed for both cases with the FM and AFM Kitaev couplings, in the parameter space of temperature T, external magnetic field h, and other non-Kitaev interactions; the origin corresponds to the QSL ground state found in the exact solution for the Kitaev model.

Schematic phase diagram
Let us first discuss the FM case shown in Fig. 14(a), which is believed to be relevant to most of the existing candidate materials. As discussed in Sect. 2.8, the Kitaev QSL state survives in a finite region in the ground state against the non-Kitaev interactions. Above the threshold, the ground state exhibits some magnetic ordering whose spin structure depends on the detailed forms of the non-Kitaev interactions. The magnetic order is expected to survive at finite T due to the spin anisotropy as well as the three dimensionality, and the critical temperature T c will rise as the non-Kitaev interactions increase.
On the other hand, while raising T in the QSL region below the threshold, the system undergoes two crossovers as briefly mentioned in Sect. 2.6 (the details will be discussed in Sect. 3). Considering the realistic value of J $ 200{300 K for A 2 IrO 3 [88][89][90][91] and J $ 100{200 K for α-RuCl 3 , 68,91-94) the high-T crossover takes place at T H $ 80{110 K and $40{ 80 K, respectively. The temperature scales are significantly higher than T c for these compounds, T c $ 15 K and ∼7 K, respectively. On the other hand, T L $ 1{2 K and $0:5{1 K are lower than T N . Thus, we expect that the candidate   If this is the case, there is a considerable T window between T H and T c , where one can expect unconventional behavior arising from the fractionalization; this will be discussed in detail in Sect. 5.1. When applying the external magnetic field, as discussed in Sect. 2.7, the QSL survives up to a nonzero field strength, but it is taken over by the forced FM state in the larger field region. An interesting question is whether the QSL behavior can be captured in the candidate materials after the magnetic order is suppressed by the magnetic field. We depict Fig. 14(a) so that there is a narrow but nonzero window for such field-induced QSL. This intriguing possibility has attracted upsurge interest in α-RuCl 3 , as will be described in Sects. 5.4, 5.5, and 5.8. Figure 14(b) represents the corresponding phase diagram for the AFM Kitaev case. The overall structure is similar to the FM case in Fig. 14(a), but there is a qualitative difference in the behavior in the magnetic field. As described in Sect. 2.7, in the AFM Kitaev case, the system appears to exhibit two successive phase transitions including the intermediate QSL phase. [69][70][71][73][74][75][76] Note that the scale of the magnetic field is almost ten times larger compared to the FM case (this is also indicated by the large difference in the magnitude of the magnetic susceptibility in Sect. 5.3). Although no realistic compounds with the AFM Kitaev coupling are at hand thus far, the peculiar phase diagram is worth investigating and will stimulate further material exploration.

Thermal Fractionalization
In this section, we discuss a distinguished thermodynamic property of the Kitaev model, which we call thermal fractionalization. 38) As discussed in Sect. 2.5, the exact QSL ground state hosts two types of quasiparticles, itinerant Majorana fermions and localized Z 2 fluxes, which have largely separated energy scales. The two energy scales show up in the thermodynamic behavior as two characteristic temperatures. The higher characteristic temperature T H is related with the itinerant Majorana fermions, which is roughly set by the COM of the fermion DOS (see Fig. 7). At T ' T H , the system exhibits a crossover irrespective of the spatial dimensions as well as the details of the model. Meanwhile, the lower one T L is related with the localized Z 2 fluxes, which is roughly set by the Z 2 flux gap [see Fig. 8 In contrast to the universal crossover at T H , the behavior at T ' T L depends on the nature of the localized Z 2 flux excitations in each system; it can be either a crossover or a phase transition. Thus, the Kitaev model, in general, exhibits three distinct states: a conventional paramagnetic (PM) state for T ≳ T H , an unconventional PM state for T L ≲ T ≲ T H , and the (asymptotic) QSL state for T ≲ T L . We call the intermediate T region the fractional PM state, where the thermal fractionalization makes the system distinct from the conventional PM state.
We discuss these intriguing behaviors by the thermal fractionalization in this section. They have been unveiled by the recently-developed numerical methods based on the Majorana representation of the Kitaev model at zero field. In Sect. 3.1, we present the results for the 2D Kitaev model on the honeycomb structure, which provides a canonical example of two successive crossovers at T H and T L . We also discuss a variant of the Kitaev model in two dimensions in Sect. 3.2, which exhibits a phase transition to a chiral spin liquid (CSL), instead of the low-T crossover at T L . In Sect. 3.3, we present the results for the Kitaev models defined on several 3D tricoordinate lattices, in which various types of the phase transitions take place between three states of matter in terms of the spin degree of freedom. Finally, in Sect. 3.4, we summarize the phase diagrams for the crossovers and phase transitions found in the 2D and 3D Kitaev models.

Successive crossovers in the 2D honeycomb case 3.1.1 Crossovers caused by thermal fractionalization
Let us begin with the original Kitaev model defined on the honeycomb structure. Figure 15 shows the T dependences of the internal energy E, specific heat C v , and entropy S per site for the isotropic Kitaev coupling J x ¼ J y ¼ J z ¼ J 38) (the results are common to the FM and AFM Kitaev couplings). The calculations were performed by using the QMC simulations based on the Majorana representation for the clusters with N ¼ 2L 2 spins (see Appendix A.1). As shown in Fig. 15(a) and its inset, the internal energy E decreases rapidly at two temperatures, T H ' 0:375J and T L ' 0:012J, while the decrease at T L is much smaller than that at T H . Correspondingly, the specific heat C v exhibits two peaks as shown in Fig. 15(b), both of which show no significant system-size dependence, indicating that these are crossovers. Interestingly, as plotted in Fig. 15(c), the entropy S is released successively by half ln 2 at each crossover. This peculiar behavior is considered to originate from the thermal fractionalization in which the original spin degree of freedom carrying the entropy of ln 2 is fractionalized into the two types of quasiparticles each carrying the entropy of half ln 2. This is confirmed by the decomposition of C v and S into the contributions from the itinerant Majorana fermions and the localized Z 2 fluxes, as shown in Figs The role of the two fractional quasiparticles in the two crossovers is shown in more explicit way by calculating the quantities associated with each quasiparticle. Figure 16(a) plots the measure of the kinetic energy of the itinerant Majorana fermions, K x ¼ Àih i j i x , where the thermal average hÁ Á Ái x is calculated on the x bond. Note that this quantity is related with the internal energy as E ¼ À 3 2 K x in the isotropic case. Also, it is equivalent to the spin correlation on the x bonds, 4hS x i S x j i x . The result indicates that the measure of the Majorana kinetic energy increases rapidly around T ¼ T H , and does not change largely in the lower-T region. This suggests that the Fermi degeneracy of the complex fermions composed of the itinerant Majorana fermions sets in at T ' T H . On the other hand, Fig. 16 displays the thermal average of the Z 2 flux, hW p i. While it becomes nonzero from high T around T H , it grows rapidly around T ¼ T L and approaches hW p i ¼ 1 (the value in the flux-free ground state) below T L .
These results clearly show that the crossover at T H is caused by the itinerant Majorana fermions, and that at T L is by the localized Z 2 fluxes. The former corresponds to the Fermi degeneracy of the complex fermions composed of the itinerant Majorana fermions, and the latter to the asymptotic freezing of the Z 2 fluxes into the flux-free state. Thus, these two crossovers are manifestations of the thermal fractionalization in thermodynamics. While decreasing T, the fractionalization of the spin degree of freedom sets in around T H with the entropy release of half ln 2 by the Fermi degeneracy, and the system enters into an unconventional PM state, dubbed the fractional PM state, below T ' T H . In the fractional PM region, the Z 2 fluxes remain disordered as the states with flipped W p are thermally excited beyond the flux gap. By approaching T L with a further decrease of T, however, the thermal excitations of the Z 2 fluxes are suppressed, and the system crosses over into the asymptotic QSL state below T ' T L with the entropy release of the rest half ln 2 by the freezing of W p . The picture of the successive crossovers will be further discussed in Sect. 3.4.

Crossovers temperature scales
What determines the values of the two crossover temperatures T H and T L ? From the above arguments, it is naturally expected that T H is set by the Fermi degeneracy temperature, which is roughly given by the COM of the fermion DOS, and that T L is set by half of the gap for the lowest excitation of the Z 2 fluxes (the flux gap is defined for a two-flux excitation); see Sect. 2.5. In the isotropic case with J x ¼ J y ¼ J z ¼ J, the COM of the fermion DOS is at '0:762J and the half of the flux gap is '0:0328J. Note that the COM of the fermion DOS is less sensitive to the flux sector, but here we use the value for the disordered flux configuration, which we call G , corresponding to the high-T limit, as the fluxes are almost disordered near T H as shown in Fig. 16(b). Considering that the specific heat peak in the two-level system with a gap of unity appears at T ¼ ' 0:417, we note that the numbers 0:762J Â ' 0:318J and 0:0328J Â ' 0:0136J are very close to T H ' 0:375J and T L ' 0:012J, respectively, which confirms the above expectation.
We can further examine these correspondences by varying the anisotropy of the Kitaev coupling. Figure 17 shows the contour plot of the entropy per site, S, normalized by ln 2 while changing J z with J x ¼ J y and J x þ J y þ J z ¼ 3. The two white regions with S=ln 2 ' 0:75 and '0:25 roughly corresponds to T H and T L , respectively. As shown in the figure, T H does not show a drastic change against J z , whereas T L does: T L has a peak around the isotropic point with J z ¼ J and rapidly decreases by increasing the anisotropy with both J z ! 0 and J z ! 3. For comparison, we plot the effective activation temperatures defined by the COM of the fermion DOS for the disorder flux configuration, G , and the Z 2 flux gap, 1 2 Á f , by the solid and dashed curves in Fig. 17. The former does not change so much for J z similar to T H (see

Majorana metal
Let us discuss the excitation spectrum of the itinerant Majorana fermions while changing T. Figure 19 shows the T dependence of the fermion DOS. While the overall structure of the DOS below T L is similar to that in the flux-free ground state as shown in Fig. 19(b), the DOS rapidly changes its form above T L . In particular, in the low-energy region, the energy-linear behavior at the lower band edge is quickly smeared out and the DOS at zero energy becomes nonzero, as shown in Fig. 19(a). 38,95) This is due to the thermal excitations of the Z 2 fluxes, which disturb the Dirac-like linear dispersion in the flux-free ground state. The nonzero DOS at the band bottom indicates that the fractional PM state above T L is regarded as a "Majorana metal", in analogy with the 2D conventional metal that has nonzero DOS at the band edges. Needless to say, the present system is an insulator with localized magnetic moments, and hence, the particles traversing the system are not electrons but the Majorana fermions. This is why we call the unconventional state the Majorana metal. An interesting consequence of this Majorana metallic state is observed in the specific heat C v . As shown in Fig. 20, C v shows T-linear dependence in the T window between T L and T H , reflecting the "metallic" nature of the system. 38) The itinerant quasiparticles also contribute heat conductions, as will be discussed in Sects. 5.6 and 5.8.

Phase transitions to 2D chiral spin liquids
Let us turn to a variant of the Kitaev model in two dimensions, which is defined on a modified lattice structure, The solid and dashed curves represent G and 1 2 Á f , respectively, where G and Á f are the COM of the fermion DOS for the disordered flux configuration corresponding to the high-T limit and the flux gap for the fluxfree ground state, respectively; ' 0:417 is the peak temperature of the specific heat in the two-level system with a gap of unity. The factor 1=2 in 1 2 Á f is introduced since Á f means the energy cost to excite two neighboring fluxes. The data of the entropy were taken from Ref. 38. The two curves are newly added.   called the triangle-honeycomb structure. The structure is obtained by replacing all the vertices of the honeycomb structure by triangles, as shown in Fig. 21. The Kitaev model can be extended straightforwardly to this tricoordinate lattice structure, 27) and one can define two different sets of the Kitaev coupling, ðJ x ; J y ; J z Þ and ðJ 0 x ; J 0 y ; J 0 z Þ, for the two types of NN bonds, intra-triangle and inter-triangle ones, respectively (see Fig. 21). 96) In the following, we consider the case The most important difference from the honeycomb model is that the lattice structure includes the elementary loops with odd number of sites. As pointed out in the seminal paper by Kitaev,27) the Kitaev model defined on the lattices with such odd cycles may break time-reversal symmetry spontaneously, as the flux operator defined on an odd-cycle plaquette describes a time-reversal pair. Indeed, Yao and Kivelson showed that the ground state of the triangle-honeycomb Kitaev model becomes a CSL with spontaneous breaking of time-reversal symmetry. 96 under open boundary conditions. The topological nature was also explained by the fact that the low-energy effective model in the limit of J 0 =J ! 0 has a similar form to the effective Hamiltonian for the honeycomb model in a magnetic field with the three-spin term in Eq. (22). 97) Thermodynamic properties of this model were studied by using the Majorana-based QMC simulations. 98) In contrast to the honeycomb case in Sect. 3.1, the model exhibits a finite-T phase transition instead of the crossover at T L . This is due to the spontaneous breaking of time-reversal symmetry by the freezing of the Z 2 fluxes; while the freezing does not break any symmetry in the honeycomb case, it breaks time-reversal symmetry in the triangle-honeycomb case because of the odd-cycle plaquettes on the triangles. Interestingly, the transition was found to be continuous in the topologicallynontrivial region for J 0 =J ≲ ffiffiffi 3 p (the estimated critical exponents are close to those of the 2D Ising universality class), but discontinuous in the topologically-trivial region for J 0 =J ≳ ffiffiffi 3 p . This suggests the existence of the tricritical point in between, while the precise location and the nature are not fully identified. The obtained phase diagram is presented in Fig. 22.
Then, what happens to the high-T crossover at T H found in the honeycomb case? It was shown that while the crossover takes place also in the triangle-honeycomb case, the amount of entropy released in the crossover can be different from the honeycomb case depending on the parameter J 0 =J. 98) In the topologically-trivial region for J 0 =J ≳ ffiffiffi 3 p , the entropy release is the same as in the honeycomb case, half ln 2. But in this case, the system exhibits another crossover at a lower T, where the entropy of 1 6 ln 2 corresponding to fluxes on the dodecagons is released. Finally, the rest of entropy 1 3 ln 2 corresponding to fluxes on the triangles is released at the phase transition to the CSL. The typical T dependences of the specific heat C v and the entropy S per site are shown in Fig. 23. On the other hand, in the topologically-nontrivial region for J 0 =J ≲ ffiffiffi 3 p , the entropy release at the high-T crossover is 1 3 ln 2; the remaining entropy 2 3 ln 2 corresponds to fourfold degeneracy in each triangle in the isolated triangle  limit (J 0 =J ! 0). In this case, the system exhibits two additional crossovers, at each of which the entropy of 1 6 ln 2 is released; see the typical behavior in Fig. 23. As mentioned before, in the isolated triangle limit, the system is effectively described by the honeycomb Kitaev model in a weak magnetic field, and hence, these two crossovers correspond to T H and T L in the honeycomb Kitaev model. Note that the lowest-T crossover appear to merge into the phase transition for J 0 =J ≳ 0:1 (see Fig. 22).
Thus, the complicated behaviors are found in the high-T crossovers depending on two types of the Kitaev coupling, J and J 0 . Nonetheless, the important point is that the highest-T crossover occurs at the temperature almost independent of J 0 =J (T Ã in Fig. 22). This originates from the Fermi degeneracy of the complex fermions composed of the itinerant Majorana fermions, similar to the honeycomb case in Sect. 3.1. Hence, the comparison between the honeycomb and triangle-honeycomb cases implies that the high-T crossover arising from the itinerant Majorana fermions is commonly seen in the variants of the Kitaev model, while the low-T one from the localized Z 2 fluxes may appear differently depending on the nature of the Z 2 flux in each system. We will further examine this conjecture in several examples in three dimensions in Sect. 3.3.

Phase transitions in three dimensions
In this section, we discuss the thermodynamic behaviors in some variants of the Kitaev model in three dimensions. In Sect. 3.3.1, we present the results for the 3D Kitaev model on the so-called hyperhoneycomb structure. In this model, the low-T crossover at T L in the 2D honeycomb case in Sect. 3.1 is replaced by a phase transition as in the triangle-honeycomb case in Sect. 3.2; however, the low-T phase is not a CSL but the Kitaev QSL in this case. We discuss the origin of the phase transition to the QSL on the basis of the distinct nature of the Z 2 flux excitations in three dimensions. In Sect. 3.3.2, we present the results for the model which exhibits a phase transition to a conventional magnetically-ordered phase in addition to that to the QSL. Finally, we discuss the phase transitions to 3D CSLs on a lattice structure with odd-cycle plaquettes, dubbed the hypernonagon lattice in Sect. 3.3.3.

Phase transition by loop proliferation: Gas-liquid transition
Mandal and Surendran discussed a variant of the Kitaev model on a 3D lattice structure, 99) which was later called the hyperhoneycomb structure shown in Fig. 24. The lattice is in a series of extensions of the honeycomb structure to three dimensions, 100) and surprisingly, it is realized in a candidate material β-Li 2 IrO 3 101) (see Sect. 4.3 for the details). Mandal and Surendran showed that the hyperhoneycomb Kitaev model retains the exact solvability and the ground state offers a 3D exact QSL. They also argued the peculiar nature of the Z 2 flux excitations, which we will discuss below.
Finite-T behavior of this 3D model was studied by the Majorana-based QMC simulations. 37) T dependences of the specific heat and entropy are presented in Fig. 25 for the isotropic case with J x ¼ J y ¼ J z ¼ J. Although the overall behaviors are similar to those for the 2D honeycomb case in Figs. 15(b) and 15(c), clear differences appear at low T; while the low-T peak in the specific heat does not depend on the system size in the 2D case, the height becomes higher and the width gets narrower in the present 3D case as increasing the system size [see also the inset of Fig. 25(a)]. This signals a phase transition instead of the crossover. A similar phase transition was also found by larger-scale simulations for the effective model in the anisotropic limit of the Kitaev coupling, called the Kitaev toric code. 102) In the present case, however, in contrast to the transition to the CSL in Sect. 3.2, the freezing of the Z 2 fluxes does not break timereversal symmetry, as the lattice structure does not include odd cycles. Then, what happens in this finite-T phase transition?
The phase transition is caused by a change of topological nature in the excitations of the Z 2 fluxes. 37,102) In the 3D case, the localized Z 2 fluxes cannot be flipped independently because of the local constraint arising from the lattice geometry. 99 What happens in the 3D case is as follows. While raising T from the flux-free QSL ground state, the localized Z 2 fluxes are thermally excited in the form of closed loops. At low T, the loop lengths are short compared to the system size. With a further increase of T, however, excitation loops with their lengths comparable to the system size are proliferated at some T because of the entropic gain, which leads to the topological transition. The critical temperature T c is set by the loop tension arising from the excitation energy proportional to the loop length. 103) Thus, the finite-T phase transition in this 3D Kitaev model is caused by the loop proliferation. The picture of this topological transition will be further discussed in Sect. 3.4.
The phase transition takes place between the high-T PM state and the low-T QSL state. The former is regarded as "gas" in terms of the spin degree of freedom, while the latter is regarded as "liquid", both of which preserve the symmetry of the system. Therefore, the phase transition is regarded as a "gas-liquid" transition in the spin degree of freedom. In contrast to the conventional gas-liquid transition, which is discontinuous in general, the numerical results in Fig. 25 do not find any discontinuity. The analysis of the effective model in the anisotropic limit concludes that the phase transition is continuous and belongs to the inverted 3D Ising universality class; the confined loops are favored in the low-T (high-T) phase in the 3D toric code (Ising model). Note that the closed loops in the 3D Ising model are composed of interacting spins, which appear in the high-T expansion and contribute to the partition function. The order parameter of this peculiar transition is not described by any local quantities but it can be identified by a global quantity called the Wilson loop, which is given by the product of all W p on the plane defined by a given loop. 37) Note that the Wilson loop measures the parity of the total number of the excited W p lines penetrating the plane. Thus, this phase transition caused by the loop proliferation evades from the conventional Landau-Ginzburg-Wilson theory for the continuous phase transitions.
A similar phase transition was found also for another 3D Kitaev model defined on the so-called hyperoctagon lattice. 42) The origin of the phase transition is common. This suggests that the loop proliferation works as a common mechanism for the gas-liquid phase transition in 3D Kitaev models. The comparative study between the hyperhoneycomb and hyperoctagon cases confirmed the correlation between T c and the loop tension. 42)

Three states of matter: Gas-liquid-solid transition
Stimulated by the finding of the gas-liquid phase transition in the spin degree of freedom, the phase transitions for three states of matter, gas, liquid, and solid, were investigated for the Kitaev toric code with additional ferromagnetic Ising interaction. 104) In this model, while increasing the Ising interaction, the QSL ground state is taken over by a FM ordered state, which is regarded as "solid". Hence, one can expect the phase transitions between the three states of matter. Figure 26(a) shows the phase diagram obtained by  extensive QMC simulations (in this case, not the Majoranabased QMC but the continuous-time world-line QMC in the original spin representation). 104) The result indicates that the gas-liquid transition described in Sect. 3.3.1 survives against the FM Ising interaction with a slight decrease of the critical temperature, but at some point, it changes into a phase transition between the high-T PM state and the low-T FM state, which is a gas-solid transition. The first-order transition line between the QSL and FM phases extends from T ¼ 0 to the tricritical points on the gas-liquid and gas-solid transition lines (the former is not identified within the numerical precision). In the PM state near the bifurcation of the phase boundaries, an interesting proximity effect was found in the flux loop excitations. 104) Similar study was conducted also for the 2D case. 104) The result is shown in Fig. 26(b). In contrast to the 3D case in Fig. 26(a), the QSL phase is limited to zero T, while there is a crossover at finite T. The crossover T decreases as the Ising interaction increases, and finally goes to zero at the quantum critical point. For larger Ising interactions, the FM state evolves with continuous growth of T c . The phase transition at T c is continuous and belongs to the 2D Ising universality class. Thus, the phase transitions for three states of matter in the spin degree of freedom look qualitatively different between the 3D and 2D cases, owing to the distinct nature of the Z 2 flux excitations.
The above study of three states of matter has been limited to the toric code corresponding to the anisotropic limit of the Kitaev coupling. The issue in a more realistic parameter region remains for future study, which is potentially relevant to understanding of the properties of 3D candidate materials for the Kitaev model (see Sect. 4.3). The 2D case is also worth investigating; 105) indeed, in a weakly anisotropic case, an interesting liquid-liquid phase transition between the Kitaev QSL and a spin-nematic quantum paramagnet was found before entering the FM ordered state. 106)

Phase transitions to 3D chiral spin liquids
In Sect. 3.2, we discussed finite-T phase transitions to 2D CSLs with spontaneous breaking of time-reversal symmetry. Similar transitions in three dimensions were studied for the 3D Kitaev model defined on the lattice structure with odd cycles, dubbed the hypernonagon structure. 107,108) In the 3D case, there is an interesting possibility of successive phase transitions, since the 3D Kitaev models can exhibit a topological transition by the loop proliferation discussed in Sect. 3.3.1, in addition to the spontaneous time-reversal symmetry breaking. Such a possibility was studied for two anisotropic limits of the Kitaev coupling in the hypernonagon Kitaev model. 107) The numerical results indicate that the system exhibits a single discontinuous phase transition with simultaneous occurrence of the loop proliferation and timereversal symmetry breaking. Interestingly, however, the low-T CSL state is not a flux-free state but shows a nonuniform spatial order of the Z 2 fluxes. The study was extended to other parameter regions apart from the anisotropic limits, and at least five distinct phases with different nonuniform flux orders were discovered. 108) Most of the studies of CSLs thus far have been limited to two dimensions since the pioneering work by Kalmeyer and Laughlin. 109

Phase diagram
As a brief summary of Sect. 3, we show schematic phase diagrams at finite T for the Kitaev models in both two and three dimensions. Figure 27(a) displays the 2D honeycomb case, 38) which will be common to other 2D cases without odd cycles in the lattice structure. In this case, the system exhibits two crossovers at T ¼ T H and T L . The former temperature scale is set by the COM of the itinerant fermion DOS, and the latter by the excitation gap for the localized Z 2 fluxes. The finite-T state is separated into three by these two crossovers: the conventional PM for T ≳ T H , the fractional PM for T L ≲ T ≲ T H , and the asymptotic QSL for T ≲ T L . The schematic picture of each region is shown in the lower panels of Fig. 27(a).
Meanwhile, Fig. 27(b) displays the 3D counterpart, inferred from the results for the 3D hyperhoneycomb 37,102) and hyperoctagon cases. 42) In this case, while the high-T crossover at T H remains in a similar manner, the low-T one is replaced by the phase transition of gas-liquid type caused by the loop proliferation. The difference arises from the distinct nature of the localized Z 2 flux excitations. The schematic picture in terms of the excited loops is shown in the lower panels of Fig. 27(b). The lower panels show the schematic pictures of the three states. The magenta spheres, the gray hexagons, and the arrows represent the itinerant Majorana fermions, the flipped localized Z 2 fluxes, and the spins, respectively (see Fig. 5). In (b), the systems undergo a "gas-liquid" phase transition at T c from the low-T QSL to the fractional PM and a crossover at T H to the conventional PM. In the schematic picture in the lower panels, the cyan and purple lines represent short and extended loops composed of the flipped localized Z 2 fluxes, respectively. In (c), the phase transition at T c occurs between the low-T CSL and the fractional PM states. Figure 27(c) presents the schematic phase diagram in the presence of odd cycles in the lattice structure. In this case, the system can show a finite-T phase transition to a CSL with spontaneous breaking of time-reversal symmetry. The transition is a single phase transition from the fractional PM state to the CSL for the 2D triangle-honeycomb 98) and the 3D hypernonagon cases. 107,108) It can be split into multiple transitions in the latter as discussed in Sect. 3.3.3, but such behavior has not been found thus far.

Material Candidates
In this section, we briefly overview candidates for materialization of the Kitaev QSL. We here focus on some of 4d-and 5d-electron compounds. Readers who are interested in more details including other candidates are referred to other review articles. 31,32,35)

Quasi-2D iridates
As discussed in Sect. 2.2, Jackeli and Khaliullin pointed out two requisites for materialization of the Kitaev coupling. They nominated A 2 BO 3 -type layered compounds as a good candidate. Following this proposal, Chaloupka and his coworkers have pointed out that this is indeed the case for the quasi-2D honeycomb iridium oxides, Na 2 IrO 3 and α-Li 2 IrO 3 . 79) These two compounds have a common quasi-2D lattice structure with the honeycomb layers composed of edge-sharing IrO 6 octahedra, [110][111][112] as shown in Fig. 2(c); the crystal symmetry belongs to space group C2=m. (We put the prefix α only for Li 2 IrO 3 since it has polymorphs as introduced in Sect. 4.3.) In these compounds, the formal valence of the Ir ions is 4+, and hence, the outermost 5d shell is partially occupied by five electrons. As described in Sect. 2.2, this leads to the low-spin 5d 5 state under the cubic crystalline electric field, and furthermore, comprises the j eff ¼ 1=2 pseudospin with the influence of the strong spinorbit coupling [see Figs. 2(a) and 2(b)]. The importance of both spin-orbit coupling and Coulomb interaction and the formation of the j eff ¼ 1=2 state have been confirmed by spectroscopic measurements. 113,114) The pseudospins are expected to interact with each other via the Kitaev coupling through the perturbation processes in the edge-sharing geometry [see Fig. 2(c)]. The predominant Kitaev coupling was experimentally confirmed for Na 2 IrO 3 by using diffuse X-ray scattering 115) and torque magnetometry. 116) It was also supported by theoretical estimates based on first-principles calculations. [88][89][90][91] Despite the presence of the predominant Kitaev coupling, these candidates do not show QSL behavior in the low-T limit; instead, they undergo a phase transition to a magnetically-ordered phase at low T. Na 2 IrO 3 exhibits a zigzagtype AFM ordering at the critical temperature T N ' 15 K, 110,117,118) while α-Li 2 IrO 3 exhibits an incommensurate spiral ordering at almost the same T. 111,112,119) The magnetic orders are considered to be induced by non-Kitaev couplings in the honeycomb layer as well as interlayer couplings, which are weaker than the Kitaev coupling. Thus, it is widely believed that the compounds are proximate to the Kitaev QSL, whereas the low-T properties are hindered by the parasitic magnetic orders.
There have been several efforts to realize the Kitaev QSL by suppressing the non-Kitaev interactions. Theoretically, it was proposed that a thin film 90) and a heterostructure 91) might be helpful for this purpose. Also, experimentally, the chemical substitutions of A-site ions locating between the honeycomb layers were attempted, and A 0 3 LiIr 2 O 6 with A 0 = Ag, Cu, and H were synthesized. 83,120,121) These compounds have a different stacking manner from Na 2 IrO 3 and α-Li 2 IrO 3 . Among them, H 3 LiIr 2 O 6 is intriguing since it does not show any magnetic ordering down to the lowest T, 83) while disorder effects have been argued, as discussed in the end of Sect. 2.8. In addition, Cu 2 IrO 3 was recently nominated as a candidate, but in this case also the effect of chemical disorder was pointed out. [122][123][124][125]

α-RuCl 3
Another candidate is a ruthenium trichloride α-RuCl 3 , which was firstly pointed out in Ref. 126. This compound has a similar quasi-2D layered honeycomb structure with edgesharing RuCl 6 octahedra, but the crystal symmetry is controversial among P3 1 12, C2=m, and " R3 depending on the samples. 94,[127][128][129][130][131][132] This might be related with the fact that the honeycomb layers are weakly coupled with each other via the van der Waals interaction. 133) The formal valence of the Ru ions is 3+, and hence, the 4d 5 electron configuration offers a playground for the Kitaev coupling similar to the iridium oxides in the previous section. The formation of the j eff ¼ 1=2 state was confirmed, e.g., by the spectroscopic measurements with the help of first-principles calculations. 68,91,126,134,135) Unfortunately, this compound also exhibits a magnetic order of zigzag type at low T. 131,132,136) The critical temperature is, however, scattered between T N ' 6:5 K and '14 K depending on the samples. It is believed that the samples with stacking faults show rather high T N ; the lowest T N ¼ 6:5 K was reported for a single crystal with " R3 symmetry. 94) One of the advantages in α-RuCl 3 is the feasibility of inelastic neutron scattering, which is a powerful tool to probe spin dynamics (note that Ir is a neutron absorber). Recently, several measurements have been done in a wide range of energy and wave vector. The results will be discussed in comparison with theoretical results for the Kitaev model in Sect. 5.4.
Another advantage is that the zigzag magnetic order in α-RuCl 3 can be suppressed by an external magnetic field of ∼8 T applied within the ab plane. 130,131) This opens an interesting possibility to realize QSL behavior in the fieldinduced PM region. We will discuss the recent development on this issue in Sects. 5.4, 5.5, and 5.8.
Last but not least, α-RuCl 3 has a unique aspect owing to the fact that this compound is a van der Waals material: The weak interlayer coupling allows to fabricate the samples in a thin film form. [137][138][139][140] More recently, interesting electronic properties were observed for heterostructures between a thin film of α-RuCl 3 and graphene. [141][142][143][144] Such fabrication of thin films and heterostructures will stimulate further studies on interesting physics arising from the potential fractional excitations in this Kitaev candidate magnet.

3D iridates
Finally, we introduce two polymorphs of Li 2 IrO 3 : β-Li 2 IrO 3 and γ-Li 2 IrO 3 . These two compounds have 3D networks of the edge-sharing IrO 6 octahedra, instead of the quasi-2D layered one in α-Li 2 IrO 3 . β-Li 2 IrO 3 has the so-  101) and γ-Li 2 IrO 3 has the stripy-honeycomb structure with space group Cccm [ Fig. 28(b)]. 100) Both structures belong to a series of the harmonic honeycomb structures. 100) In both cases, the local coordination is common to α-Li 2 IrO 3 , and the Ir ions comprise tricoordinate lattices, for which the Kitaev model can be extended in a straightforward manner. Thus, these polymorphs have attracted attention as candidates for the 3D Kitaev QSL discussed in Sect. 3.3.1. 145,146) However, they show spiral magnetic ordering at rather high temperature T N $ 40 K. 100,101,147,148) Interestingly, the magnetic orders can be suppressed by applying relatively small magnetic fields 149,150) as well as external pressure. 151,152)

Comparative Study between Theory and Experiment
In this section, we discuss the signatures of thermal fractionalization in the Kitaev QSL through the comparison between theory and experiment. On the theoretical side, we concentrate on the Kitaev model in Eq. (4) defined on the honeycomb structure, neglecting other additional interactions discussed in Sect. 2.8, as it allows to obtain reliable results by well-controlled numerical techniques. All the following results are for the isotropic Kitaev coupling J x ¼ J y ¼ J z ¼ J. Meanwhile, on the experimental side, we present the data for three candidates: the honeycomb iridium oxides, Na 2 IrO 3 and α-Li 2 IrO 3 , and the ruthenium trichloride α-RuCl 3 .

Specific heat and entropy
Let us first begin with the comparison for the specific heat and entropy. Figure 29(a) displays the experimental data for a candidate material for the Kitaev model, Na 2 IrO 3 . 153) The specific heat exhibits a broad peak around 110 K, in addition to a sharp anomaly at T N ' 15 K associated with the magnetic ordering. While decreasing T, the entropy is released corresponding to the high-T broad peak, and shows an interesting T dependence with inflection points; the decrease becomes slow around 60 K, where the entropy is roughly half R ln 2 (R is the gas constant). With a further decrease of T, the entropy is continuously released, and finally, decreases rapidly at the magnetic phase transition at T N ' 15 K. Qualitatively similar behaviors were observed for the related compound α-Li 2 IrO 3 153) and another candidate α- At first glance, these experimental data look similar to the theoretical results for the Kitaev model presented in Sect. 3.1.1, except for the sharp anomaly at the magnetic transition temperature. Then, it is natural to ask whether the similarities provide experimental evidence for the thermal fractionalization arising from the Kitaev QSL. The answer is that although they look consistent with theory, it is difficult to admit them as strong evidence. On one hand, the broad peak in the specific heat at high T is in fact commonly seen in frustrated magnets; the suppression of magnetic ordering by the frustration leaves development of short-range spin

012002-19
© 2020 The Physical Society of Japan correlations, which gives rise to the entropy release in the high-T region. This is also the case in the Kitaev model: As shown in Fig. 16(a), the crossover at T ¼ T H is related with the growth of NN spin correlations. Hence, the broad peak of the specific heat alone cannot be evidence of the thermal fractionalization. On the other hand, the approximately half R ln 2 entropy at the shoulderlike feature also looks consistent with the theoretical result, but this is again not conclusive, considering that in general it is not easy to precisely estimate the lattice contributions in experiments. Also, theoretically, it is difficult to predict how non-Kitaev interactions, which are inevitably present in real materials, affect the behavior of the entropy at low T. 154,155) Then, what could be evidence in these thermodynamic quantities? A specific feature to the Kitaev QSL is the low-T crossover at T ¼ T L by the freezing of the localized Z 2 fluxes. Unfortunately, in the candidate materials shown above, T L ' 0:012J is considered to be around 1 K, which is lower than the critical temperatures. Thus, the interesting behavior associated with the Z 2 fluxes, if any, is hindered by the parasitic magnetic ordering caused by non-Kitaev interactions. A potential route to unveil the crossover behavior is to suppress the magnetic order by applying an external magnetic field, as discussed in Sect. 2.7. Such an experiment was indeed performed for α-RuCl 3 , and a peak was observed in the region where the magnetic order is suppressed by the magnetic field. [156][157][158] Meanwhile, the specific heat in the magnetic field was recently calculated for the Kitaev model by using a newly-developed CTQMC method; 78) a similar peak was obtained in the high-field region, while the data at low T and low field are lacked because of the negative sign problem. Further detailed comparison is necessary to identify the signature of the Z 2 fluxes.

Spin correlation
The equal-time spin correlation was indirectly obtained for α-RuCl 3 by an optical measurement. 159) In this experiment, several peaks were identified in the optical conductivity above the Mott gap ∼1 eV, as shown in Fig. 30(a). Among them, the lowest-energy peak just above the Mott gap, denoted as α in Fig. 30(a), shows considerable T dependence. As this excitation reflects virtual motions of electrons beyond the Mott gap, the T dependence is considered to contain the information on the development of spin correlations originating from the virtual exchange processes. The T dependence of the spectral weight of the peak α is shown in Fig. 30(b). The data show that, while decreasing T, the spectral weight grows down to ∼40 K, whereas it almost saturates in the lower-T region, even below the critical temperature T N . This behavior resembles the T dependence of the NN spin correlations in the Kitaev model plotted in Fig. 16(a), where T H $ 35 K by assuming J $ 8 meV. Nevertheless, the change of the weight in Fig. 30(b) is rather small (∼10%), which might be due to contributions from other excitations in the optical spectrum and the T independent contributions in the virtual exchange processes. It is desired to make further quantitative comparison and also to perform more direct measurement of the spin correlations.

Magnetic susceptibility
The magnetic susceptibility for the Kitaev model was calculated by Majorana-based numerical techniques, [39][40][41] by using the formula where ¼ 1=ðk B TÞ is the inverse temperature (we set the Boltzmann constant k B ¼ 1), and hS z i ðÞS z j i is the dynamical spin correlation in the (2 þ 1)-dimensional space, where S i ðÞ ¼ e H S i e ÀH (τ is the imaginary-time). Figure 31 shows the results obtained by the combined technique between the Majorana-based QMC and CTQMC methods 41) (see Appendix A.1 and A.3). Note that all the off-diagonal components with ≠ vanish in the Kitaev model, 62) and xx ¼ yy ¼ zz ¼ in the isotropic case.
As shown in Fig. 31, although the T dependence as well as the overall magnitude is different between the cases with FM and AFM Kitaev coupling, the two cases share the following features. (i) At sufficiently high T, χ obeys the Curie-Weiss law as in other magnets, but it starts to deviate below T $ J. The Curie-Weiss behavior is given by ¼ 1=ð4T À JÞ for the FM case and ¼ 1=ð4T þ JÞ for the AFM case. (ii) χ exhibits a peak in the fractional PM region between T L and  (iv) In the low-T limit, χ approaches a nonzero value. Similar asymptotic behavior is commonly seen in the magnetic systems in which the total spin is not conserved. In the present case, owing to the Dirac-like linear dispersion in the fermionic excitations, the asymptotic behavior is expected to be proportional to T 3 up to a constant, but it is hard to extract such behavior from the present numerical results. 160) For comparison, we showcase the experimental data for the candidate materials in Fig. 32. The data for Na 2 IrO 3 in Fig. 32(a) shows that the susceptibility obeys the Curie-Weiss law above ∼150 K, but starts to deviate at lower T. 110) Similar behavior was observed also for α-Li 2 IrO 3 111) [see also Fig. 32(b) 112) ]. In both cases, a peak appears at a slightly higher T than the critical temperature T N ' 15 K. Below the peak, the susceptibility turns to decrease and exhibits an inflection point around T N , and finally approaches a nonzero constant at the lowest T. These behaviors appear to be at least qualitatively similar to the theoretical results in Fig. 31, although one cannot compare the data below T N . Similar behaviors were observed for α-RuCl 3 130,136) [see Fig. 32(c)]. Meanwhile, a readily-seen discrepancy between theory and experiment is the magnetic anisotropy in the experimental data, as shown in Fig. 32. The theoretical results are isotropic for the isotropic case, and it is also difficult to explain the magnetic anisotropy by the anisotropy in the Kitaev coupling. 40) The importance of additional non-Kitaev interactions as well as the anisotropy of the g factor was pointed out for the magnetic anisotropy. 68,[161][162][163] It remains as a future issue to quantitatively explain the T dependence of the anisotropic susceptibility and to determine the magnitude and sign of the Kitaev coupling. We will comment on the sign of the Kitaev coupling in Sect. 5.4.
Can we say that the comparison between theory and experiment for the magnetic susceptibility provides evidence for the proximity to the Kitaev QSL? As in the case of the specific heat and entropy in Sect. 5.1, the similarity found in the T dependence is suggestive but not sufficient to draw

012002-21
© 2020 The Physical Society of Japan conclusions. This is because the deviation from the Curie-Weiss behavior and the broad peak structure at a lower T are commonly observed in a wide class of frustrated magnets as a consequence of the growth of short-range spin correlations under the frustration. A more decisive feature would be an experimental observation of the rapid decrease around T L with the inflection point originating from the freezing of the localized Z 2 fluxes. This is, however, hindered again by the magnetic ordering in the real compounds.

Inelastic neutron scattering
Inelastic neutron scattering is a powerful experimental tool to probe the spin dynamics. The scattering intensity is proportional to the dynamical spin structure factor which includes the information on the spin dynamics as a function of wave vector q and frequency ω. Figure 33(a) displays the theoretical results for the QSL ground state of the Kitaev model for both FM and AFM cases. 164,165) Here, the dynamical spin structure factor is calculated by where r ij is the vector connecting the sites i and j, and S i ðtÞ is the Heisenberg representation of S i . In the following,
There are several interesting features in the results shown in Fig. 33(a). (i) The q dependence is weak. This is due to the fact that the Kitaev model possesses extremely short-ranged spin correlations, as discussed in Sect. 2.4. (ii) The intensity vanishes below the rather strong response at low energy ! $ 0:4J z which corresponds to ! $ 0:1J in our definition (see the figure caption). This is due to the gap opening in the flux excitations. As discussed in Sect. 2.5, the spins are fractionalized into the Majorana fermions and the Z 2 fluxes, meaning that the elementary spin-flip excitation is given by a composite of the Majorana fermion excitation and the Z 2 flux excitation. Hence, the spin excitation spectrum in Sðq; !Þ reflects the gap opening in the fractional excitations of the Z 2 fluxes. At the same time, this suggests that the strong intensity above the gap predominantly originates from the Z 2 flux excitations. (iii) In addition to the low-energy response, the spectrum has a broad incoherent intensity in the highenergy region extending up to ! $ 6J z corresponding to ! $ 1:5J in our definition. This reflects mainly the itinerant Majorana fermion excitations, which has the bandwidth $1:5J as shown in Sect. 2.5.
Recently, the inelastic neutron scattering measurements have been intensively performed for α-RuCl 3 . In an early experiment for powder samples, an unusual incoherent intensity was observed in the energy range of ! ¼ 6{8 meV, in both below and above the critical temperature T N , as shown in Figs. 34(a) and 34(b), respectively. 93) This is clearly distinguished from the strong response at a lower energy only appearing below T N [indicated by the white arrow in Fig. 34(a)] which is regarded as the spin-wave excitations in the ordered phase. The incoherent response at high energy has a resemblance to that in the theoretical result at T ¼ 0 shown in Fig. 33(a). More interestingly, it remains visible up to ∼70 K, which is much higher than T N , as shown in Fig. 34(c). 93) Nevertheless, at this stage, there was no theory for the T dependence for comparison.
The T dependence of the dynamical spin structure factor Sðq; !Þ was computed by using the combined techniques between the Majorana-based CDMFT and CTQMC, 39,40) and the Majorana-based QMC and CTQMC methods 41) (see Appendix). The calculations were done by where S i;j ð!Þ is obtained by solving by using the maximum entropy method with the Legendre polynomial expansion (see Ref. 40 for the details). The results obtained by the Majorana-based QMC and CTQMC method are shown in Fig. 33(b). There are several interesting features. (i) At sufficiently high T in the conventional PM region above T H [lowest two panels in Fig. 33(b)], Sðq; !Þ has an almost q-independent broad peak centered around ! ¼ 0. (ii) While approaching T H with a decrease of T, however, an incoherent response gradually grows at high energy centered at ! $ J [middle-lower two panels in Fig. 33(b)]. This high-energy incoherent response persists down to the lowest T, gradually developing a weak q dependence. (iii) With a further decrease of T toward T L , a quasi-elastic response grows in the low-energy region [center and middle-upper panels in Fig. 33(b)]. (iv) Below T L , this quasi-elastic response is shifted to the ! > 0 region with opening of a small gap [upper two panels in Fig. 33(b)], and the entire spectrum smoothly converges to the T ¼ 0 results shown in Fig. 33(a).
The contrasting T dependence between the high-energy incoherent response and the low-energy quasi-elastic response reflects the distinct nature between the two types of fractional excitations arising from the thermal fractionalization. As discussed in Sect. 3.1, the crossovers at T ¼ T H and T L are caused by the itinerant Majorana fermions and the localized Z 2 fluxes, respectively. Therefore, the growth of the high-energy incoherent response in Sðq; !Þ below T ' T H

012002-23
© 2020 The Physical Society of Japan is considered to be dominated by the itinerant Majorana fermions, while that of the quasi-elastic response toward T ' T L as well as the gap opening below T L is by the localized Z 2 fluxes. This is consistent with the assignment discussed above for the T ¼ 0 results in Fig. 33(a).
After the theoretical studies for T > 0, inelastic neutron scattering experiments were performed for single crystals of α-RuCl 3 . 94,166) An example is shown in Fig. 35. As observed in the theoretical results in Fig. 33(b), an unconventional incoherent response appears in a wide range of energy up to ∼12 meV below ∼100 K. (The differences in the energy and T scales from Fig. 34 might be ascribed to the sample difference. 94,167) ) The detailed comparison with theory in Fig. 35 indicates that, in the wide-T range from the conventional PM region to just above T N , the overall q and ω dependences of the spectra can be accounted for by the theoretical results for the Kitaev model with isotropic FM coupling. Although the growth of the quasi-elastic response toward T L as well as the gap opening below T L predicted by theory was not observed in experiments because of the magnetic ordering at T N , the agreement strongly suggests that the candidate material α-RuCl 3 is in proximity to the Kitaev QSL. 168) Despite the overall good agreement, there remain some discrepancies between theory and experiment, especially at low T and low ω. A representative feature is a star shape in the q dependence of the scattering intensity at low ω above T N . 94,166) The numerical results for the Kitaev model show a round shape, unlike the star one. 94) The coexistence of such a low-energy feature and the high-energy incoherent response was discussed by considering the role of additional non-Kitaev interactions. 166,[170][171][172][173] Let us briefly comment on the sign of the Kitaev coupling. The comparison in Fig. 35 indicates that the FM Kitaev coupling well accounts for the weak q dependence in the experimental data. In the earlier studies, 93,166) however, the AFM Kitaev coupling was deduced from the comparison between experiment and theory for the weak q dependence of the high-energy continuum. The AFM Kitaev coupling was also suggested by theory based on first-principles calculations. 174) On the other hand, other theoretical studies based on quantum chemistry electronic-structure calculations 68) and first-principles calculations 91) suggest the FM Kitaev coupling. We note that, in a later experimental study by the same group, 175) the FM Kitaev coupling was deduced from the careful analyses of the spectral weights.
More recently, inelastic neutron scattering experiments have been done in a magnetic field. 175,176) The experimental data for a powder sample of α-RuCl 3 are shown in Fig. 36. The results show that the spin excitation spectrum in the region where the magnetic order is suppressed by the magnetic field [ Fig. 36(b)] is qualitatively similar to that above T N at zero field [ Fig. 36(a)]; the low-energy contribution from magnon excitations observed in the ordered phase [ Fig. 36(c)] is absent, and the high-energy incoherent response is commonly observed in Figs. 36(a) and 36(b). This suggests that an unconventional state potentially described by the Kitaev QSL is realized in the field-induced PM region [see the phase diagram in Fig. 36(d); see also Theoretically, it is hard to obtain reliable results in a magnetic field since the exact solvability is lost and the Majorana-based numerical techniques cannot be applied straightforwardly, as described in Sect. 2.7. However, the spin dynamics was recently obtained by a CTQMC method in a wide range of field and T. 78) It was shown that Sðq; !Þ preserves the unconventional features reflecting the fractional excitations in the wide-field range before entering the forced FM region in the high field and low T. This may explain the unconventional spectrum in the field-induced PM state discovered in the experiment above. Moreover, the theoretical result unveiled a crossover behavior from the fractional quasiparticle picture to the conventional magnon picture while increasing the magnetic field, which is one of the confinement-deconfinement phenomena. 78) Similar issue was studied by the exact diagonalization of a 24-site cluster for a model including non-Kitaev interactions. 177) While an experiment was performed recently, 176) further detailed comparison between theory and experiment is highly desired for these interesting issues.

Nuclear magnetic resonance
In addition to the inelastic neutron scattering, the NMR is an important probe of the spin dynamics. The NMR relaxation rate is a measure of the dynamical spin susceptibility through the formula 178) 1 where A q is the hyperfine coupling constant, ? ðq; ! 0 Þ is the dynamical susceptibility for the spin component perpendicular to the field direction, and ! 0 is the resonance frequency in the NMR measurement. Note that the dynamical susceptibility ðq; !Þ is related with the dynamical spin structure factor discussed in the previous section through the fluctuation-dissipation theorem, as In the NMR measurements, ! 0 is in general negligibly small compared to the typical energy scale of the system, J in the present case. Thus, by taking the limit of ! 0 ! 0 in Eq. (32) and using Eq. (33), one can obtain where S ? ðq; !Þ is the dynamical spin structure factor for the spin components perpendicular to the field. The NMR relaxation rate 1=T 1 was calculated for the Kitaev model by using the Majorana-based numerical techniques. [39][40][41] The calculations were done in the limit of zero field, which correspond to the nuclear quadrupole resonance (NQR) in experiments. Considering the fact that the Kitaev model has nonzero spin correlations only for the same site and between the NN sites (see Sect. 2.4), 1=T 1 for the magnetic field along the z direction is computed by Eq. (34) as 1=T z 1 ¼ a 0;x S xx j; j ð! ¼ 0Þ þ a 0;y S yy j; j ð! ¼ 0Þ þ a 1;x S xx NN ð! ¼ 0Þ þ a 1;y S yy NN ð! ¼ 0Þ; ð35Þ where S NN ð!Þ represents the NN component on the μ bond [see also Eqs. (30) and (31)]. In Eq. (35), the coefficients a 0;x , a 0;y , a 1;x , and a 1;y are determined by the hyperfine coupling constant A q depending on the details of the actual compounds. Figure 37 shows the results for (a) the onsite and (b) NNsite components separately, defined as respectively; here we omit the coefficients a 0;x , a 0;y , a 1;x , and a 1;y . In Eq. (37), the sign is + (−) for the FM (AFM) case [S NN ð! ¼ 0Þ changes sign but the absolute value is the same for both cases]. Note that 1=T x 1 ¼ 1=T y 1 ¼ 1=T z 1 ¼ 1=T 1 for the isotropic case. Comparison to experiments can be made for the superpositions with appropriate coefficients determined by A q . The results in Fig. 37   spin correlations are almost the same for the two components after the fractionalization sets in. (iii) While decreasing T below T H , both components grow in the fractional PM region and show a broad peak at T ' 0:04J. (iv) Both components are rapidly suppressed around T L . This is ascribed to the gap opening in the Z 2 flux excitations, as observed in Sðq; !Þ in Sect. 5.4. Indeed, the low-T behaviors are well fitted by the activation-type function proportional to expfÀaÁ f =ðk B TÞg, where Á f is the flux gap and a is a coefficient. 160) An interesting feature among these behaviors is the growth of 1=T 1 in the fractional PM region below T H . This means that the dynamical spin correlations are developed in this T region. On the other hand, as discussed in Sects. 3.1.1 and 5.2, the equal-time spin correlations are almost saturated below T H and do not show significant T dependence. These observations indicate that the Kitaev model exhibits distinct T dependences between the dynamical and static spin correlations below T H where the thermal fractionalization sets in. Such dichotomy is hardly seen in conventional magnets, except for critical behaviors in magnetic ordering. Thus, the strong enhancement with the broad peak in 1=T 1 under the saturated static spin correlations would be an indication of the thermal fractionalization in the Kitaev QSL.
Related to this enhancement, let us make a remark on the Korringa law. As introduced in Sect. 2.5, the system is described by noninteracting Majorana fermions coupled to the Z 2 fluxes. Indeed, the T-linear specific heat is observed in the fractional PM region, which we call the Majorana metal in Sect. 3.1.3. From this picture, one might expect the Korringa law, 1=ðT 1 T 2 Þ $ constant, which holds for free fermion systems, in the same T region. The numerical data, however, do not support this expectation. 40) This might be due to the fact that the spin-flip excitation is a composite of both itinerant Majorana fermions and localized Z 2 fluxes, as discussed in Sect. 5.4. NMR measurements have been done for α-RuCl 3 by several groups. [180][181][182][183] The representative data are shown in Fig. 38. In the low-field region for ≲9 T where the magnetic ordering takes place at low T, 1=T 1 grows gradually while decreasing T, and shows a sharp anomaly at the critical temperature T N , followed by a rapid decrease below T N . On the other hand, in the higher-field region where the magnetic order is suppressed, 1=T 1 grows gradually but turns to decrease after showing a broad peak, as shown in Fig. 38(a). While increasing the magnetic field, the peak height is gradually decreased and the peak temperature is shifted to higher T. The low-T decrease is well fitted by the activationtype function, as shown in Fig. 38(b), while 1=T 1 appears to approach a nonzero constant or show a slight increase at the lowest T measured in this experiment. We note, however, that the low-T behaviors of 1=T 1 are scattered among the data from different groups. For instance, in Ref. 181, the powerlaw T dependence was observed in some range of the magnetic field, from which the existence of gapless excitations was concluded. Meanwhile, from the measurement down to 0.4 K in Ref. 183, another exponential decrease was found at the lower-T region than measured in the previous studies, from which two gap structure was identified.
As mentioned above, the theoretical results in Fig. 37 are obtained in the zero-field limit, which correspond to NQR, and hence, the direct comparison with the experimental data is not straightforward. 184) Nevertheless, it is interesting to point out that the experimental data in the high-field PM region look similar to the theoretical results in the points (iii) and (iv) raised above, while the low-T asymptotic behaviors are controversial in experiments. This suggests the possibility that the fractional PM state is realized in the magnetic field. Indeed, in Ref. 182, the authors proposed an empirical function for the T dependence of 1=T 1 by analyzing the theoretical results at zero field, and showed that it fits well the experimental data in a wide range of the magnetic field, as presented in Fig. 38(c). Interestingly, the estimates of the gap by this fitting procedure appear to be consistent with the prediction from the perturbation theory in Sect. 2.7: The gap is proportional to h 3 up to a constant.
Theoretical analysis was recently extended to nonzerofield regions by a CTQMC method. 78) The results indicate that the overall behavior of 1=T 1 is retained in a wide range of T and field; in particular, the broad peak structure is preserved with a decrease of the peak height and a shift of the peak temperature to a higher T while increasing the magnetic field. These behaviors are apparently consistent with the experimental data shown in Fig. 38. The agreement suggests that the Kitaev model qualitatively explain the behavior of 1=T 1 in the field-induced PM region, and furthermore, that the fractional PM state appears to extend to a wide-field region in the candidate material α-RuCl 3 .

Thermal conductivity
In Sects. 5.4 and 5.5, we have discussed the signatures of the thermal fractionalization in spin dynamics. As mentioned above, however, a spin-flip excitation is a composite  One suitable probe is thermal transport. This is because in the Kitaev QSL heat is carried solely by the itinerant Majorana fermions, as the Z 2 fluxes are completely localized. The thermal response is measured as the thermal conductivity , which is defined by J Q ¼ r T, where J Q is the thermal current flowing in the α direction induced by the thermal gradient applied to the β direction, r T. Here, ; ¼ x; y, which correspond to the a and b directions of the Cartesian coordinate shown in Fig. 1.
In order to capture the itinerant nature of Majorana fermions, the longitudinal component of the thermal conductivity, ¼ , was measured for α-RuCl 3 . 185) The results are shown in Fig. 39(a) for several samples. Figure 39(b) shows the results after careful subtraction of the contributions from phonons. The data indicate that there are additional contributions in a wide-T range centered at ∼100 K.
Theoretical results were obtained for the Kitaev model almost at the same time by using the Majorana-based QMC method. 44) In the calculations, the thermal current J Q is defined by the time derivative of the energy polarization P E as where P E is introduced from the Hamiltonian by replacing the exchange constant J on the bond hiji to J R ij with R ij ¼ 1 2 ðr i þ r j Þ. Using the above definitions, the longitudinal thermal conductivity was computed by the Kubo formula given as where J Q ðtÞ is the Heisenberg representation of J Q , and V is the volume of the system. Note that the thermal current   operator J Q commutes with all the Z 2 bond variables r in the Hamiltonian in Eq. (13), and therefore, Eq. (39) can be calculated by using the sign-free Majorana-based QMC technique in Appendix A.1. Figure 40 shows the longitudinal thermal conductivity in the isotropic case of the honeycomb Kitaev model. Note that xx ¼ yy and the result is common to the FM and AFM cases. The result indicates that the thermal conductivity exhibits a broad peak around T H . This is a direct consequence of the thermal fractionalization; the itinerant Majorana fermions appear in the system when the thermal fractionalization sets in by approaching T H from high T, but their thermally-activated population decreases with a further decrease of T because of the Fermi degeneracy. The theoretical result resembles qualitatively the experimental data in Fig. 39.

Raman scattering
The comparison of the thermal conductivity in the previous section suggests the existence of heat carriers in the insulating material besides phonons. However, it is not straightforward to conclude that the carriers are Majorana fermions. In this section, we discuss another measurement, the Raman scattering, which could probe the Majorana fermions more directly.
The Raman scattering is a powerful tool to identify the magnetic excitations by using light. Theoretically, the intensity of the Raman scattering spectrum is calculated as 186) Here, R is the Loudon-Fleury operator 187) given by where in and out are the polarization vectors of the incoming and outgoing lights, and d is the vector connecting a NN μ bond for the sites i and j. Note that, in the isotropic case with J x ¼ J y ¼ J z ¼ J assumed here, there is no polarization dependence. 186) Figure 41 shows the Raman scattering intensity Ið!Þ calculated for the exact QSL ground state of the Kitaev model. 186) The spectrum includes a broad incoherent response in a wide-energy range up to about 3J [note that the energy scale in Fig. 41 is four-times larger than the present definition, as in Fig. 33(a)]. Equations (40) and (41) indicate that on the basis of the Loudon-Fleury approach 187) the Raman response originates solely from the itinerant Majorana fermions in the Kitaev QSL, as S i S j are written by the Majorana operators i and j and do not affect the Z 2 variable configurations f r g. 186) Hence, the broad response in Fig. 41 is a direct consequence of the fermionic excitations with the wide bandwidth shown in Sect. 2.5. Figure 42(a) displays the experimental result measured for α-RuCl 3 at 5 K. 92) In addition to the sharp peaks around 14 and 20 meV, which are presumably from phonon excitations through the spin-lattice coupling, the spectrum exhibits a broad incoherent response ranging up to ∼25 meV, as indicated by the blue shade in Fig. 42(a). This incoherent response is similar to that found in the theoretical result at T ¼ 0 in Fig. 41, suggesting the existence of the itinerant Majorana fermions. Figure 42(b) displays the magnetic contributions for several T, and Fig. 42(c) plots the T dependence of the intensity integrated between 2.5 and 12.5 meV. 92) In conventional magnets, T dependence of the intensity is usually well fitted by using the Bose-Einstein distribution function nð!Þ, since the excitations are given by magnons and phonons, both of which obey the Bose-Einstein statistics. The result plotted in Fig. 42(c) shows that this is not the case for α-RuCl 3 : There are additional contributions that cannot be fitted by using nð!Þ in the wide-T range. This peculiar T dependence could be evidence of the Majorana fermions, but there was no theoretical result at finite T at this stage.
Finite-T behaviors of the Raman scattering intensity for the Kitaev model were obtained by using the Majorana-based QMC method. 43) Note that this dynamical quantity can also be calculated by the sign-free QMC technique in Appendix A.1, since the Loudon-Fleury operator R commutes with all r as  the thermal current operator J Q in Sect. 5.6. Figure 43(a) displays the results for the T and ω dependence. While increasing T from the ground state, the incoherent nature of the spectrum is retained, but the weight distribution changes gradually; the low-ω weight increases continuously up to T ' T H and saturates above T H , while the weight around ! ¼ J shows a slight increase up to T ' 0:05J, which is slightly above T L , but turns to decreases at higher T. 43) Figure 43(b) presents the comparison between theory and experiment. In this comparison, the experimental spectrum is assumed to be a simple summation of the magnetic contribution from the Kitaev model and that from (unidentified) bosonic excitations. The result in Fig. 43(b) shows that the T dependence of the Raman intensity integrated in the middle-energy range from 5 to 12.5 meV is well reproduced by the theoretical result in a wide-T range by assuming J ¼ 10 meV. Furthermore, the theoretical calculations showed that the dominant contribution in this T range comes from pair creations and annihilations of the emergent fermions composed of the Majorana fermions. 43) This is indeed seen from the fact that the theoretical result is well reproduced by a simple function ð1 À f Þ 2 , where f is the Femi-Dirac distribution function, as indicated by the green dashed curve in Fig. 43(b). The good agreement between theory and experiment strongly suggests the existence of fermionic excitations in the experimental data in the wide PM region above T N , which are absent in conventional magnets.
After this surprising result, similar analyses were performed for another candidates, iridium oxides β-and γ-Li 2 IrO 3 188) (see Sect. 4.3). Despite the 3D honeycomblike structures in these compounds, the Raman scattering intensity exhibits similar T dependence, which is well fitted by ð1 À f Þ 2 . This suggests that the fermionic excitations are commonly present in the candidate materials for the Kitaev QSL. Also, we note that contributions of non-Kitaev interactions 189,190) and an external magnetic field 191) were recently discussed.

Thermal Hall conductivity
The unconventional contribution in the Raman intensity strongly suggests the existence of fermionic excitations, but it is still difficult to conclude that the excitations are nothing but the Majorana fermions, especially solely from the experimental data. To prove the existence of the Majorana fermions, one needs to explicitly identify the consequence from their peculiar nature, e.g., the equivalence between the particle and its anti-particle. In this section, we discuss one of such direct consequences discovered in the recent measurements of the thermal Hall transport.
As discussed in Sect. 2.7, Kitaev showed by using the perturbation theory that a weak magnetic field induces a gapped topologically-nontrivial state showing the halfquantized thermal Hall effect due to the chiral Majorana edge mode 27) (see Fig. 12). As the half quantization is a direct consequence of the fact that the Majorana fermions carry half degrees of freedom of the electrons, its measurement provides a smoking gun for the Majorana nature.
Prior to experiments, T dependence of xy was numerically calculated by using the Majorana-based QMC method for the effective model derived by the perturbation theory given by Eqs. (13) and (23). 44) In the calculations, a contribution from "the gravitational magnetization" was taken into account in addition to the Kobo formula similar to the longitudinal case given in Eq. (39). 192,193) Figure 44 shows the results. While decreasing T, xy =T increases gradually below T $ J, and approaches rapidly the half quantized value =12 below T L . The low-T asymptotic behavior is fitted by / expðÀÁ f =TÞ, where Á f is the flux gap. Interestingly, xy =T shows nonmonotonic T dependence in the intermediate-T region, originating from thermal fluctuations of the localized Z 2 fluxes which scatter the itinerant Majorana fermions. 44) The corresponding experiment was performed for α-RuCl 3 . 194) The results are shown in Fig. 45(a). The T dependence above T N is qualitatively similar to that in the theoretical results replotted in Fig. 45(b); xy =T becomes nonzero below ∼80 K and shows a broad peak above T N . While further decreasing T, however, the experimental data decrease and change the sign to negative below T N . In this experiment, the magnetic field was applied along the c axis, which cannot suppress the magnetic order in the field range measured, and hence, the half quantization, if any, is hindered by the magnetic ordering.
Recently, xy =T was measured in the magnetic field tilted from the c axis, which can suppress the magnetic ordering. 195) Note that α-RuCl 3 has strong easy-plane anisotropy as shown in Fig. 32(c). The typical experimental data for the T dependence are shown in Fig. 46(a). As shown in the inset, xy =T increases from zero below ∼60 K and once overshoots the half quantized value below ∼20 K. With a further decrease of T, xy =T turns to decrease, and below ∼5.5 K, it becomes almost T independent as shown in the main panel; the asymptotic constant value indeed coincides with the half quantized value within the experimental errors. The field dependences at low T are presented in Fig. 46(b). The results clearly show that the half quantization appears in a narrow but finite range of the magnetic field. These results strongly suggest the existence of the chiral Majorana edge mode in the topologically-nontrivial state in the field-induced PM region. We note, however, that the T dependence of xy =T is different from the theoretical results in Fig. 44 both quantitatively and qualitatively. The experimental data exhibits the overshoot above the half quantization value, which is not obtained in the theoretical results. Moreover, the set-in temperature of the half quantization is considerably high compared to the theoretical prediction: The asymptotic convergence in theory appears well below T L , which roughly corresponds to ∼1 K, as shown in Fig. 44. One of the reasons for such discrepancies is that the theoretical results were obtained for the effective model which is justified in the weak-field limit. More sophisticated theory beyond the perturbation is highly desired. Furthermore, non-Kitaev interactions may play an important role in the topological phenomena. Indeed, it was pointed out that a symmetric offdiagonal interaction contributes to the stabilization of the gapped topological state. 72,196) Another caveat is the contribution from phonons. The large value of the longitudinal thermal conductivity xx at low T suggests the dominant phonon contribution. 197) The possibility of the observation of quantized xy even in such a situation was theoretically discussed. 198,199) Figure 47 summarizes the field-T phase diagram elaborated by the experiments. In the field region between ∼7 and ∼9 T after the magnetic order is suppressed (red area in Fig. 47), the half quantization of xy =T is observed below ∼5 K. This is the region where the Majorana topological state is suggested to be realized. Thus, the results offer strong evidence of the Kitaev-type QSL with a gapped excitation in the field-induced PM state [see also the schematic phase diagram in Fig. 14(a) in Sect. 2.9].
The half quantization of xy =T has attracted great attention since it can be regarded as the direct evidence of the spin fractionalization in the Kitaev system, especially the Majorana fermionic nature. Moreover, it is intriguing as the set-in temperature is rather high compared to other topological phenomena like the quantized anomalous Hall effect in magnetic topological insulators. 200,201) Theoretically, it was pointed out that nonabelian anyons emergent in the topological state can be utilized for fault-tolerant quantum computation. 16,17) Thus, the experimental finding of the half quantization may offer a first step toward topological quantum computing based on peculiar quasiparticles in magnets.

Summary and Perspectives
In this article, we have overviewed the recent development in the research of the Kitaev quantum spin liquids and their experimental realization. We have reviewed finite-T properties of the Kitaev model, including the spin dynamics, which have been revealed by the Majorana-based numerical techniques developed by the authors and their collaborators. In the Kitaev model, the spin degree of freedom is fractionalized into two different types of quasiparticles: itinerant Majorana fermions and localized Z 2 fluxes. They have largely different energy scales and affect the thermodynamics and spin dynamics in a peculiar manner, which we call thermal fractionalization. We have discussed a number of fingerprints of the thermal fractionalization in experimentally observable quantities, and compare them with available experimental data for the candidate materials. Let us summarize the main points, focusing on the 2D honeycomb case: • The Kitaev model exhibits two characteristic temperatures corresponding to the two types of quasiparticle excitations. They define two crossovers at largely different temperatures T H and T L (T H ) T L ), signaled by two broad peaks in the specific heat and corresponding successive releases of the entropy by half ln 2 (Sect. 3.1.1). Similar behavior corresponding to the high-T crossover was observed experimentally in candidate materials, Na 2 IrO 3 , α-Li 2 IrO 3 , and α-RuCl 3 (Sect. 5.1). • The high-T crossover at T ¼ T H is caused by the itinerant Majorana fermions, while the low-T one at T ¼ T L is by the localized Z 2 fluxes. The temperature scales T H and T L are set by the center of mass of the density of states for the complex fermion band and the Z 2 flux gap, respectively (Sect. 3.1.2). • The two crossovers define three distinct regimes: the conventional paramagnetic state for T ≳ T H , the frac-  tional paramagnetic state for T L ≲ T ≲ T H , and the asymptotic quantum spin liquid state for T ≲ T L [see Fig. 27 Weiss law below T $ J (J is the Kitaev coupling), shows a peak in the intermediate-T region between T L and T H . Similar behaviors were observed for the candidate materials. Theoretically, the susceptibility shows a rapid decrease around T ¼ T L and approaches a nonzero value in the low-T limit, but these behaviors are hindered by the magnetic ordering in the real compounds (Sect. 5.3). • The dynamical spin structure factor Sðq; !Þ shows a characteristic T dependence. Below T ' T H , Sðq; !Þ develops an incoherent response at ! ' J with less q dependence, which persists down to the lowest T. The overall T, q, and ω dependences agree well with the experimental data by inelastic neutron scattering for α-RuCl 3 . Theoretically, while approaching T L , an additional quasielastic response grows rapidly, and it is gapped out below T L reflecting the gap opening in the flux excitations. But, these behaviors are not observed in experiments due to the magnetic ordering (Sect. 5.4). • The NMR relaxation rate 1=T 1 increases below T H , where the static spin correlations are almost saturated. This dichotomy between dynamical and static spin correlations is a possible indication of the thermal fractionalization. 1=T 1 exhibits a broad peak above T L , and decreases exponentially below T L reflecting the flux gap opening. Similar behaviors were observed in experiments under a magnetic field, suggesting the potential realization of the Kitaev quantum spin liquid in the field-induced paramagnetic state (Sect. 5.5 We have also discussed interesting signatures of the thermal fractionalization for the Kitaev models with some extensions from the original honeycomb one. There appear a variety of phase transitions and crossovers, as schematically summarized in Fig. 27 in Sect. 3.4. We list the key aspects in the following, which await for the experimental confirmation: • In the 3D Kitaev model, the nature of the Z 2 flux excitations is qualitatively different from that in two dimensions. Because of the local constraint on the Z 2 fluxes, the excitations are allowed only in the form of closed loops in three dimensions. This changes the low-T crossover in the 2D cases into a phase transition. This transition takes place between the high-T paramagnet and the low-T quantum spin liquid, which can be regarded as a gas-liquid transition in terms of the spin degree of freedom of insulating magnets (Sect. 3.3.1). • When extending the Kitaev model by adding non-Kitaev interactions, the system may undergo phase transitions among three states of mattergas, liquid, and solid. The phase diagram is distinct between the 2D and 3D cases, reflecting the different nature of the Z 2 flux excitations (Sect. 3.3.2). • When the Kitaev model is defined on the lattice structures with odd-site loops, the ground state can be a chiral spin liquid. In this case, the low-T crossover is replaced by a finite-T phase transition with breaking of time-reversal symmetry caused by Z 2 flux ordering (Sects. 3.2 and 3.3.3). Despite the clarification of many intriguing aspects of the thermal fractionalization in the Kitaev model and the successful comparison with experimental data, there remain a number of open issues in this rapidly growing field. We hope that the present review will be helpful for studying the following issues in future studies.
• Further theoretical understanding of the Kitaev model and its extensions: -It is highly desired to clarify the effect of the external magnetic field on the phase diagram, the topological properties of the elementary excitations, and the excitation spectra. This is crucially important for comparison with experimental data, especially the remarkable properties discovered in the field-induced paramagnetic state in α-RuCl 3 . -In the magnetic field, the case of the antiferromagnetic Kitaev coupling is also intriguing, since an additional topological phase was predicted theoretically, as discussed in Sect. 2.7. It is also important to find the candidate materials for the antiferromagnetic Kitaev coupling, by pushing forward the recent efforts introduced in Sect. 2.7. -Development of new theoretical techniques is a key to breakthrough in understanding of the effects of the magnetic field and non-Kitaev interactions listed above. -It is also important to clarify the effects of non-Kitaev interactions which exist in real compounds, as mentioned in Sect. 2.8. In particular, it is crucial to study such effects in sufficiently large system sizes with high resolution in both energy and momentum, since the subdominant interactions can lead to keen competitions between different phases and fine structures in the excitation spectra. -It is worth extending the analyses to other lattices, especially in three dimensions. Besides the hyperhoneycomb, hyperoctagon, and hypernonagon structures discussed in Sect. 3.3, a variety of extensions were discussed for other lattices. 202) Interestingly, depending on the underlying lattice structures, the itinerant Majorana fermions form Majorana Fermi surfaces, nodal lines, or topologically-protected Weyl nodes. In addition, the Z 2 flux configurations can be suffered from frustration. 203) A comprehensive study of finite-T properties for such extensions will deepen our understanding of the Kitaev quantum spin liquids and fractionalization. -It would also be interesting to consider extensions of the Kitaev model to larger spins. The local conserved quantity on each plaquette exists also in the larger spin cases. 45) Recently, thermodynamic properties were studied numerically. [204][205][206][207] While the realization of such systems was theoretically proposed, 208) the search for the candidate materials has just begun. 209,210) • Further quantitative comparison with experiments: -Further experimental identification of fractional quasiparticles is an important issue. In particular, the Z 2 flux excitations have not been identified clearly thus far. It would be helpful to further study the fieldinduced paramagnetic state in α-RuCl 3 at lower T. -Regarding the potential topological quantum spin liquid in the field-induced paramagnetic state, the crucial questions are what kind of the gap protects the topological state, how large the gap is, and how it depends on the field. Extensive experiments have been done for the excitation gap in the magnetic field, for instance, the specific heat, [156][157][158] NMR, [180][181][182][183] electron spin resonance, 211) terahertz spectroscopy, [212][213][214] inelastic neutron scattering, 175,176) and thermal conductivity measurements, 197,215,216) but there still remains controversy, even among the results obtained by the same experimental probes. Although the theoretical study in the field is also very difficult, close comparison between experiment and theory on this gap issue will be crucial to deeper understanding of the field-induced state. -It is important to precisely estimate the additional non-Kitaev interactions for each candidate material by further comparison between theory and experiment. This issue has been addressed by the analyses of, e.g., the magnon spectra in the ordered phases. 217,218) The gap problem above would also be helpful to this issue. Also, further detailed analysis on the magnetic anisotropy would play an important role, as stated in Sect. 5.3.
-It is also important to discuss the effect of disorder, which is inevitably present in real compounds, on the physical observables at finite T. This includes nonmagnetic=magnetic impurities, [219][220][221] dislocations, 222) chemical inhomogeneity, and so on. • Coupling to other degrees of freedom: -Given the fractional quasiparticles, it will be very interesting to consider the coupling to other degrees of freedom, for instance, the electric charge. The dynamics of a single hole doped into the Kitaev quantum spin liquid was studied. 223,224) It was also predicted that carrier doping to the Kitaev model and its extensions may lead to topological superconductivity, reflecting the exotic nature of the Kitaev quantum spin liquid. [225][226][227] Theoretical studies beyond the mean-field calculations as well as the experimental realization are highly desired. -It will also be intriguing to study the proximity effect to other magnets, metals, and superconductors. Recent development in the heterostructure of α-RuCl 3 and graphene, which was introduced in Sect. 4.2, is a good example in this direction. The coupling between the fractional quasiparticles and other degrees of freedom, such as mobile electrons, Cooper pairs, magnons, and phonons, may lead to unprecedented physics. Indeed, the coupling to mobile electrons was discussed for the Kitaev-Kondo model, and topological superconductivity was predicted. 228,229) In addition, effects of lattice strain are also worth investigating as a source of exotic states. 230,231) • Further materialization of Kitaev quantum spin liquids: -As partly reviewed in this article, the candidate materials for the Kitaev spin liquids are still limited. Further exploration is needed. In particular, highly desired are candidates which show the Kitaev spin liquid nature at zero or weaker magnetic field. Materials with the AFM Kitaev coupling are also desired, as mentioned above. -Material design for new lattice structures is important.
In particular, 3D materials are desired for studying the intriguing physics listed above. Interesting proposals were made by using metal organic frameworks. 232,233) In addition, quasi-one-dimensional candidates, e.g., with a ladder structure, are also interesting to further clarify the nature of fractional quasiparticles. -It would also be important to explore candidates in the form of thin films and heterostructures, especially for studies of the proximity effects mentioned above. • Control of fractional quasiparticles: -In the topologically nontrivial phase under the magnetic field, each excited flux in the bulk accompanies a Majorana zero mode, which obeys nonabelian statistics [see Fig. 12(b)]. Toward topological quantum computation by using the nonabelian anyons, it is a crucial task to invent a way for controlling them, e.g., braiding and fusion. A potential way will be to use local geometry of the system, such as defects, dislocations, edges, and interfaces. Another way would be local perturbations, e.g., by using the scanning tunneling microscope. -Along this direction, it will be quite important to clarify nonequilibrium dynamics of the fractional quasiparticles, as the topological quantum computing will be implemented by the time evolution of the quasiparticles. Although there were several attempts for clarifying the nonequilibrium dynamics by theory [234][235][236][237][238][239][240][241][242] and also in experiments, [243][244][245][246] but further studies are desired.
where A f r g is an N Â N Hermite matrix with pure imaginary matrix elements, and therefore, A f r g ij ¼ ÀA f r g ji . This is diagonalized as where f y and f are the creation and annihilation operators of the complex fermion corresponding to the energy E f r g (> 0), and E f r g is the ground-state energy. Here and hereafter, the sum P is taken for positive energies ( ¼ 1; 2; . . . ; N=2) although all the eigenvalues of A f r g appear in pairs as AE 1 2 E . The complex fermions ff g are introduced such that where U f r g j is the jth component of the eigenvector associated with the eigenvalue 1 2 E f r g of the matrix A f r g . To calculate thermodynamic quantities, we introduce the partition function by Tr f i g e ÀH frg : ðA : 4Þ This is rewritten as ðA : 5Þ where F f r g is the free energy of the Majorana fermion system for the configuration of f r g, which is given by ln½Tr f i g e ÀH frg : ðA : 6Þ Using the eigenvalues of the matrix A f r g , the partition function of the Majorana fermion system is evaluated as Similar to the Hamiltonian, an operator commuting with all f r g can be labeled by f r g as O f r g . The thermal average of such an operator can be calculated by Tr f i g ½Oe ÀH fr g ¼ h " O f r g i ; ðA : 8Þ where we introduce the expectation value of O for the configuration of f r g as On the other hand, one cannot straightforwardly calculate thermal averages of the operators not commuting with f r g, such as dynamical spin correlations. We will introduce a way to calculate such quantities in Appendix A.3. Using Eqs. (A·5) and (A·10), finite-T properties of the Kitaev model can be calculated by using the MC sampling on the configurations of f r g. At a certain temperature, we calculate the free energy F f r g and " O f r g for a given configuration f r g in a finite-size cluster by exact diagonalization of the Hermite matrix A f r g . Using the Markov-chain MC simulation, the sequence ðf r g 1 ; f r g 2 ; f r g 3 ; . . . ; f r g N MC Þ is successively generated so as to reproduce the probability distribution e ÀF frg f =Z. In the sequence of f r g, the thermal average of an operator O is evaluated by replacing hÁ Á Ái by hÁ Á Ái MC as In Sects. 3.1.1 and 3.3.1, this technique is applied to calculate the internal energy, specific heat, entropy per site, and the DOS for the complex fermion band. The internal energy per site is calculated as where : ðA : 13Þ The specific heat per site can also be calculated as ðA : 14Þ From the specific heat, the entropy per site is obtained as In addition, the contributions from itinerant Majorana fermions and localized Z 2 fluxes are separately calculated as ðA : 16Þ respectively. The corresponding contributions to the entropy are calculated in a similar manner to Eq. (A·15). The fermion DOS is computed by ðA : 18Þ which depends on temperature T. Using this expression, E and C v are written as The same method is applied to compute the thermal conductivity and the Raman scattering intensity in Sects. 5.6 and 5.7, respectively. These are feasible as the thermal current operator and the Raman operator commute with all f r g. In Sect. 5.8, the thermal Hall conductivity is calculated in the same manner, but in this case, for the Hamiltonian including the effect of the Zeeman coupling effectively in Eq. (23). For this effective Hamiltonian, f r g are still conserved and the thermal current operator commutes with f r g.

A.2 Cluster dynamical mean-field theory
In the Majorana representation, one can also apply the CDMFT, which has been developed for interacting fermion systems. 247) In the case of the Kitaev model, the system can be regarded as a noninteracting fermion system coupled with localized classical variables, similar to the Falicov-Kimball and the double-exchange models, as mentioned in Sect. 2.6. For this category of the models, the impurity problem in the where pðf r gÞ is the weight of the configuration of f r g, which is given by ðA : 32Þ Z f r g is calculated from Green's functions as det½ÀG f r g ði! n Þ: ðA : 33Þ Finally, the self-energy for the impurity problem is obtained as