Topological Nematic Phase Transition
in Kitaev Magnets Under Applied Magnetic Fields

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 Kitaev (2006) 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 Kitaev (2006, 2003). After the Kitaev’s seminal paper, various phenomena associated with Majorana fermions in Kitaev magnets have been extensively explored Baskaran et al. (2007); Feng et al. (2007); Knolle et al. (2014); Nasu et al. (2015); Hermanns et al. (2015a, b); O’Brien et al. (2016); Song et al. (2016); Nasu et al. (2017a); Udagawa (2018); Yamada et al. (2017a); Yamada (2021); Takikawa and Fujimoto (2019, 2020). Since Jackeli and Khaliullin Jackeli and Khaliullin (2009); Yamaji et al. (2014); Rau et al. (2014); Kim and Kee (2016); Winter et al. (2016); Yamada et al. (2017b) 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 $\mathrm{Na_{2}IrO_{3}}$ Lee et al. (2010); Singh and Gegenwart (2010), $\alpha$-$\mathrm{Li_{2}IrO_{3}}$ Singh et al. (2012), and $\alpha$-RuCl₃ Plumb et al. (2014). $\alpha$-RuCl₃ under an applied magnetic field is one of the best candidates for the Kitaev magnet Baek et al. (2017). 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 Kasahara et al. (2018); Vinkler-Aviv and Rosch (2018); Ye et al. (2018). Moreover, according to the recent specific heat measurements Tanaka et al. (2020), 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 ¹¹1Although the crystal structure of $\alpha$-RuCl₃ may be monoclinic rather than trigonal Johnson et al. (2015); Kim and Kee (2016), the field angle dependence of heat capacity observed in $\alpha$-RuCl₃ by Tanaka et al. Tanaka et al. (2020) shows the $C_{3}$ symmetry within error bars for magnetic fields $<9$ T. This approximated threefold rotational symmetry drastically breaks down to obvious twofold one for strong magnetic fields $>9$ T, which implies the possibility of nematic phase transition in $\alpha$-RuCl₃ Tanaka et al. (2020). Remarkably, the nematic phase transition accompanies the disappearance of the half-quantized thermal Hall effect Kasahara et al. (2018), 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, $\Gamma$, and $\Gamma^{\prime}$ terms.

In this work, we theoretically discuss the possibility of a nematic phase transition due to many-body interactions among itinerant Majorana fermions in Kitaev magnets. Effects of interactions between Majorana fermions have been extensively studied for Kitaev magnets as well as topological superconductors Stoudenmire et al. (2011); Niu et al. (2012); Thomale et al. (2013); Sticlet et al. (2014); Hermanns et al. (2015b); Rahmani et al. (2015); Affleck et al. (2017); Li and Franz (2018); Wamer and Affleck (2018); Rahmani et al. (2019); Rahmani and Franz (2019); Sannomiya and Katsura (2019); Li et al. (2019); Roy et al. (2020); Tummuru et al. (2021); Vishveshwara and Weld (2020); Sachdev and Ye (1993); Maldacena and Stanford (2016); Rahmani and Franz (2019); Yamada (2020). However, physics arising from Majorana interactions in real Kitaev materials is still elusive. On the other hand, various numerical studies on extended Kitaev model utilizing the spin representation, which effectively include non-perturbative effects of Majorana interactions, have been carried out Chaloupka et al. (2010); Okamoto (2013); Rau et al. (2014); Chaloupka and Khaliullin (2016); Nasu et al. (2017b); Gohlke et al. (2018); Rusnačko et al. (2019); Gordon et al. (2019); Jiang et al. (2019); Kaib et al. (2019); Lee et al. (2020); Gohlke et al. (2020); Liu et al. (2021); Buessen and Kim (2021), and nematicity is discussed in recent studies Nasu et al. (2017b); Gohlke et al. (2018); Lee et al. (2020); Gohlke et al. (2020); Liu et al. (2021). Rich phase diagrams of Kitaev magnets are obtained from these studies based on the spin Hamiltonian, while a further theoretical study is needed to understand the results of the recent specific heat measurements mentioned above, which we believe is from the gapped QSL phase to another QSL phase Tanaka et al. (2020).

To focus on the nature of a QSL phase in Kitaev magnets from a different point of view, we assume the flux-free state in our model. Although the flux-free assumption may be violated by strong additional spin-exchanges like Heisenberg, $\Gamma$, and $\Gamma^{\prime}$ 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. Gohlke et al. (2018); Knolle et al. (2018), 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, 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.

The Hamiltonian of the pure Kitaev model is expressed as,

where $\sigma_{i}$ represents the spin operator on the $i$th site and $\langle ij\rangle$ 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 $\pm 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 vortex-free 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 $\bm{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 “$\rightarrow$” means taking the standard gauge in the 0-flux ground state, and $\langle\!\langle ij\rangle\!\rangle$ 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 $\pm 1$, and the others are four-Majorana interactions. Here the site summation for $\sum_{\mathrm{Y}}$ is taken as shown in the upper Y-shaped diagram of Fig. 1(c) and $\sum_{\mathrm{Y^{\prime}}}$ 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, $H_{\Gamma^{\prime}}=\Gamma^{\prime}\sum_{\langle ij\rangle_{\alpha},\beta\neq\alpha}\left[\sigma_{i}^{\alpha}\sigma_{j}^{\beta}+\sigma_{i}^{\beta}\sigma_{j}^{\alpha}\right]$ where $\alpha,\,\beta$ 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}{\Gamma^{\prime}}^{2})$. Besides, other non-Kitaev interactions, such as Heisenberg and $\Gamma$ 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 $\kappa$ and $g$ are independent for real candidate materials of Kitaev magnets, as mentioned above. The Majorana operator is represented as $c_{i}$ for the $i$th site, and the operators obey $c_{i}^{\dagger}=c_{i}$ and $\{c_{i},c_{j}\}=2\delta_{ij}$. The hopping amplitude $t=J$ is set unity without loss of generality.

First, we apply the mean-field (MF) analysis to the four-body interaction term $\mathcal{H}_{4}$. Then, the general form of the MF Hamiltonian is,

where the index $a=1,\dots,6$ specifies the directions of hopping with the amplitude $\tau_{a}$, as described in Fig. 1(d). The phase factors $\eta_{ij}=\pm 1$ is necessary to make $\mathcal{H}_{\mathrm{MF}}$ antisymmetric. The sign of $\eta_{ij}$ is defined as illustrated in Fig. 1(d) such that the direction of the arrow indicates the positive hopping from the $j$th site to the $i$th 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 $\ket{\Psi_{\mathrm{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 $\Delta_{a}$ ($a=1,\dots,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 Pujari et al. (2015),

If $\Delta_{1}=\Delta_{2}=\Delta_{3}$ ($\Delta_{4}=\Delta_{5}=\Delta_{6}$), $\phi$ ($\psi$) must be 0 from the definition, but otherwise will have a nonzero value. In other words, $\phi$ and $\psi$ characterizes breaking of the threefold rotational symmetry.

As seen in Figs. 2(a) and (c), for large value of $|g|$, a nematic phase with $|\phi|\neq 0$, $|\psi|\neq 0$, characterized by two-fold rotational symmetry, appears. We also find in Fig. 2(b) that $\arg(\phi)$ is $\pi/3$ in most regions of the $|\phi|\neq 0$ area. The argument of $\phi$ has informations about the direction of the nematic order. Since $\Delta_{a}$ takes a negative value in this calculation, $\arg(\phi)=\pi/3$ corresponds to the nematic order with $|\Delta_{3}|>|\Delta_{1}|,|\Delta_{2}|$. On the other hand, as seen in Fig. 2(d), $\arg(\psi)=\pi/3$ for $g<0$ in most of the nematic region, while $\arg(\psi)=-2\pi/3$ for $g>0$. This is because that $\psi\rightarrow-\psi$ under time-reversal operation, while $\phi\rightarrow\phi$. In the calculations, we intentionally introduce tiny anisotropy of bonds to obtain symmetry-breaking 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 Fukui et al. (2005). The results are shown in Fig. 2(e). In the chiral QSL phase without the nematic order, the Chern number $\nu=\pm 1$. Remarkably, on the other hand, in the nematic phase with $|\phi|\neq 0$ and $|\psi|\neq 0$, we find $\nu=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 ($\tau_{1}$, $\tau_{2}$, $\tau_{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 Kitaev (2006), the gapful Majorana nematic phase with the Chern number $\nu=0$ is identified with the toric code phase with Abelian anyons. It is noted that this nematic phase transition is the first-order type with the discontinuous changes of the order parameters $\phi$ and $\psi$, and hence, there is no gap closing at the transition point. We stress that although the nonzero $\psi$ apparently causes the anisotropy of the Majorana hopping terms, both of the two nematic orders $\phi$ and $\psi$ cooperatively stabilize the toric code phase, since $\phi$ and $\psi$ are nonlinearly coupled with each other via the MF equations. The details are shown in Appendix B.

At the end of this section, we should mention another phase realized in small regions with $\nu=\pm 1$ inside nematic phases shown in Fig.2(e). Although this phase is also a nematic phase with $|\phi|\neq 0$, $\arg(\phi)=-2\pi/3$ and hence, $|\Delta_{1}|,|\Delta_{2}|>|\Delta_{3}|$ is realized. We call this phase “zigzag nematic phase”. This zigzag nematic phase is in $B$ phase, i.e. with the Chern number $\nu=\pm 1$.

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), $\phi$ and $\psi$ 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 $\{\phi,\phi^{\dagger}\}/2$ and $\{\psi,\psi^{\dagger}\}/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 $|\phi|$ and $|\psi|$, and the arguments of $\phi$ and $\psi$. 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(\phi)=\pi+2n\pi/3\,(n\in\mathbb{Z})$, and the other is a zigzag nematic phase with two strong bonds so that $\arg(\phi)=2n\pi/3\,(n\in\mathbb{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 number $\nu=\pm 1$, where $\arg(\phi)=-2\pi/3$, $\arg(\psi)=-2\pi/3$ ($\arg(\phi)=-2\pi/3$, $\arg(\psi)=\pi/3$) for $g<0$ ($g>0$), 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.

Using the MF analysis and the exact diagonalization method, we find that the topological nematic phase transition occurs in strongly correlated regions $|g|\sim t$, $g\sim\kappa$. Since non-Kitaev interactions such as off-diagonal exchange interactions significantly renormalize $t$, $\kappa$, and $g$, it is possible to realize this condition for the extended Kitaev model. Particularly, the $\Gamma$ term reduces the hopping amplitude $t$, and drives the system into strongly correlated regions. We note that the Y(Y’)-shaped interactions of order $\sim O(h_{x,y,z}{\Gamma^{\prime}}^{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 off-diagonal exchange interaction, i.e. the $\Gamma$ 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 $\alpha$-RuCl₃ under parallel magnetic fields Tanaka et al. (2020) 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 Zhang et al. (2019, 2020). From an experimental side, an electric probe through the hyperfine interaction Yamada and Fujimoto (2020) may be useful for the detection of the Majorana nematic phase.

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.