Asymmetric modulation of Majorana excitation spectra and nonreciprocal thermal transport in the Kitaev spin liquid under a staggered magnetic field

The quantum spin liquid (QSL) is a spin state in insulating magnets which does not break any symmetry of the systems down to absolute zero temperature Balents ; SB ; KM . This exotic state is expected to be realized when magnetic frustration and quantum fluctuation are in action. Among many proposals, the Kitaev honeycomb model has drawn attention since its ground state provides a rare example of the QSL in more than one dimension Kitaev2006 . In the Kitaev spin liquid, the elementary excitations are described by itinerant Majorana fermions and localized $Z_{2}$ fluxes as a consequence of the fractionalization of the spin degree of freedom Kitaev2006 . The quantum entanglement between the fractional excitations would be utilized for decoherence-free topological quantum computing Kitaev2003 . Since it was pointed out that a class of materials with strong spin-orbit coupling and electron correlation is a good platform for the realization of the Kitaev model JK , extensive investigation has been performed toward the identification of the Kitaev spin liquid and the fractional excitations for the candidate materials, such as $A_{2}{\rm IrO_{4}}$ $(A={\rm Li\ and\ Na})$ and ${\alpha}$-RuCl₃ TTJKN ; MN ; Winter ; Trebst .

To capture the evidence of the fractionalization, the effect of a magnetic field has been studied both theoretically and experimentally. In theory, the magnetic field opens a gap at the Dirac-like nodal points at the Brillouin zone edges in the excitation spectrum of the itinerant Majorana fermions, and the system becomes a Majorana Chern insulator with a chiral Majorana edge mode in the weak field limit, which leads to the half-quantized thermal Hall effect Kitaev2006 ; NYM . Experimentally, the external magnetic field suppresses the parasitic long-range order stabilized by residual non-Kitaev interactions and may realize the Kitaev spin liquid in the field Wolter ; Jansa ; Banerjee2018 ; Zheng ; Nagai . Indeed, the half-quantized thermal Hall effect was observed in the field-induced quantum paramagnetic region in one of the candidates, ${\alpha}$-RuCl₃ Kasahara ; Yokoi ; Yamashita ; Bruin , which is regarded as strong evidence of the chiral Majorana edge modes in the gapped topological Majorana state.

Recently, modulations of the Majorana excitation spectrum by the magnetic field have been examined in a wider scope. For instance, an asymmetric modulation of the Majorana band structure was discussed as a combined effect of a site-dependent internal magnetic field and a non-Kitaev interaction TF2 . A similar but different type of asymmetric modulation was discussed by introducing an electric field in addition to the magnetic field CMR . In both cases, the sublattice symmetry is broken, which activates the valley degree of freedom and may bring about the Majorana Fermi surfaces depending on the parameters. Such “Majorana band engineering” would be useful for identifying the fractional excitations in the Kitaev spin liquid, but the systematic study remains yet to be investigated, especially the consequence on the thermal transport which is sensitive to the modulations of the Majorana excitation spectrum.

In this paper, we perform a systematic study of the Majorana excitation spectrum and the thermal transport in the Kitaev model, by introducing both uniform and staggered magnetic fields. Within the perturbation theory in terms of the magnetic fields, we show that the Majorana band structure is modulated in an asymmetric way depending on the magnitudes and directions of the magnetic fields, which allows us to control the valley degree of freedom. Furthermore, it becomes gapless with the formation of the Majorana Fermi surfaces in some parameter regions of the magnetic fields. For these situations, we calculate the longitudinal thermal transport coefficient by using the Boltzmann transport theory. We find that the thermal transport becomes nonreciprocal due to the asymmetric modulation of the Majorana band structure. By the comprehensive study while changing the uniform and staggered fields, we show that the linear component of the thermal current probes the Majorana excitation gap, while the nonlinear one is sensitive to the asymmetry of the Majorana band structure. We discuss our results with the estimates of the thermal currents expected in real materials. We also discuss other additional effects of the magnetic field, including the interaction between the Majorana fermions. Finally, we mention about contributions from the conventional magnons, which would be important in heterostructures with an antiferromagnet to realize an internal staggered magnetic field.

This paper is organized as follows. In Sec. II, we introduce the model that we use in this paper. In Sec. III, we discuss the Majorana excitation spectra while changing the uniform and staggered magnetic fields, with emphasis on the Majorana excitation gap, the effective Majorana density, and the asymmetry of the Majorana band structure. In Sec. IV, we present the results of the linear and nonlinear thermal transport coefficients for the modulated Majorana bands. In Sec. V, we discuss quantitative estimates of the thermal currents, the effect of the second-order perturbation and the Majorana interaction, and the contributions from the conventional magnons. Section VI is devoted to the summary.

We start from the Kitaev model on the honeycomb lattice under the magnetic field which includes both uniform and staggered components. The Hamiltonian is given by ${\cal H}={\cal H}_{\rm K}+{\cal H}_{\rm Z}$, where

	$\displaystyle{\cal H}_{\rm K}$	$\displaystyle=-\sum_{\alpha=x,y,z}\sum_{\langle i,j\rangle_{\alpha}}J^{\alpha}S_{i}^{\alpha}S_{j}^{\alpha},$		(1)
	$\displaystyle{\cal H}_{\rm Z}$	$\displaystyle=-\sum_{i}{\bm{H}}_{i}\cdot{\bm{S}}_{i}.$		(2)

Here, $S_{i}^{\alpha}={\sigma}_{i}^{\alpha}/2$ is the ${\alpha}$ component of the spin-1/2 operator at site $i$ defined by the Pauli matrix ${\sigma}_{i}^{\alpha}$ ($\alpha=x,y,z$; we set $\hbar=1$ here and for a while), $J^{\alpha}$ is the coupling constant for the Kitaev exchange interaction on the ${\alpha}$ bond, and ${\bm{H}}_{i}$ denotes the magnetic field at site $i$ given by

	$\displaystyle{\bm{H}}_{i\in{\rm A}}$	$\displaystyle={\bm{H}}-{\bm{H}}_{\rm st}=H\left({\bm{h}}-r_{\rm st}{\bm{h}}_{\rm st}\right)\equiv H{\bm{h}}_{\rm A},$		(3)
	$\displaystyle{\bm{H}}_{i\in{\rm B}}$	$\displaystyle={\bm{H}}+{\bm{H}}_{\rm st}=H\left({\bm{h}}+r_{\rm st}{\bm{h}}_{\rm st}\right)\equiv H{\bm{h}}_{\rm B},$		(4)

for site $i$ belonging to the sublattice A and B, respectively (see Fig. 1). In Eqs. (3) and (4), ${\bm{H}}$ and ${\bm{H}}_{\rm st}$ denote the uniform and staggered components, respectively; $H=|{\bm{H}}|$, $|{\bm{h}}|=|{\bm{h}}_{\rm st}|=1$, and $r_{\rm st}=|{\bm{H}}_{\rm st}|/|{\bm{H}}|$.

Following the perturbation theory with respect to the uniform magnetic field up to third order Kitaev2006 , the Zeeman coupling term ${\cal H}_{\rm Z}$ in Eq. (2) can be effectively described by three-spin interactions as

$\displaystyle{\cal H}^{\prime}$

$\displaystyle=-\delta\sum_{\{i,j,k\}}h_{i}^{x}h_{j}^{y}h_{k}^{z}S_{i}^{x}S_{j}^{y}S_{k}^{z},$

(5)

where ${\bm{h}}_{i}\equiv{\bm{H}}_{i}/H={\bm{h}}_{\rm A(B)}$ for $i\in{\rm A(B)}$, and $\delta\sim H^{3}/\Delta_{\rm f}^{2}$ with $\Delta_{\rm f}$ being a (mean) flux excitation energy. In Eq. (5), the summation is taken for neighboring three sites $\{i,j,k\}$ as exemplified in Fig. 1(a). By using the Jordan-Wigner transformation CH ; FZX ; CN , both ${\cal H}_{\rm K}$ and ${\cal H}^{\prime}$ can be written in terms of the Majorana fermions as

	$\displaystyle{\cal H}_{\rm K}$	$\displaystyle=\frac{i}{4}\sum_{{{\bm{r}}}}\Bigl{[}J^{x}c_{{{\bm{r}}},{{\rm A}}}c_{{{\bm{r}}}+{\bm{a}}_{1},{{\rm B}}}+J^{y}c_{{{\bm{r}}},{{\rm A}}}c_{{{\bm{r}}}+{\bm{a}}_{2},{{\rm B}}}+J^{z}c_{{{\bm{r}}},{{\rm A}}}c_{{{\bm{r}}},{{\rm B}}}\Bigr{]},$		(6)
	$\displaystyle{\cal H}^{\prime}$	$\displaystyle=-\frac{i\delta}{8}\sum_{{\bm{r}}}\Bigl{[}\tilde{h}_{\rm A}^{yzx}c_{{{\bm{r}}}-{\bm{a}}_{3},{{\rm A}}}c_{{{\bm{r}}},{{\rm A}}}+\tilde{h}_{\rm A}^{zxy}c_{{{\bm{r}}},{{\rm A}}}c_{{{\bm{r}}}-{\bm{a}}_{2},{{\rm A}}}+\tilde{h}_{\rm A}^{xyz}c_{{{\bm{r}}}-{\bm{a}}_{1},{{\rm A}}}c_{{{\bm{r}}},{{\rm A}}}$
		$\displaystyle\quad\quad\quad\ +\tilde{h}_{\rm B}^{yzx}c_{{{\bm{r}}},{{\rm B}}}c_{{{\bm{r}}}-{\bm{a}}_{3},{{\rm B}}}+\tilde{h}_{\rm B}^{zxy}c_{{{\bm{r}}}-{\bm{a}}_{2},{{\rm B}}}c_{{{\bm{r}}},{{\rm B}}}+\tilde{h}_{\rm B}^{xyz}c_{{{\bm{r}}},{{\rm B}}}c_{{{\bm{r}}}-{\bm{a}}_{1},{{\rm B}}}\Bigr{]},$		(7)

where $c_{{{\bm{r}}},{{\rm A}}({{\rm B}})}$ is the Majorana fermion operator at the A(B) sublattice site on the $z$ bond at the position $\bm{r}$, and $\tilde{h}_{\rm A(B)}^{{\alpha}{\beta}{\gamma}}=h_{\rm A(B)}^{\alpha}h_{\rm B(A)}^{\beta}h_{\rm A(B)}^{\gamma}$; ${\bm{a}}_{1}=(1/2,\sqrt{3}/2)$ and ${\bm{a}}_{2}=(-1/2,\sqrt{3}/2)$ are the primitive translational vectors and ${\bm{a}}_{3}\equiv{\bm{a}}_{2}-{\bm{a}}_{1}$ (we set the lattice constant as unity) [see Fig. 1(a)]. Hereafter, we consider the case with the isotropic ferromagnetic Kitaev couplings, $J^{x}=J^{y}=J^{z}=J>0$, while the results are the same for the antiferromagnetic case within the above perturbation theory. We note that the second-order perturbation with respect to the magnetic field modulates $J^{\alpha}$, and the third-order perturbation gives rise to additional three-spin terms to Eq. (5) which correspond to interactions between the Majorana fermions. We will discuss these additional effects in Sec. V.2.

In this section, we discuss the Majorana excitation spectra in the effective Hamiltonian ${\cal H}_{\rm eff}={\cal H}_{\rm K}+{\cal H}^{\prime}$. The Fourier transform of Eqs. (6) and (7) is summarized into

$\displaystyle{\cal H}_{\rm eff}=\sum_{{{\bm{k}}}\in{\rm BZ}(k_{a}>0)}{\bm{c}}_{{\bm{k}}}^{\dagger}\left(\begin{array}[]{cc}\Delta_{{{\bm{k}}},{{\rm A}}}&if_{{\bm{k}}}\\ -if_{{\bm{k}}}^{*}&\Delta_{{{\bm{k}}},{{\rm B}}}\end{array}\right){\bm{c}}_{{\bm{k}}},$

(10)

where the summation of the momentum ${\bm{k}}=(k_{a},k_{b})$ is taken over a half of the first Brillouin zone, and ${\bm{c}}_{{\bm{k}}}^{\dagger}=\left(c_{{{\bm{k}}},{{\rm A}}}^{\dagger},c_{{{\bm{k}}},{{\rm B}}}^{\dagger}\right)$ and ${\bm{c}}_{{\bm{k}}}={}^{t}\left(c_{{{\bm{k}}},{{\rm A}}},c_{{{\bm{k}}},{{\rm B}}}\right)$ are the Fourier components of the Majorana fermion operators defined as

	$\displaystyle c_{{{\bm{k}}},{{\rm A}}({{\rm B}})}^{\dagger}=\frac{1}{\sqrt{2N_{\rm u}}}\sum_{{\bm{r}}}c_{{{\bm{r}}},{{\rm A}}({{\rm B}})}{{\rm e}}^{i{{\bm{k}}}\cdot{{\bm{r}}}},$		(11)
	$\displaystyle\ \ c_{{{\bm{k}}},{{\rm A}}({{\rm B}})}=\frac{1}{\sqrt{2N_{\rm u}}}\sum_{{\bm{r}}}c_{{{\bm{r}}},{{\rm A}}({{\rm B}})}{{\rm e}}^{-i{{\bm{k}}}\cdot{{\bm{r}}}},$		(12)

where $N_{\rm u}$ is the number of the unit cell, and

	$\displaystyle f_{{\bm{k}}}$	$\displaystyle=\frac{J}{2}\left({{\rm e}}^{i{{\bm{k}}}\cdot{\bm{a}}_{1}}+{{\rm e}}^{i{{\bm{k}}}\cdot{\bm{a}}_{2}}+1\right),$		(13)
	$\displaystyle\Delta_{{{\bm{k}}},{{\rm A}}}$	$\displaystyle=\frac{\delta}{2}\left(\tilde{h}_{\rm A}^{xyz}\sin{{\bm{k}}}\cdot{\bm{a}}_{1}-\tilde{h}_{\rm A}^{zxy}\sin{{\bm{k}}}\cdot{\bm{a}}_{2}+\tilde{h}_{\rm A}^{yzx}\sin{{\bm{k}}}\cdot{\bm{a}}_{3}\right),$		(14)
	$\displaystyle\Delta_{{{\bm{k}}},{\rm B}}$	$\displaystyle=-\frac{\delta}{2}\left(\tilde{h}_{\rm B}^{xyz}\sin{{\bm{k}}}\cdot{\bm{a}}_{1}-\tilde{h}_{\rm B}^{zxy}\sin{{\bm{k}}}\cdot{\bm{a}}_{2}+\tilde{h}_{\rm B}^{yzx}\sin{{\bm{k}}}\cdot{\bm{a}}_{3}\right).$		(15)

Note that the Majorana operators Eqs. (11) and (12) satisfy $c_{{{\bm{k}}},{{\rm A}}({{\rm B}})}^{\dagger}=c_{-{{\bm{k}}},{{\rm A}}({{\rm B}})}$. Then, the eigenenergy is given by the diagonalization of Eq. (10) as

$\displaystyle E_{{{\bm{k}}}\pm}=\frac{1}{2}\left(\Delta_{{{\bm{k}}}+}\pm\sqrt{4|f_{{\bm{k}}}|^{2}+\Delta_{{{\bm{k}}}-}^{2}}\ \right),$

(16)

with $\Delta_{{{\bm{k}}}\pm}\equiv\Delta_{{{\bm{k}}},{{\rm A}}}\pm\Delta_{{{\bm{k}}},{{\rm B}}}$. Note that $E_{{{\bm{k}}}\pm}$ is always particle-hole symmetric with respect to zero energy, namely $E_{{{\bm{k}}}\pm}=-E_{-{{\bm{k}}}\mp}$, because of the nature of the Majorana fermions, and the only positive energy parts contribute to the physical properties.

Let us first discuss the overall behavior of the dispersion relation in Eq. (16) while changing the uniform and staggered components of the magnetic fields. In the absence of the magnetic field $H=0$, $E_{{{\bm{k}}}\pm}$ has the Dirac-like linear dispersions with the gapless nodal points at the K and K’ points, as shown in Fig. 2(a) Kitaev2006 . When introducing a uniform magnetic field, the nodal points are gapped out, as shown in Fig. 2(b). In this situation, the insulating state becomes topologically nontrivial having a nonzero Chern number Kitaev2006 . Meanwhile, when we apply a staggered magnetic field only, the nodal points remain gapless, but their energies at the K and K’ points are shifted in an opposite energy direction, as shown in Fig. 2(c). This makes the two valleys at the K and K^′ points inequivalent and $E_{{{\bm{k}}}\pm}$ asymmetric in momentum space, because of the breaking of sublattice symmetry. As a consequence, the staggered magnetic field yields the Majorana Fermi surfaces around the K and K’ points, whose sizes are increased with increasing the strength of the staggered magnetic field, as shown in the right panel of Fig. 2(c). Finally, when we apply both uniform and staggered fields, the gap opening and the asymmetric energy shift occur simultaneously, as shown in Fig. 2(d). In this case, the system can be insulating or gapless with the Majorana Fermi surfaces depending on the magnitudes of the uniform and staggered components, as demonstrated in the right panel of Fig. 2(d).

Figure 3 displays the Majorana gap $\Delta_{\rm M}$ while varying the direction of the uniform component of the magnetic field, ${\bm{h}}$, with a staggered component ${\bm{h}}_{\rm st}$ along the $a$ direction. When the staggered field is zero or small, the Majorana excitation spectrum is gapped in most cases, as shown in Figs. 3(a) and 3(b). In the gapped region, the system is the Majorana Chern insulator as in the case with the uniform magnetic field only. With an increase of the staggered field, the asymmetric energy shift at the K and K’ points becomes larger, leading to an increase of the gapless region, as shown in Figs. 3(c) and 3(d).

Figure 4 shows the asymmetric modulation of $E_{{{\bm{k}}}\pm}$ by plotting the constant energy cuts at $E_{{{\bm{k}}}\pm}=0.35J$ around the K and K’ points. In Fig. 4(a), we present the results while rotating ${\bm{h}}$ within the $ac$ plane with ${\bm{h}}_{\rm st}$ along the $a$ direction. When ${\bm{h}}$ is parallel to ${\bm{h}}_{\rm st}$ ($\theta=\pi/2$), the effective magnetic field $\tilde{h}_{A(B)}^{\alpha\beta\gamma}$ in Eq. (7) becomes zero, and hence, the dispersion is equivalent to that in the absence of the magnetic field and the constant energy cuts are symmetric between the K and K’ points. The rotation of the direction of ${\bm{h}}$ modulates the dispersion in an asymmetric way; when the constant energy cut around the K’ point is enlarged, that around the K point is shrunk. We note that the asymmetry is always present in the $k_{a}$ direction except for $\theta=\pi/2$, while the dispersion becomes symmetric in the $k_{b}$ direction for $\theta=0$, $\pi/4$, and $3\pi/4$. On the other hand, when ${\bm{h}}$ is rotated within the $ab$ plane, the dispersion becomes asymmetric in both $k_{a}$ and $k_{b}$ directions, except for the symmetric case with $\varphi=0$, as shown in Fig. 4(b). Finally, Fig. 4(c) shows the results while rotating ${\bm{h}}$ within the $ab$ plane with ${\bm{h}}_{\rm st}$ along the $c$ direction $(\theta=0)$, which demonstrates that the asymmetry is changed by the rotation of ${\bm{h}}_{\rm st}$.

In this section, we calculate the longitudinal thermal conductivity in the presence of both uniform and staggered magnetic fields, by using the Boltzmann transport theory. We apply a weak thermal gradient to the $a$ or $b$ direction. Then, the thermal current can be expanded in terms of the thermal gradient as

$\displaystyle j_{Q,\xi}=\kappa_{\xi}^{(1)}\left(-\frac{\partial T}{\partial\xi}\right)+\kappa_{\xi}^{(2)}\left(\frac{\partial T}{\partial\xi}\right)^{2}+\cdots$

(17)

where $\xi=a$ or $b$; $\kappa_{\xi}^{(1)}$ and $\kappa_{\xi}^{(2)}$ are the linear and nonlinear thermal conductivities given by

	$\displaystyle\kappa_{\xi}^{(1)}$	$\displaystyle=\frac{2\tau}{\Omega}\sum_{n}\sum_{{{\bm{k}}}\in{\rm BZ}(E_{{{\bm{k}}}n}>0)}\frac{\partial f(E_{{{\bm{k}}}n})}{\partial T}E_{{{\bm{k}}}n}v_{{{\bm{k}}}n\xi}^{2},$		(18)
	$\displaystyle\kappa_{\xi}^{(2)}$	$\displaystyle=\frac{2\tau^{2}}{\Omega}\sum_{n}\sum_{{{\bm{k}}}\in{\rm BZ}(E_{{{\bm{k}}}n}>0)}\frac{\partial^{2}f(E_{{{\bm{k}}}n})}{\partial T^{2}}E_{{{\bm{k}}}n}v_{{{\bm{k}}}n\xi}^{3},$		(19)

respectively, where the summation is taken over the momentum ${\bm{k}}$ in the first Brillouin zone for which $E_{{\bm{k}}n}>0$, $\Omega$ is the system size, $n=\pm$ is the band index, $f(E_{{{\bm{k}}}n})=({{\rm e}}^{E_{{{\bm{k}}}n}/T}+1)^{-1}$ is the Fermi distribution function, $v_{{{\bm{k}}}n\xi}\equiv\partial E_{{{\bm{k}}}n}/\partial k_{\xi}$ is the group velocity, and $\tau$ is the relaxation time of Majorana fermion introduced phenomenologically. We set the Boltzmann constant $k_{\rm B}$ unity and take $J\tau=1$ in this section.

The thermal current is generated by thermally excited Majorana fermions. To discuss the field dependence of the linear response $\kappa_{\xi}^{(1)}$, we introduce the “effective Majorana density” as

$\displaystyle D_{\rm M}=\int_{0}^{\infty}\sum_{n}\sum_{{{\bm{k}}}\in{\rm BZ}(E_{{{\bm{k}}}n}>0)}\delta(\varepsilon-E_{{{\bm{k}}}n})\left(-\frac{df(\varepsilon)}{d\varepsilon}\right)d\varepsilon,$

(20)

in addition to the Majorana gap $\Delta_{\rm M}$; see Sec. IV.1. Meanwhile, $\kappa_{\xi}^{(2)}$ describes the second-order nonlinear response arising from the asymmetry of the Majorana band dispersions. To measure the asymmetry, we define the “band asymmetry parameter” as

$\displaystyle{\Lambda}_{\xi}=-\frac{1}{\Omega}\sum_{{{\bm{k}}}\in{\rm BZ}(k_{\xi}>0)}\sum_{n}\frac{\partial f_{{{\bm{k}}}n}}{\partial E_{{{\bm{k}}}n}}(v_{{{\bm{k}}}n\xi}+v_{-{{\bm{k}}}n\xi}),$

(21)

where the sum of momentum is taken in half of the first Brillouin zone with $k_{\xi}>0$ to extract the asymmetry. We use this quantity to discuss the field dependences of $\kappa_{\xi}^{(2)}$ in Sec. IV.2. We perform the following calculations for the system with periodic boundary conditions by taking $864\times 864$ $k$ points in the first Brillouin zone.

To summarize, we have studied the thermal transport in the Kitaev model under uniform and staggered magnetic fields by using the Boltzmann transport theory. The staggered field breaks the sublattice symmetry and activates the valley degree of freedom. Depending on the amplitudes and directions of the uniform and staggered magnetic fields, the band structure of the Majorana fermions, which are yielded by the fractionalization of spins in the quantum spin liquid state in the Kitaev model, is asymmetrically modulated in momentum space, hosting the gapless or gapped nodal points at finite energy and the Majorana Fermi surfaces. We found that the linear thermal current shows a characteristic field dependence, which correlates with the magnitude of the Majorana gap or the effective Majorana density in the presence of the Majorana Fermi surfaces. On the other hand, we showed that the asymmetric band modulation gives rise to the nonlinear thermal transport, which also depends on the field amplitude and directions. We presented the quantitative estimates by referring to the material parameters of ${\alpha}$-RuCl₃, and found that the linear thermal response is observable while the improvement of the experimental accuracy is required for the nonlinear component. We also discussed how other additional effects of the magnetic fields arising from the second- and third-order perturbation modify the behavior of the thermal transport. Furthermore, we discussed contributions from magnons in a heterostructure to realize the internal staggered magnetic field, which would be distinguished from the Majorana contributions.

Our findings illuminate the Majorana band modulations by the uniform and staggered magnetic fields and their observation through the thermal transport measurements, but there remain several open questions. The present analysis has been done by the perturbation theory in terms of the magnetic fields, in which the quantum spin liquid state without the localized $Z_{2}$ fluxes (flux-free state) is retained. The flux-free state, however, does not hold beyond the perturbative regime. The $Z_{2}$ fluxes are also excited by raising temperature, which was shown to affect the thermal transport NYM . Another interesting unperturbed effect is a field-induced gapless spin liquid state in the antiferromagnetic Kitaev model ZKSF ; LJCLW ; GMP ; NKKM ; HT ; RVT ; JDJ , which would also affect the thermal transport significantly TZSSS . In addition, in the candidate materials like $\alpha$-RuCl₃, the effect of non-Kitaev interactions, such as the Heisenberg and off-diagonal interactions, cannot be neglected. Despite some recent attempts YNKM ; MKM , it is still a challenge to access the interesting parameter region at low temperature in applied magnetic fields. Further development is required to clarify these important issues.

We acknowledge S. Okumura, R. Sano, and A. Tsukazaki for fruitful discussions. This work is supported by JST CREST (JP-MJCR18T2).