Topological Nematic Phase Transition in Kitaev Magnets Under Applied Magnetic Fields

We propose a scenario of realizing the toric code phase, which can be potentially utilized for fault-tolerant quantum computation, in candidate materials of Kitaev magnets. It is demonstrated that four-body interactions among Majorana fermions in the Kitaev spin liquid state, which are induced by applied magnetic fields as well as non-Kitaev-type exchange interactions, trigger a nematic phase transition of Majorana bonds without magnetic orders, accompanying the change of the Chern number from $\pm 1$ to zero. This gapful spin liquid state with zero Chern number is nothing but the toric code phase. Our result potentially explains the topological nematic transition recently observed in $\alpha$-RuCl$_3$ via heat capacity measurements [O. Tanaka et al., arXiv:2007.06757].


I. INTRODUCTION
A quantum spin liquid (QSL) is an exotic phase of matter, where spins do not have long-range order even at zero temperature. The Kitaev model [1] is an exactly solvable model of QSLs, which is defined on the honeycomb lattice with bond-dependent Ising interactions. A remarkable feature of the Kitaev model is that the system is described by Majorana fermions interacting with Z 2 gauge fields, allowing the realization of Abelian and non-Abelian anyons. For almost isotropic bond-dependent Ising interactions, the systems exhibits a chiral QSL phase with non-Abelian anyons, when an applied magnetic field opens an energy gap of Majorana fermions. On the other hand, for highly anisotropic Ising interactions, the toric code phase with Abelian anyons is realized. It is noted that both of these phases can be utilized for topological quantum computation [1,2]. After the Kitaev's seminal paper, various phenomena associated with Majorana fermions in Kitaev magnets have been extensively explored [3][4][5][6][7][8][9][10][11][12][13][14][15][16]. Since Jackeli and Khaliullin [17][18][19][20][21][22] pointed out that isotropic bond-dependent Ising interactions arise in some honeycomb layered materials, the experimental search for candidate materials has been an important subject. Several materials including 4d 5 or 5d 5 transition metal ions with a strong spin-orbit coupling have been proposed for Kitaev magnets, such as Na 2 IrO 3 [23,24], α-Li 2 IrO 3 [25], and α-RuCl 3 [26]. α-RuCl 3 under an applied magnetic field is one of the best candidates for the Kitaev magnet [27]. For this material, a signature of Majorana fermions in the chiral QSL phase is observed via the measurement of the half-integer thermal quantum Hall effect [28][29][30]. Moreover, according to the recent specific heat measurements [31], approximated threefold rotational symmetry of the honeycomb lattice is drastically broken to obviously twofold one at high fields, implying that nematic-type phase transition occurs [32]. Remarkably, the nematic phase transition accompanies * takahashi@blade.mp.es.osaka-u.ac.jp the disappearance of the half-quantized thermal Hall effect [28], though experimental results suggest that the system may be still in a QSL phase. To understand these observations, we need another ingredient missing in previous research for Kitaev magnets, which includes various additional spin-spin interactions, such as Heisenberg, Γ, and Γ terms.
To focus on the nature of a QSL phase in Kitaev magnets from a different point of view, we assume the fluxfree state in our model. Although the flux-free assumption may be violated by strong additional spin-exchanges like Heisenberg, Γ, and Γ terms, it allows us not only to clarify properties of a QSL phase directly but also to explore topological phase transition mentioned above. Besides, under the flux-free assumption we can combine the mean-field analysis and exact diagonalization, which is impossible in the previous framework of the Majorana mean-field theory (see Refs. [55,64], for example), so that we can study rotational symmetry breaking in the extended Kitaev model from various aspects by computing the Majorana band structure, the Majorana bond order, and the Chern number, which are hard to be handled by the direct treatment of the spin model. Furthermore, arXiv:2101.05959v3 [cond-mat.str-el] 23 May 2021 the scenario that the Majorana four-body interactions give rise to rotational symmetry breaking is applicable to any flux sectors, which is described by interacting Majorana systems, and thus, it is a good starting point to focus on one flux sector. From now on, we elucidate that the four-body interactions can induce nematic-type Majorana bond ordered phase which is a gapped QSL with zero Chern number, i.e. the toric code phase.

II. MODEL
The Hamiltonian of the pure Kitaev model is expressed as, where σ i represents the spin operator on the ith site and ij means a nearest-neighbor (NN) bond. x, y, and z are indices for the bond direction. The phase diagram in the parameter space (J x , J y , J z ) obtained in Kitaev's original paper is shown in Fig. 1(a). For strongly anisotropic Kitaev interactions, the so-called A phase, which is the toric code phase, is realized. On the other hand, for almost isotropic Kitaev interactions, another phase (B phase), which possesses gapless Dirac bands of Majorana fermions, appears. Once the time-reversal symmetry is broken in B phase, however, the spectrum acquires an energy gap, leading to a chiral QSL state with the Chern number ±1. For both phases, the system is described by Majorana fermions interacting with Z 2 gauge fields. In the following, we restrict our analysis within the vortexfree states. Then, in the case with applied magnetic fields, from the perturbation theory in the third order, effective interactions of three-spin operators are obtained: where h = (h x , h y , h z ) is an external magnetic field and we here assume that J x = J y = J z = J. The site summation i, j, k is taken as shown in Fig. 1(b). In the Majorana representation, operators are rearranged to three types of terms: where "→" means taking the standard gauge in the 0-flux ground state, and ij means a next-nearest-neighbor (NNN) bond. The first term is an NNN hopping that gives rise to an energy gap, leading to the Chern number ±1, and the others are four-Majorana interactions.
Here the site summation for Y is taken as shown in the upper Y-shaped diagram of Fig. 1(c) and Y is taken for the lower one of Fig. 1(c). In the following, we focus on these types of Majorana interactions for concreteness.
However, we believe that our main results are generic also for other types of neighboring four-body interactions. It is noted that in the case with the off-diagonal exchange interaction, where α, β run for x, y, and z, the Y(Y')-shaped interactions also arise from the third-order perturbation terms of order O(h x,y,z Γ 2 ). Besides, other non-Kitaev interactions, such as Heisenberg and Γ terms, contribute to bring additional NN and NNN hoppings as shown in Appendix C. Thus, generally, the coefficients of the second and the third terms of Eq. (3) are different from that of the first one, and we end up with the following Majorana Hamiltonian, with We, here, stress again that the parameters κ and g are independent for real candidate materials of Kitaev magnets, as mentioned above. The Majorana operator is represented as c i for the ith site, and the operators obey The hopping amplitude t = J is set unity without loss of generality.

III. MEAN-FIELD ANALYSIS
First, we apply the mean-field (MF) analysis to the four-body interaction term H 4 . Then, the general form of the MF Hamiltonian is, where the index a = 1, . . . , 6 specifies the directions of hopping with the amplitude τ a , as described in Fig. 1(d).
The phase factors η ij = ±1 is necessary to make H MF antisymmetric. The sign of η ij is defined as illustrated in Fig. 1(d) such that the direction of the arrow indicates the positive hopping from the jth site to the ith site. Utilizing the Hellmann-Feynman theorem, we can derive self-consistent equations as follows (see Appendix A for more details): where the definition of bond order parameters is with |Ψ MF being the ground state of the MF Hamiltonian.
We solve the MF Hamiltonian by numerical iteration. First, we focus on a nematic transition. The bond order parameters ∆ a (a = 1, . . . , 6) play an important role for nematic transitions with rotational symmetry breaking. The Majorana nematic phase considered here is defined as a phase with a bond order in a specific direction. We define nematic order parameters as [65], If ∆ 1 = ∆ 2 = ∆ 3 (∆ 4 = ∆ 5 = ∆ 6 ), φ (ψ) must be 0 from the definition, but otherwise will have a nonzero value. In other words, φ and ψ characterizes breaking of the threefold rotational symmetry. As seen in Figs. 2(a) and (c), for large value of |g|, a nematic phase with |φ| = 0, |ψ| = 0, characterized by two-fold rotational symmetry, appears. We also find in Fig. 2(b) that arg(φ) is π/3 in most regions of the |φ| = 0 area. The argument of φ has informations about the direction of the nematic order. Since ∆ a takes a negative value in this calculation, arg(φ) = π/3 corresponds to the nematic order with |∆ 3 | > |∆ 1 |, |∆ 2 |. On the other hand, as seen in Fig. 2(d), arg(ψ) = π/3 for g < 0 in most of the nematic region, while arg(ψ) = −2π/3 for g > 0. This is because that ψ → −ψ under time-reversal operation, while φ → φ. In the calculations, we intentionally introduce tiny anisotropy of bonds to obtain symmetrybreaking solutions, lifting three-fold degeneracy. Now, we examine topological features of the nematic phase. For this purpose, we calculate the Chern number numerically using the Fukui-Hatsugai-Suzuki method [66]. The results are shown in Fig. 2(e). In the chiral QSL phase without the nematic order, the Chern number ν = ±1. Remarkably, on the other hand, in the nematic phase with |φ| = 0 and |ψ| = 0, we find ν = 0. Thus, the nematic phase transition accompanies the topological phase transition. This implies that the nematic transition which induces strong anisotropy of the Majorana hopping terms (τ 1 , τ 2 , τ 3 in Eq. (6)) drives the system from B phase to A phase in the phase diagram shown in Fig.1(a). In fact, in this nematic phase, since Z 2 vortices are still suppressed, the system has no magnetic long-range order, and is in the QSL state. Moreover, the nematic phase is described by the quadratic Hamiltonian of Majorana fermions Eq. (6) with an energy gap. According to Kitaev's argument on the 16-way of anyons composed of Z 2 vortices and free fermions [1], the gapful Majorana nematic phase with the Chern number ν = 0 is identified with the toric code phase with Abelian anyons. It is noted that this nematic phase transition is the firstorder type with the discontinuous changes of the order parameters φ and ψ, and hence, there is no gap closing at the transition point. We stress that although the nonzero ψ apparently causes the anisotropy of the Majorana hopping terms, both of the two nematic orders φ and ψ cooperatively stabilize the toric code phase, since φ and ψ are nonlinearly coupled with each other via the MF equations. The details are shown in Appendix B.

IV. EXACT DIAGONALIZATION
To confirm the results of the MF calculation, we employ the exact diagonalization calculations. We apply the Lanczos method to the system size up to 32 sites. First, we discuss the case of 18 sites. When we use periodic boundary conditions (PBCs), φ and ψ must be zero, because the system completely preserves the rotational symmetry for finite size systems. Thus, instead, we consider the fluctuation of the order parameters {φ, φ † }/2 and {ψ, ψ † }/2 which can have a nonzero value even for the PBC. The calculated results are shown in Figs. 3(a) and (g). There are two bright yellow regions, which implies that the nematic phase transition occurs discontinuously for large values of g. This result confirms the MF calculation. To examine the nematic order more directly, and to determine the direction of bond orders, we switch to an antiperiodic boundary condition (APBC) to intentionally break the threefold rotational symmetry so that we can compute the nonzero |φ| and |ψ|, and the arguments of φ and ψ. The nematic phase transition is verified again as seen in Figs. 3(c) and (i). Moreover, another phase boundary appears within the nematic regions. This implies that there are two different nematic phases. One is a nematic phase which has one strong bond so that arg(φ) = π + 2nπ/3 (n ∈ Z), and the other is a zigzag nematic phase with two strong bonds so that arg(φ) = 2nπ/3 (n ∈ Z) (see Fig. 3(e)). However, the shape of the phase boundary between two nematic phases is quite different from that of the MF result shown in Figs. 2(b) and (d). In fact, this phase boundary crucially depends on the system size, as seen below.
The results of the 32-site calculation are shown in Figs. 3(b), (d), (f), (h), (j), and (l). Although the region of the nematic phase depicted by orange color becomes narrower compared to the 18-site calculation, the nematic phase is definitely realized for large values of |g|. We also find two types of nematic phases as in the MF calculations and the 18-site calculations. However, the shape of the phase boundary is quite different. The region of the zigzag nematic phase with the Chern num- , becomes very small for the 32-site calculation. To clarify the fate of the zigzag nematic phase, we need calculations for larger system sizes. Although there is a significant finite-size effect in the calculations, we can clearly observe the tendency towards the nematic phase transition. In the end, small non-zero values in non-nematic phases seen in (c), (d), (i), and (j) are finite-size effects, and decrease toward zero as the system size increases.

V. DISCUSSION
Using the MF analysis and the exact diagonalization method, we find that the topological nematic phase transition occurs in strongly correlated regions |g| ∼ t, g ∼ κ. Since non-Kitaev interactions such as off-diagonal exchange interactions significantly renormalize t, κ, and g, it is possible to realize this condition for the extended Kitaev model. Particularly, the Γ term reduces the hopping amplitude t, and drives the system into strongly correlated regions. We note that the Y(Y')-shaped interactions of order ∼ O(h x,y,z Γ 2 ) cancel out for magnetic fields parallel to the honeycomb plane. Nevertheless, other four-body interactions with different configurations of Majorana fermions, which do not disappear even for parallel magnetic fields, can be generated from the Heisenberg exchange interaction and the other offdiagonal exchange interaction, i.e. the Γ term. The details of pertubative calculations with respect to non-Kitaev interactions are shown in Appendix C and we expect that the nematic phase transition recently observed for α-RuCl 3 under parallel magnetic fields [31] may be explained by taking into account these contributions. To establish this scenario, we need further both theoretical and experimental investigations. It is an interesting future issue to extend our analysis to general flux configurations, for which various Majorana phases are predicted [67,68]. From an experimental side, an electric probe through the hyperfine interaction [69] may be useful for the detection of the Majorana nematic phase.

VI. SUMMARY
We have demonstrated that, in the Kitaev spin liquid state, four-body interactions among Majorana fermions which are induced by applied magnetic fields and non-Kitaev exchange interactions give rise to the topological nematic phase transition from the chiral QSL phase to the toric code phase.
In this section, we reveal the relation between the MF hopping strength τ a and bond order ∆ a . The Fourier transform of a Majorana operator is given by, where N unit is the total number of unit cells and j represents an each site. λ is a position type inside the unit cell. In other words, inverse Fourier transform is expressed as to satisfy c † k,λ = c −k,λ and {c k,λ , c † q,µ } = δ k,q δ λ,ν . Using δ(k) = j e irj ·k /N unit , we can calculate each term of the MF Hamiltonian as, where we choose a basis (n 1 , n 2 ) of the translation group used in Ref. [1] (see Fig. 4). So the MF Hamiltonian in where D 2 (k) = τ 4 sin(−n 2 · k) + τ 5 sin(n 1 · k) +τ 6 sin((n 2 − n 1 ) · k). (A6) Then the energy spectrum and the occupied energy under the Fermi level are written as where the summation k is taken over the first Brillouin zone. Then, we evaluate the mean field energy by using the MF ground state |Ψ MF and the MF Hamiltonian H MF , where N is the total number of sites, and we use, respectively. Utilizing the Hellmann-Feynman theorem, we obtain, ∂E MF ∂τ a (a = 4, 5, 6).

(A10)
We can explicitly perform the derivatives using Eq. (A8) and get a set of equations, On the other hand, the total ground state energy which we will take variations with respect to {τ a } can be easily written as and one can express the first term and the second term with ∆ a as In addition, applying Wick's theorem, the final term of Eq. (A11) can be written as In short, the total energy density can be expressed as, Now, we can obtain the self-consistent equations by minimizing the energy with respect to variational parameters {τ a }, or more explicitly, Finally, comparing Eq. (A10) and Eq. (A14), we end up with, (A15) In this section, we discuss how to distinguish between a nematic phase and a zigzag nematic phase in strongly interacting regions. Since ∆ a (a = 1, ..., 6) has a negative value in most of parameter regions (see Fig. 5(b)), the nematic order of φ is schematically expressed as shown in Fig. 5(a); arg(φ) reflects the direction of bond orders with |φ| = 0. For instance, arg(φ) = π/3 implies |∆ 1 | = |∆ 2 | < |∆ 3 |, which is a direct signal of a nematic phase with one strong bond along the a = 3 direction. On the other hand, if arg(φ) = −2π/3, there are two strong bonds in a = 1, 2 directions. In general, arg(φ) = π+2nπ/3 (n ∈ Z) represents a nematic phase in a one-strong-bond order and arg(φ) = 2nπ/3 (n ∈ Z) indicates a zigzag nematic phase that has two strong bonds (see Fig. 5(c)). Furthermore, Fig. 5(b) shows that the symmetry breakings of φ and ψ occur at the same transition point, where we intendedly realize anisotropy of bond orders by substituting several patterns of initial values and take the lowest energy state in each calculation. We also note symmetrical properties of φ and ψ under time reversal operation. Since time-reversal operation for Majorana operators depends on the sublattice degrees of freedom, ∆ a for a = 1, 2, 3 preserves the time-reversal symmetry, while ∆ a for a = 4, 5, 6 changes its sign under time-reversal operation. Therefore, the time-reversal operation acts on φ and ψ as described below: Because of these symmetrical properties, arg(φ) is not changed by the transformation g → −g and κ → −κ, which correspond to the time-reversal operation, while arg(ψ) changes its value as π/3 ↔ −2π/3 under this transformation, as seen in Figs. 2 and 3 in the main text. In short, for g < 0, the region with arg(ψ) = π/3 (arg(ψ) = −2π/3) corresponds to ∆ 4 = ∆ 5 < ∆ 6 (∆ 4 = ∆ 5 > ∆ 6 ), and for g > 0, the region with arg(ψ) = −2π/3 (arg(ψ) = π/3) corresponds to ∆ 4 = ∆ 5 < ∆ 6 (∆ 4 = ∆ 5 > ∆ 6 ).

Appendix C: DERIVATION OF FOUR-MAJORANA INTERACTIONS ARISING FROM NON-KITAEV INTERACTIONS
In this section, we present the details of the derivation of the four-Majorana interactions caused by non-Kitaev interactions on the basis of perturbative calculations. In these calculations, we need to evaluate matrix elements of perturbation terms for the eigenstates of Z 2 vortices. For this purpose, we, first, examine the operation of spin operators on the eigenstate, and then, explore for nonzero matrix elements of perturbation terms generated by symmetric off-diagonal exchange interactions, i.e. the Γ term, and the Γ term, and also by the Heisenberg exchange interaction. It is found that the Y-shaped four-Majorana interactions are generated by the Γ term combined with the Zeeman term. On the other hand, the Γ term and the Heisenberg interaction give rise to armchair-shaped interaction, and zigzag-shaped interactions, respectively (see below). Up to the third-order perturbation, this analysis exhausts all four-body interactions among neighboring Majorana fermions generated by the combination of applied magnetic fields and the Heisenberg or the symmetric off-diagonal exchange interactions.

Flip of Z2 gauge fields due to the operation of spin operators
In the following, we use the Majorana representation of s = 1/2 spin operators introduced by Kitaev [1], The pure Kitaev model is expressed in terms of itinerant Majorana fields c j interacting with Z 2 gauge fieldsû α jk = ib α j b α k on αbonds (α = x, y, z). Here, we consider effects of the operation of spin operators on the eigenstate of flux operators, w, where |Ψ w and |Ψ u are, respectively, the eigenstates of flux operators, and that of Z 2 gauge fieldsû, and We also define a state |φ that is an excited state with a flipped Z 2 gauge field. Note that the following relations hold: Using these relations, one can see that |φ is actually a state with a flipped eigenvalue ofû z : This means that the sign of the eigenvalue of the Z 2 gauge field for the |φ state is opposite to that of the |Ψ u state. The flip of the Z 2 gauge field occurs also in the case that the spin operator σ z acts on the other site of the z-bond,

Perturbative calculations with respect to non-Kitaev interactions
We start with the following Hamiltonian for candidate materials of the Kitaev magnet on a honeycomb lattice such as α-RuCl 3 and Na 2 IrO 3 [18,19,21,70], where σ α i is an α = x, y, z component of an s = 1/2 spin operator at a site i. H K is the Kitaev interaction, H J H is the Heisenberg exchange interaction between the nearest neighbor sites, and H Γ and H Γ are symmetric off-diagonal exchange interactions. Here, ij α denotes that the ith site and the jth site are connected via a nearest-neighbor α-bond on the honeycomb lattice. We note that the definitions of J H , Γ, and Γ are different from conventional ones used in first principle calculations [19,21,70] by a factor 1/4. We also take into account the Zeeman term due to an external magnetic field h = (h x , h y , h z ), Putting V = H J H + H Γ + H Γ + H Z , we carry out the perturbative expansions with respect to V around the vortex-free ground state in the same spirit as the Kitaev's paper [1]. In this perturbation analysis, intermediate excited states have vortex excitations (visons) with finite energy gaps. On the basis of the results derived in the previous section, we examine nonzero matrix elements in each perturbation term which generates four-Majoarana interactions.
We, first, consider the Y-shaped interaction, which is derived from the third-order perturbation with respect to the Γ term, H Γ , and the Zeeman term H Z . To be concrete, we show an example of the perturbation processes in Fig. 6. In this example of the third-order perturbations, the operations of Γ σ z 1 σ x 2 and Γ σ x 2 σ y 3 and h x σ x 4 on the vortex-free ground state result in the final state, which is also the vortex-free ground state. Thus, these perturbation processes are allowed within the ground state sector. On the other hand, H Γ and H J H do not give the perturbation processes within the ground state sector which result in the Y-shaped interaction. Thus, we, here, omit these two terms in V . Then, the third-order perturbation term which leads to the Y-shaped interaction is given by, where α, β, γ = x, y, z, and Π 0 is a projection to the vortex-free spin liquid state. Also, ∆ is the energy gap of two visons, which is given by ∆ ∼ 0.26J [1] for the configurations shown in Figs. 6(b), (c). To analyze this term more precisely, we use the fact that in the Majorana fermion representation of the Kitaev spin liquid state, gauge Majorana fields b α i should be paired on the α-bond connecting two sites i and j to form Z 2 gauge fieldsû α ij = ib α i b α j , since the Kitaev spin liquid state is expressed by the eigen state of the Z 2 gauge fields. Then in the case of α = x, Eq. (C15) is recast into, This amounts to the Y-shaped interaction of four itinerant Majorana fermions on the sites i, j, k, l shown in Fig. 6(e). Secondly, we consider the armchair-shaped interaction, which is generated by the third-order perturbation with respect to the Γ term and the Zeeman term. An example of the perturbation processes is shown in Fig. 7. In this example of the third-order perturbations, the operations of Γσ x 1 σ y 2 and h y σ y 3 and h x σ x 4 on the vortex-free ground state generate the final state, which is also the vortex-free ground state. Thus, these perturbation processes are allowed within the ground state sector. On the other hand, H Γ and H J H do not generate perturbation processes which lead to the the armchair-shaped interaction satisfying the above condition. Thus, we omit H Γ and H J H in V in this perturbative calculation. Then, the third-order perturbation term is given by, where ∆ ∼ 0.23J is the energy gap of two visons for the configuration shown in Fig. 7 This results in the armchair-shaped interaction of four neighboring Majorana fermions on the sites i, j, k, l shown in Fig. 8(e). Finally, we consider the zigzag-shaped interaction, which is generated by the third-order perturbation with respect to the Heisenberg term H J H and the Zeeman term H Z . An example of the perturbation processes is shown in Fig. 8  on the vortex-free ground state generate the final state, which is also the vortex-free ground state. Thus, these perturbation processes are allowed within the ground state sector. On the other hand, H Γ and H Γ do not contribute to the generation of the zigzag-shaped interaction, and are omitted in the following calculations. Then, the third-order perturbation term is given by, where ∆ represents the energy gap of four visons as illustrated in Fig. 8(b). In the Majorana representation, Eq. (C17) is recast into, H This amounts to the zigzag-shaped interaction of four neighboring Majorana fermions on the sites i, j, k, l shown in Fig. 8(e). We, here, stress again that up to the third-order perturbations, the above analysis exhausts all four-body interactions among neighboring Majorana fermions which are generated by the combination of applied magnetic fields and the Heisenberg or symmetric off-diagonal exchange interactions.
We also take into account the renormalization of the nearest-neighbor and next-nearest-neighbor hopping amplitudes of itinerant Majorana fermions due to non-Kitaev interactions and magnetic fields. Then, finally, we obtain the effective Hamiltonian for itinerant Majorana fermions up to the third order perturbation, where the definitions of the indices of the hopping amplitudes, a = 1 ∼ 6, are similar to those shown in Fig. 4(b), and jk a means an a-bond connecting next-nearest-neighbor sites j and k, and, with κ ≡ 6h x h y h z /∆ 2 . In the main text, we put t 1 = t 2 = t 3 = t, and t 4 = t 5 = t 6 = κ to highlight spontaneous rotational-symmetry breaking due to four-Majorana interactions. It is worth mentioning that, as seen in Eqs.(C19)-(C21), the normalization due to the Heisenberg interaction and the Γ term with Γ > 0, reduces the nearest-neighbor hopping amplitudes in the case of the ferromagnetic Kitaev interaction, J > 0. This remarkable feature implies that for real candidate materials such as α-RuCl 3 , where the magnitudes of Γ and J are in the same order, the band width of itinerant Majorana fermions is substantially reduced, and thus, the systems may be in strongly correlated regions with g ∼ t in Eq.(5) of the main text, for which the nematic transition occurs.