Transport in honeycomb lattice with random $\pi$-fluxes: Implications for low-temperature thermal transport in the Kitaev spin liquids

Quantum spin liquids (QSLs) are an exotic state of matter with no local symmetry breaking, deconfined fractionalized quasiparticles, and emergent gauge excitations [1, 2, 3]. The Kitaev model, an exactly solvable honeycomb model, provides a typical framework for describing such a phase, where spins fractionalize to Majorana fermions and $Z_{2}$ gauge excitations (visons) [4]. Thanks to the proposal by Jackeli and Khaliullin [5], which suggests that the Kitaev model can be realized in certain strongly spin-orbit coupled systems, a growing number of Kitaev candidate materials have been discovered in the past decade [6, 7].

Despite the theoretical proposal of QSLs for several decades, their experimental identification remains challenging. Thermal transport experiments are a promising technique for characterizing QSLs, as they allow the detection of charge-neutral, mobile quasiparticles. For example, in the Kitaev candidate material $\alpha-\text{RuCl}_{3}$, a half-quantized thermal Hall conductance has been observed, indicating the existence of a chiral Majorana mode at the edge of the non-Abelian Kitaev spin liquid [8, 9, 10]. To fully understand the experimental results, it is crucial to predict the thermal transport signatures of QSLs, particularly those proximate to the exact Kitaev model, in various different regimes [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21].

The longitudinal thermal conductivity $\kappa_{\text{xx}}(\omega,T)$ of the Kitaev model has been numerically investigated using Kubo’s formalism, and the DC thermal conductivity $\kappa_{\text{K}}(T)$ can be obtained by extrapolation of $\kappa_{\text{xx}}(\omega,T)$ to the $\omega\rightarrow 0$ limit [11, 12, 13, 14, 15]. When the temperature $T$ is comparable to or smaller than the vison gap $\Delta_{v}$, such extrapolation appears to give results with significant error bars and calculations with larger system sizes are needed [11, 15], leaving the low-temperature behavior of $\kappa_{\text{K}}$ inconclusive. Neglecting the gauge excitations in the Kitaev model, it is possible to regard an undoped graphene as two stacks of the Kitaev model, and thus expect that the thermal transport of Kitaev spin liquid and the electric transport of graphene may have similar low-temperature behavior. It was proposed that in graphene there exists a universal minimum electrical conductivity $\sigma_{\text{min}}=e^{2}/\pi h$ per valley per spin, regardless of the concentration of disorders [22, 23, 24]. It is thus natural to ask whether a similar ‘minimum thermal conductivity’ (MTC) also exists in a pure Kitaev model. A previous quantum Monte Carlo study seems to support the existence of MTC in the low-temperature regime, i.e. $\lim_{T\rightarrow 0}\kappa_{\text{K}}/T=k_{B}^{2}/12\hbar$, but further justification is required due to the significant uncertainties and finite-size effects in the numerical simulation [11]. Moreover, at low temperatures, the phase coherence length may exceed the dimensions of a realistic sample, within which the quantum interference effects are essential. Consequently, the thermally excited random $Z_{2}$ flux background may strongly localize Majorana particles [25, 14, 26], provided that the localization length is smaller than the system size. Such localization effects are uncaptured by previous calculations performed for relatively small systems, and overlooking these effects may crucially affect the interpretation of experimental results.

To understand the low-temperature thermal conductivity of the Kitaev model in the absence of disorder and magnetic fields, in this work we numerically study the transport of complex fermions in a random $\pi$-flux honeycomb model (RPFHM) near half-filling. This approach is motivated by the observation that the Kitaev model (see Fig. 1) can be seen as half of the RPFHM (see Fig. 1), and the thermal conductivity of the former is linked to the electrical conductivity of the latter via the Wiedemann-Franz (WF) law. Furthermore, the RPFHM is of academic interest on its own, as this random $Z_{2}$-flux problem is less studied compared to random $U(1)$-flux models [27, 28, 29, 30, 31].

We focus on the zero-temperature longitudinal DC conductivity of the RPFHM and its localization properties as a function of flux density and chemical potential. We extract the dependence of DC conductivity on flux density and Fermi wavevector in the semi-classical diffusive regime. Our analysis reveals Anderson localization due to random $\pi$-fluxes, and we obtain 2-D localization lengths. Our results indicate that the low-temperature thermal conductivity of a clean Kitaev model neither exhibits MTC nor vanishes, but instead diverges as $\kappa_{\text{K}}\sim T^{3}e^{\Delta_{v}/k_{B}T}$, which is unexpected based on previous studies. Additionally, when $k_{B}T\sim\Delta_{v}$, the thermal conductivity could be significantly suppressed due to strong Majorana localization.

The low-temperature behavior of thermal conductivity $\kappa_{\text{K}}$ of the Kitaev model can be inferred using earlier results. We first briefly introduce the Kitaev model [4]

where Pauli matrices $\vec{\sigma_{i}}=(\sigma^{x}_{i},\sigma^{y}_{i},\sigma^{z}_{i})$ describe the spin degrees of freedom on site $i$, $\langle ij\rangle_{\lambda}$ denotes the $\lambda$-type nearest-neighbor bond ($\lambda=x,y,z$) between site $i$ and $j$ and each bond is only summed once, as shown in Fig. 1. In this work, we are only interested in the isotropic case so we set $J=J_{x}=J_{y}=J_{z}$ hereafter. The Kitaev model can be exactly solved by mapping spins to Majorana fermions and the Hamiltonian becomes

where the Majorana fermions satisfy anticommutation relations $\{\gamma_{i},\gamma_{j}\}=2\delta_{ij}$, and the bond operator $\tilde{u}_{ij}=-\tilde{u}_{ji}=\pm 1$ acts as a $Z_{2}$ gauge field.

After making a gauge transformation $f_{j}^{\dagger}\rightarrow-if_{j}^{\dagger},j\in A$ sublattice, followed by a Majorana transformation $f_{i}^{\dagger}=(\gamma_{i}^{1}+i\gamma_{i}^{2})/2$ to Eq. (1), one can see the RPFHM (1) is exactly two copies of the Kitaev model (10), with $t=2J$. Therefore the thermal conductivity of the Kitaev model (9) is half of that of the RPFHM (1), namely $\kappa_{\text{G}}=2\kappa_{\text{K}}$, provided both have the same flux configuration. At thermal equilibrium, $\kappa_{\text{K}}$ can be expressed as the weighted average of thermal conductivity $\kappa_{\text{K}}^{\vbox{\hbox{\scalebox{0.7}{$\varhexagon$}}}}$ over every flux configuration $\varhexagon$, namely $\kappa_{\text{K}}=(\sum_{\vbox{\hbox{\scalebox{0.7}{$\varhexagon$}}}}Z_{\vbox{\hbox{\scalebox{0.7}{$\varhexagon$}}}}\kappa_{\text{K}}^{\vbox{\hbox{\scalebox{0.7}{$\varhexagon$}}}})/(\sum_{\vbox{\hbox{\scalebox{0.7}{$\varhexagon$}}}}Z_{\vbox{\hbox{\scalebox{0.7}{$\varhexagon$}}}})$. Here $Z_{\vbox{\hbox{\scalebox{0.7}{$\varhexagon$}}}}=\text{Tr}e^{-\beta H_{\vbox{\hbox{\scalebox{0.7}{$\varhexagon$}}}}}$ is the partition function for specific flux configuration $\varhexagon$ described by the Hamiltonian $H_{\vbox{\hbox{\scalebox{0.7}{$\varhexagon$}}}}$, and $\sum_{\vbox{\hbox{\scalebox{0.7}{$\varhexagon$}}}}$ denotes sum over all flux configurations. In the Kitaev model the ground state contains no flux [62], while at temperatures much lower than the vison gap $\Delta_{v}\approx 0.1536J$ [4], thermally excited visons are dilute and may be regarded as uncorrelated on different plaquettes, with flux density approximately given by

Here we emphasize that although the calculation of $\kappa_{\text{K}}$ seems a quenched disorder problem, it is an artifact as the bond variables $u_{ij}$ commute with the Hamiltonian. In fact, it is still an annealed average problem, as the free energy is given by $\mathcal{F}=-k_{B}T\ln\sum_{\vbox{\hbox{\scalebox{0.7}{$\varhexagon$}}}}Z_{\vbox{\hbox{\scalebox{0.7}{$\varhexagon$}}}}$, where the logarithm is taken after the sum over all flux configurations is computed.

The thermal conductivity $\kappa_{\text{G}}$ of RPFHM is linked to the electrical conductivity $\sigma_{G}$ by the Wiedemann–Franz (WF) law $\kappa_{\text{G}}=\sigma_{\text{G}}\mathcal{L}T$ where $\mathcal{L}$ is the Lorentz number. For a non-interacting Fermi liquid, the Lorentz number is a universal constant $\mathcal{L}_{0}=\pi^{2}k_{B}^{2}/3e^{2}$, while around the Dirac point, the Lorentz number may be a few times larger than $\mathcal{L}_{0}$ even for the non-interacting case [63, 64]. Combining the above arguments, one obtains

where $\sigma_{\text{G}}(n_{v},T,\mu)$ denotes the DC conductivity of RPFHM, which has flux density $n_{v}$, Fermi level $\mu$ and temperature $T$.

We shall now demonstrate the low-temperature behavior of $\kappa_{\text{K}}$ using Eq. (12). As $k_{B}T\gg|\mu|$, the conduction is mainly contributed by the thermally excited quasiparticles with energy $E\sim k_{B}T$. When $k_{B}T\ll\Delta_{v}$, we observe that the thermal de Broglie wavelength $\lambda_{\text{th}}\sim T^{-1}$ is much shorter than the mean flux distance $r_{m}\sim e^{\Delta_{v}/2k_{B}T}$. This indicates that the transport is away from the Dirac point and in the semiclassical regime. As the localization length is exponentially large $\xi\approx l_{0}e^{\pi k_{F}l_{e}}$ in this regime, the localization effects are negligible for realistic system size, or when the phase coherence length is much smaller than the localization length $L_{\phi}\ll\xi$ due to the coupling to environments. Substituting the thermal de Broglie wavevector $k_{\text{th}}=2\pi/\lambda_{\text{th}}$ to Eq. (4) (see appendix A for the justification), which is applicable when $k_{\text{th}}r_{m}\gg 1$, and using Eqs. (11)(12), one obtains diverging thermal conductivity $\kappa_{\text{K}}\sim T^{3}e^{\Delta_{v}/k_{B}T}$ at low temperatures. When $k_{B}T\sim\Delta_{v}\approx 0.15J$, which corresponds to $\lambda_{\text{th}}\gtrsim r_{m}$, the thermal transport is in the quantum regime and localization effects are non-negligible. We assume both the system size $L$ and phase coherence length $L_{\phi}$ is much larger than the thermal localization length $\xi_{T}=\xi(n_{v}(T),k_{\text{th}}(T))$, which is estimated to be a few hundred of lattice spacing from Fig. 5. In this case, the Majorana fermions are strongly localized and significantly suppressed thermal conductivity $\kappa_{\text{K}}/T\ll k_{B}^{2}/2\pi\hbar$ may be observable. This is in contrast with the low-temperature case $k_{B}T\ll\Delta_{v}$ where $\kappa_{\text{K}}/T\gg k_{B}^{2}/2\pi\hbar$ and analogous to the large resistivity $\rho\gg h/e^{2}$ observed in graphene due to the strong localization [65, 66, 67]. We note that it is difficult to obtain $\kappa_{\text{K}}$ at higher temperature $k_{B}T\sim J\gg\Delta_{v}$ using our calculations, as the thermal transport is also contributed by quasiparticles at high energies that are not considered in this work. Besides, in this regime visons can no longer be regarded as dilute, and the validity of Eq. (11) becomes questionable.

In this work, we have utilized a combination of numerical methods to investigate the transport of RPFHM near the Dirac point. Our results reveal that when the wavelength of Dirac fermion is much shorter than the average flux spacing, the semiclassical DC conductivity is quadratic in the Fermi momentum and inversely proportional to the flux density. We have also demonstrated that intervalley scattering by the $\pi$-fluxes leads to weak (strong) localization of Dirac fermions at high (low) Fermi energy.

Our results imply that the thermal transport of the Kitaev model is semiclassical at low temperatures $k_{B}T\ll\Delta_{v}$, and that the thermal conductivity diverges as $\kappa_{\text{K}}\sim T^{3}e^{\Delta_{v}/k_{B}T}$. When $k_{B}T\sim\Delta_{v}$, the itinerant Majorana fermions may be strongly localized, leading to a significantly suppressed thermal conductivity $\kappa_{\text{K}}/T\ll k_{B}^{2}/2\pi\hbar$. Our predictions are expected to hold even in the presence of disorder as long as $\tau_{v}^{-1}\gg\tau_{d}^{-1}$, where $\tau_{v}^{-1}(\tau_{d}^{-1})$ is the scattering rate due to vison (disorder). When the temperature is much smaller than the crossover temperature $T_{d}$ at which $\tau_{d}^{-1}=\tau_{v}^{-1}$, the disorder becomes the dominant source of elastic scattering and the transport is similar to that of an undoped graphene. Consequently, $\kappa_{\text{K}}/T$ should reduce with decreasing temperature and finally saturates to a value comparable or smaller to $1/2\pi\hbar$ [35, 52, 53]. Such temperature dependence, as shown in Fig. 6, was unseen in previous numerical simulations and may be observable in a clean sample at very low temperatures. In realistic Kitaev materials, non-Kitaev-type interactions also allow the hopping of visons and may have important consequences [16, 18], which are left for future studies.