# Magnetic-field Induced Topological Transitions and Thermal Conductivity
in a Generalized Kitaev Model

Heqiu Li  Department of Physics, University of Toronto, Toronto, Ontario M5S
1A7, Canada  Yong Baek Kim<EMAIL_ADDRESS>Department of Physics,
University of Toronto, Toronto, Ontario M5S 1A7, Canada  School of Physics,
Korea Institute for Advanced Study, Seoul 02455, Korea  Hae-Young Kee
<EMAIL_ADDRESS>Department of Physics, University of Toronto,
Toronto, Ontario M5S 1A7, Canada  Canadian Institute for Advanced Research,
CIFAR Program in Quantum Materials, Toronto, Ontario M5G 1M1, Canada

###### Abstract

Recent experiments on Kitaev spin liquid candidate materials reported non-
monotonic behavior of thermal conductivity as a function of magnetic field,
which lead to conflicting interpretations of its origin. Motivated by this
development, we study the magnetic field dependence of thermal conductivity of
a generalized Kitaev model, which allows the phase transitions between
different flux sectors as a function of the magnetic field. The thermal
conductivity due to Majorana fermions shows dip-bump structures as the
magnetic field increases, which is caused by either the transitions between
different flux sectors of Kitaev spin liquids or the topological transitions
that change the Majorana Chern number within the same flux sector. It is shown
that the change of Chern number is closely related to the four-Majorana-
fermion interaction induced by the magnetic field. The non-monotonic behavior
in thermal conductivity emerges at finite temperature, and it becomes weaker
when temperature decreases toward zero. Our model provides a generic mechanism
for the Kitaev spin liquids to develop non-monotonic magnetic-field dependence
of thermal conductivity while the comparison to realistic materials remains an
open question for future investigation.

##  I Introduction

Quantum spin liquid is an exotic phase of quantum magnets that has drawn great
attention since it was initially proposed Anderson (1973). Contrary to
conventional magnetic phases, quantum spin liquids are not magnetically
ordered down to zero temperature, and the spins are fractionalized into spinon
quasiparticles coupling to emergent gauge fields, exhibiting a lot of
interesting collective phenomena including fractionalized excitations, ground
state degeneracy and long range entanglement Wen (2002); Balents (2010);
Savary and Balents (2016); Zhou  _et al._ (2017); Hermanns  _et al._ (2018);
Knolle and Moessner (2019); Broholm  _et al._ (2020). The celebrated Kitaev
honeycomb model Kitaev (2006) is exactly solvable with quantum spin liquid
ground state. Kitaev spin liquids arise from the bond-dependent spin
interactions that frustrate the spin orientation at each site, which can be
generated from $t_{2g}$ orbitals with $j_{\rm{eff}}=\frac{1}{2}$ in materials
with strong spin-orbit coupling Jackeli and Khaliullin (2009). The Kitaev
honeycomb model provides a plausible platform for the experimental realization
of quantum spin liquid, which leads to significant experimental and
theoretical efforts in studying Kitaev materials Jackeli and Khaliullin
(2009); Rau  _et al._ (2014); Rau and Kee (2014); Rau  _et al._ (2016); Trebst
(2017); Kasahara  _et al._ (2018); Gohlke  _et al._ (2018); Takagi  _et al._
(2019); Motome and Nasu (2020); Ye  _et al._ (2020); Fuchs  _et al._ (2020);
Li and Okamoto (2021); Feng  _et al._ (2021), and many Kitaev spin liquid
candidates including (Na,Li)2IrO3 Singh and Gegenwart (2010); Chaloupka  _et
al._ (2010); Liu  _et al._ (2011); Singh  _et al._ (2012); Choi  _et al._
(2012); Modic  _et al._ (2014); Kim  _et al._ (2016) and $\alpha$-RuCl3 Plumb
_et al._ (2014); Sandilands  _et al._ (2015); Kim  _et al._ (2015); Sears  _et
al._ (2015); Majumder  _et al._ (2015); Sandilands  _et al._ (2016); Banerjee
_et al._ (2016); Do  _et al._ (2017); Hirobe  _et al._ (2017) are identified.

Recent thermal transport experiments on Kitaev spin liquid candidate
$\alpha$-RuCl3 have reported non-monotonic dependence of longitudinal thermal
conductivity on magnetic field similar to quantum oscillations Czajka  _et
al._ (2021); Bruin  _et al._ (2021). There exist several explanations for the
reason of this non-monotonic behavior, e.g., from the quantum oscillation from
spinon Fermi surface Czajka  _et al._ (2021) or from phase transitions in the
spin liquid Bruin  _et al._ (2021); Suetsugu  _et al._ (2022). Contrary to
ordinary quantum oscillations in fermionic systems which arise from quantized
Landau levels, the spinons in quantum spin liquids are charge neutral and they
cannot directly couple to orbital magnetic field to give Landau levels.
Therefore, the mechanism of this non-monotonic behavior in thermal
conductivity remains an open question.

Figure 1:  Illustration of terms in the Hamiltonian. A and B sublattices are
labeled as red and blue dots. The dashed arrows in (a),(b),(c) represent
$R_{3}$, $R_{3}^{\prime}$ and $R_{2}$ respectively. (d) shows the relative
locations of the sites in $R_{I}$. The sites at A and B sublattices are
labeled by red and blue dots.

Motivated by these experiments, in this work we explore what mechanism can
lead to non-monotonic dips and bumps in thermal conductivity as a function of
magnetic field. We consider a generalized Kitaev model and demonstrate that
dips and bumps in thermal conductivity can arise naturally due to phase
transitions induced by the magnetic field. Two important ingredients are
included in this generalized Kitaev model, the frustrated further-neighbour
Majorana hopping from spin interactions Zhang  _et al._ (2019, 2020) and the
four-Majorana interaction induced by the magnetic field.

The frustrated further-neighbour Majorana hopping allows different flux
sectors to be the ground state. In Kitaev spin liquids spinons are coupled to
an emergent $\mathbb{Z}_{2}$ gauge field and each plaquette can have
$\mathbb{Z}_{2}$ flux $0$ or $\pi$, where $\pi$ flux corresponds to a vison.
Based on the density and distribution of visons, the quantum spin liquid can
be classified into different flux sectors. Each flux sector has distinct
spinon spectrum and free energy. The original Kitaev honeycomb model has its
ground state in the $0$-flux sector where all plaquettes have $0$ flux by
Lieb’s theorem Lieb (1994). If further-neighbour Kitaev spin interactions are
considered, frustrated third neighbour Majorana hopping terms Zhang  _et al._
(2019, 2020) can arise. These terms commute with the $\mathbb{Z}_{2}$ gauge
field hence they preserve the exact-solubility of the model. With these terms,
different flux sectors can be stabilized as the ground state, which opens up
the possibility for phase transitions between flux sectors. The magnetic field
perturbatively generates both Majorana hopping terms and four-Majorana
interaction terms. We consider both types of terms and applied self-consistent
mean field theory to decouple the Majorana interaction. We found that the
Majorana interaction can effectively rescale the frustrated hopping and lead
to transitions that can change Majorana Chern number.

We take into account these frustrated hopping terms and the four-Majorana
interaction in the generalized Kitaev model. We show that external magnetic
field can lead to transitions between different flux sectors as well as
topological transitions within each flux sector that change the Majorana Chern
number. Because different flux sectors have distinct spinon spectrum and
thermal conductivity, and because thermal conductivity is sensitive to the
spinon band gap which must close when Chern number changes, these phase
transitions will lead to dips and bumps in thermal conductivity as a function
of magnetic field.

This paper is organized as follows. In Section II we establish the generalized
Kitaev model that includes the four-Majorana interaction terms induced by
magnetic field and the frustrated Majorana hopping terms. In Section III we
utilize self-consistent mean field theory to compute Majorana interaction. In
Section IV we present our computation of thermal conductivity and show that it
is a non-monotonic function of magnetic field with dips and bumps. In Section
V we show that this dip-bump feature in thermal conductivity originates from
phase transitions induced by magnetic field. Further relation between our
results and the experiments are discussed in the discussion section.

##  II Generalized Kitaev model with multiple flux sectors

The Hamiltonian of the original Kitaev honeycomb model Kitaev (2006) with
nearest neighbour (NN) spin interactions is given by

$H_{1}=-K_{1}\sum_{\langle
jk\rangle}\sigma_{j}^{\alpha}\sigma_{k}^{\alpha}=iK_{1}\sum_{\langle
jk\rangle}u_{jk}c_{j}c_{k},$ (1)

where $j,k$ denote lattice sites, the spin operator at each site is
represented by $\sigma_{j}^{\alpha}=ib_{j}^{\alpha}c_{j}$, and the
$\mathbb{Z}_{2}$ gauge field $u_{jk}=ib_{j}^{\alpha}b_{k}^{\alpha}$ where
$\alpha\in\\{x,y,z\\}$ is fixed by the direction of bond $jk$. We choose the
gauge such that in the $0$-flux sector $u_{jk}=+1$ if $j$ ($k$) is at $A$ (B)
sublattice, where $A,B$ sublattices are labeled in Fig.1 (a). The original
Kitaev model $H_{1}$ considers nearest neighbour interaction and the ground
state is in the 0-flux sector. In real materials, further-neighbour spin
interactions exist in general. These longer-range interactions are allowed by
symmetry, and they can stabilize different flux sectors as the ground state
and lead to phase transitions between different flux sectors. Therefore, we
consider a generalized Kitaev model with additional third-neighbour spin
interaction terms proposed in Refs Zhang  _et al._ , 2019, 2020:

$\displaystyle H_{3}$ $\displaystyle=$ $\displaystyle K_{3}\sum_{\langle
jklm\rangle\in
R_{3}}\sigma_{j}^{\alpha}\sigma_{k}^{\gamma}\sigma_{l}^{\alpha}\sigma_{m}^{\gamma}$
(2) $\displaystyle-K_{3}^{\prime}\sum_{\langle jklm\rangle\in
R_{3}^{\prime}}\sigma_{j}^{\alpha}\sigma_{k}^{\beta}\sigma_{l}^{\beta}\sigma_{m}^{\alpha}$

Here $R_{3}$ and $R_{3}^{\prime}$ denote the sets of vectors shown in Fig.1
(a) and (b) respectively, and $\langle jklm\rangle\in R_{3}$ means the sites
$j$ and $m$ are connected by the zigzag line passing through $jklm$ with
$\mathbf{r}_{m}-\mathbf{r}_{j}\in R_{3}$. The summation is over all such
zigzag lines. $\alpha$ is along the direction of bond $jk$, $\gamma$ is along
bond $lm$, and $\beta$ is distinct from the directions of $jk$ and $kl$. In
terms of Majorana fermions, $H_{3}$ can be rewritten as

$\displaystyle H_{3}$ $\displaystyle=$ $\displaystyle-iK_{3}\sum_{\langle
jklm\rangle\in R_{3}}u_{jk}u_{kl}u_{lm}c_{j}c_{m}$ $\displaystyle-
iK_{3}^{\prime}\sum_{\langle jklm\rangle\in
R_{3}^{\prime}}u_{jk}u_{kl}u_{lm}c_{j}c_{m}$

These terms commute with gauge field $u_{jk}$, hence they preserve the exact
solubility of the model and the ground state can still be labeled by flux
sectors as the original Kitaev model. Note that each vector in $R_{3}$
corresponds to two zigzag lines (marked by colors in Fig.1 (a)), and the
contribution from these two paths will cancel out if there is a $\pi$-flux
inside the hexagon. Hence the $\pi$-flux sector in which every hexagon has
flux $\pi$ is not affected by the $K_{3}$ term. It has been shown that by
varying $K_{3}$ and $K_{3}^{\prime}$, different flux sectors can be stabilized
as the ground state Zhang  _et al._ (2019, 2020).

The external magnetic field can generate three-spin interactions in the Kitaev
model as follows Kitaev (2006):

$H_{h}={-}h\sum_{\langle jlk\rangle\in
R_{2}}\sigma_{j}^{\alpha}\sigma_{k}^{\beta}\sigma_{l}^{\gamma}{-}h\sum_{\langle
jlk\rangle\in R_{I}}\sigma_{j}^{\alpha}\sigma_{k}^{\beta}\sigma_{l}^{\gamma}$
(3)

Here $R_{2}$ and $R_{I}$ are shown in Fig.1(c) and (d) respectively. The
choice of $\alpha,\beta,\gamma$ for the first term is along $jl,lk,lm$
respectively and for the second term is along $jm,km,lm$ respectively. These
terms arise from the leading order perturbation of magnetic field and
$h\sim\frac{B_{x}B_{y}B_{z}}{K_{1}^{2}}$. When these terms are rewritten as
Majorana fermions, the first term becomes a NNN hopping term and the second
term becomes a four-Majorana interaction term:

$\displaystyle H_{h}$ $\displaystyle=$ $\displaystyle H_{2}+H_{I},\ \
H_{2}={-}ih\sum_{\langle jlk\rangle\in R_{2}}u_{jl}u_{lk}c_{j}c_{k}$ (4)
$\displaystyle H_{I}$ $\displaystyle=$
$\displaystyle{-}h\sum_{\Ydown}u_{jm}u_{km}u_{lm}c_{m}c_{j}c_{k}c_{l}$ (5)
$\displaystyle{-}h\sum_{\Yup}u_{jm}u_{km}u_{lm}c_{m}c_{j}c_{k}c_{l}$

Here the relative locations of sites in $H_{I}$ are shown in Fig.1(d). The
total Hamiltonian for the generalized Kitaev model is

$H=H_{1}+H_{2}+H_{3}+H_{I}$ (6)

This Hamiltonian is symmetric under space-inversion and threefold rotation. It
also has a gauge symmetry indicated by $c_{j}\rightarrow(-1)^{n_{j}}c_{j}$,
$u_{jk}\rightarrow(-1)^{n_{j}-n_{k}}u_{jk}$ where $n_{j}\in\\{0,1\\}$. Time-
reversal maps $h\sim\frac{B_{x}B_{y}B_{z}}{K_{1}^{2}}$ to $-h$, hence time-
reversal symmetry is broken by $H_{2}$ and $H_{I}$ terms induced by the
magnetic field. Without the magnetic field, the Hamiltonian describes a system
of free Majorana fermions for each given flux sector, and different flux
sectors can be stabilized as the ground state when we vary parameters $K_{3}$
and $K_{3}^{\prime}$. The magnetic field breaks time-reversal symmetry and
brings in interactions between Majorana fermions. We show below that the
magnetic field can induce first-order transitions between different flux
sectors, and within the same flux sector the increasing magnetic field can
also induce topological transitions that changes the Chern number of occupied
bands. Both effects can lead to jumps and bumps in physical observables
including thermal conductivity.

##  III Mean field theory for Majorana interactions

We handle the four-Majorana interactions $H_{I}$ in Eq.(5) by self-consistent
mean field theory. For each flux sector, $H_{I}$ contains terms like
$\pm|h|c_{m}c_{j}c_{k}c_{l}$. Denote the center of each $\Yup$ or $\Ydown$ by
site $m$. Note that $\pm|h|c_{m}c_{j}c_{k}c_{l}=-\frac{|h|}{2}(\pm
ic_{m}c_{j}+ic_{k}c_{l})^{2}+|h|$, we can make a Hubbard-Stratonovich (H-S)
transformation to decouple the square term. To preserve the threefold
rotational symmetry we introduce three real auxiliary fields
$\Delta_{jk},\Delta_{kl},\Delta_{lj}$. Define
$\eta_{mjkl}=-\text{sign}(h)u_{jm}u_{km}u_{lm}=\pm 1$, then the H-S
transformation is given by:

$\displaystyle-
hu_{jm}u_{km}u_{lm}c_{m}c_{j}c_{k}c_{l}=|h|\eta_{mjkl}c_{m}c_{j}c_{k}c_{l}$
$\displaystyle\rightarrow\frac{3\Delta_{kl}^{2}}{2|h|}-i\Delta_{kl}(\eta_{mjkl}c_{m}c_{j}+c_{k}c_{l})$
$\displaystyle\ \ \
+\frac{3\Delta_{lj}^{2}}{2|h|}-i\Delta_{lj}(\eta_{mjkl}c_{m}c_{k}+c_{l}c_{j})$
$\displaystyle\ \ \
+\frac{3\Delta_{jk}^{2}}{2|h|}-i\Delta_{jk}(\eta_{mjkl}c_{m}c_{l}+c_{j}c_{k})+|h|$
(7)

Here the factor of 3 is to ensure integrating out the $\Delta$’s recovers the
original four-fermion interaction, and $\eta_{mjkl}$ keeps track of the
overall sign of the four-Majorana term. We apply Eq.(7) to $H_{I}$ with
distinct $\Delta$’s for each term, then we have decoupled the Majorana
interaction into free Majorana coupled to auxiliary fields. We further require
the terms related by a translation of (extended) lattice vector have the same
$\Delta$, then the decoupling preserves the translational symmetry. The number
of independent $\Delta$’s is three times the number of sites in the extended
unit cell for each flux sector. For example, the 0-flux sector has two sites
in the unit cell and it needs six $\Delta$’s. The $\pi$-flux sector is shown
in Fig.2(a), where the signs of the $\mathbb{Z}_{2}$ gauge field $u_{jk}$ at
the thick black bonds are flipped. The extended unit cell of the $\pi$-flux
sector has four sites and it needs twelve $\Delta$’s.

The mean field solution is obtained by finding the set of $\Delta$’s that
minimizes the free energy. This is also equivalent to finding $\Delta$’s that
satisfy the self-consistency equations:

$\displaystyle\Delta_{kl}$ $\displaystyle=$ $\displaystyle\frac{|h|}{3}\langle
i\eta_{mjkl}c_{m}c_{j}+ic_{k}c_{l}\rangle$ $\displaystyle\Delta_{lj}$
$\displaystyle=$ $\displaystyle\frac{|h|}{3}\langle
i\eta_{mjkl}c_{m}c_{k}+ic_{l}c_{j}\rangle$ $\displaystyle\Delta_{jk}$
$\displaystyle=$ $\displaystyle\frac{|h|}{3}\langle
i\eta_{mjkl}c_{m}c_{l}+ic_{j}c_{k}\rangle$ (8)

Here the average denotes the thermal average at the same temperature as the
temperature under which the free energy is computed. It is evident from Eq.(8)
that the $\Delta$’s are real. Our mean field theory is different from Ref Li
and Franz, 2018 in that our decoupling ensures the mean field solution
obtained by solving the self-consistency equations and minimizing free energy
agree with each other. We have computed the mean field solutions for different
flux sectors with flux
$0,\pi,\frac{1}{2}\pi,\frac{1}{3}\pi,\frac{2}{3}\pi,\frac{1}{4}\pi,\frac{3}{4}\pi$.

Our solutions show that the $\Delta$’s respect the symmetries of each flux
sector. For example, the $0$-flux sector has uniform flux distribution in
space, and all the $\Delta$’s in Eq.(7) have the same magnitude in the mean
field solution. We can label each $\Delta$ term by a colored bond with half
the length of a NN bond as in Fig.2(c). The sign of $\Delta$’s are distributed
in such a way that it modifies the NN hopping strength from $K_{1}$ to
$K_{1}+2\Delta$, where the factor of 2 is because each NN bond is shared by
two $\Delta$’s as in Fig.2(c), and the NNN hopping strength is modified from
$h$ to $h+\Delta$. The uniform magnitude of $\Delta$’s also appears in the
$\pi$-flux sector. For other flux sectors whose flux distribution in space is
not uniform, e.g., the $\frac{1}{4}\pi$-flux phase in Fig.2(b), only the
symmetry-related $\Delta$’s have the same magnitude. Here all the $\Delta$’s
at the blue bonds in Fig.2(b) have the same value $\Delta_{a}$, and all
$\Delta$’s at the red bonds have another value $\Delta_{b}$. As the magnetic
field increases, $\Delta$ also increases and the ratio $\Delta/h$ becomes a
constant at large $h$, as shown Fig.2(d). With this mean field solution, we
are able to compute the Majorana spectrum and the evolution of thermal
conductivity with magnetic field.

Figure 2:  (a): $\pi$-flux sector in which every hexagon has flux $\pi$. The
$\mathbb{Z}_{2}$ gauge field $u_{jk}$ flips sign at the thick black bonds to
give rise to $\pi$ flux in the hexagons near it. The extended unit cell is
denoted by the parallelogram. (b): $\frac{1}{4}\pi$-flux sector where
white(gray) hexagons have $0$($\pi$) flux inside it. $u_{jk}$ flips sign at
the thick black bonds. The extended unit cell is denoted by the parallelogram.
In the mean field solution all the $\Delta$’s labeled by the blue (red) bonds
have the same magnitude, which preserves the symmetry of the flux sector. (c):
Different terms in the mean field Hamiltonian with $\Delta$ and
$\Delta^{\prime}$ are labeled by blue and red bonds respectively. (d):
$\Delta$ as a function of $h$ for $0$-, $\pi$\- and $\frac{1}{4}\pi$-flux
sectors computed at $K_{1}=1,K_{3}=0.3,K_{3}^{\prime}=0.34,T=\frac{1}{24}$.
For $\frac{1}{4}\pi$-flux sector the $\Delta$ at the blue bonds in (b) are
used. $\Delta$ increases with $h$ as magnetic field increases, and the
$\Delta$ becomes linear in $h$ when $h$ is large.

##  IV Thermal conductivity

Thermal conductivity can be computed from the mean field Hamiltonian for each
flux sector. The Hamiltonian in momentum space is
$H=\sum_{\mathbf{k}}\psi_{\mathbf{k},i}^{\dagger}H^{\Delta}(\mathbf{k})_{ij}\psi_{\mathbf{k},j}$,
where
$\psi_{\mathbf{k},j}=\frac{1}{{\sqrt{2N}}}\sum_{j}e^{i\mathbf{k}\cdot\mathbf{r}_{j}}c_{j}$
is the Fourier transform of Majorana operators, $H^{\Delta}(\mathbf{k})$ is
the Hamiltonian matrix that also depends on $\Delta$’s. The system has zero
chemical potential, and the heat current operator is given by Katsura  _et
al._ (2010)

$\mathbf{J}^{Q}=\sum_{\mathbf{k}}\psi_{\mathbf{k},i}^{\dagger}\left(\frac{1}{2}\partial_{\mathbf{k}}(H^{\Delta}(\mathbf{k})^{2})\right)_{ij}\psi_{\mathbf{k},j}$
(9)

The current-current correlation in Matsubara frequency is

$\Pi_{\mu\nu}(i\Omega_{n})=-\int_{0}^{\beta}d\tau e^{i\Omega_{n}\tau}\langle
T_{\tau}J_{\mu}^{Q}(\tau)J_{\nu}^{Q}(0)\rangle$ (10)

The thermal conductivity $\kappa_{xx}$ can be computed by Mahan (2013); Durst
and Lee (2000); Wang  _et al._ (2021)

$\displaystyle\kappa_{xx}$ $\displaystyle=$
$\displaystyle-\lim_{\Omega\rightarrow 0}\frac{\text{Im
}\Pi_{xx}(\Omega+i\delta)}{\Omega T}$ (11) $\displaystyle=$
$\displaystyle\lim_{\Omega\rightarrow
0}\pi\int\frac{d^{2}k}{(2\pi)^{2}}\int_{-\infty}^{\infty}d\omega\frac{n_{F}(\omega)-n_{F}(\omega+\Omega)}{\Omega
T}$ $\displaystyle\ \
\text{Tr}\left[J_{x}^{Q}(\mathbf{k})A(\mathbf{k},\omega)J_{x}^{Q}(\mathbf{k})A(\mathbf{k},\omega+\Omega)\right],$

where $n_{F}$ is the Fermi distribution function, $J_{x}^{Q}(\mathbf{k})$ is
the matrix of heat current operator,
$A(\mathbf{k},\omega)=-\frac{1}{\pi}\text{Im}G^{\text{ret}}(\mathbf{k},\omega)=\frac{\gamma}{(\omega-H^{\Delta}(\mathbf{k}))^{2}+\gamma^{2}}$
is the matrix of spectral function and we have added a small impurity
scattering rate $\gamma=0.03K_{1}$. Eq.(11) involves the derivative of $n_{F}$
which is peaked at the Fermi energy with a width proportional to temperature
$T$, hence only the states with energy close to the Fermi level within the
scale of $T$ can contribute to $\kappa_{xx}$. This indicates thermal
conductivity will be small if the density of states near the Fermi energy is
low.

The thermal conductivity as a function of $h/K_{1}$ for $0$\- and $\pi$-flux
sectors with parameters $K_{1}=1,K_{3}=0.3,K_{3}^{\prime}=0.34,T=\frac{1}{24}$
are plotted in Fig.3(a). We only show these two flux sectors because they have
lower free energy than the other flux sectors, and the results for the other
flux sectors are shown in the appendix. Fig.3(a) shows that the $\kappa_{xx}$
curves develop dips and bumps as a function of magnetic field. The detailed
shape of the curves depend on the parameters, but the appearance of dip-bump
features in each flux sector is generic. The free energy with the same
parameters is shown in Fig.3(b). It shows the ground state is in the
$\pi$-flux sector for small $h$ and it goes through a transition to the
$0$-flux sector at large $h$. This first-order transition can lead to a jump
in physical observables. Fig.3(d) shows the thermal conductivity of the flux
sector with the lowest free energy. The red dashed line represents the phase
transition from the $\pi$\- to $0$-flux sector, which is accompanied by a dip
in $\kappa_{xx}$. The temperature dependence of thermal conductivity for $0$\-
and $\pi$ flux sectors are shown in Fig.4, where in (b) and (d) the curves are
normalized by $\kappa_{xx}$ at $h=0$. It shows that the dips and bumps become
weaker at low temperature, because thermal conductivity is a transport
property at finite temperature which will decrease when temperature is low.

Figure 3:  Physical quantities as a function of $h/K_{1}$ for $0$\- and
$\pi$-flux sectors with parameters
$K_{1}=1,K_{3}=0.3,K_{3}^{\prime}=0.34,T=\frac{1}{24}$. We only show the
results for these two flux sectors because they have lower free energy than
other flux sectors. (a) Thermal conductivity. (b) Free energy. (c) Chern
number. (d) Thermal conductivity of the flux sector with the lowest free
energy. The maxima of $\kappa_{xx}$ in (a) coincide with the changes of Chern
number in (c). The red dashed lines in (d) represent the transition from
$\pi$\- to $0$-flux sector and the topological transition that changes Chern
number. These transitions give rise to dip-bump features in thermal
conductivity as a function of magnetic field.  Figure 4:  Thermal conductivity
as a function of $h/K_{1}$ for $0$\- and $\pi$-flux sectors at different
temperature with $K_{1}=1,K_{3}=0.3,K_{3}^{\prime}=0.34$. (a), (b) are for
$0$-flux sector and (c), (d) are for $\pi$-flux sector. (b) and (d) are
rescaled by $\kappa_{xx}$ at $h=0$. The dips and bumps become weaker at low
temperature.

##  V Dips and bumps in thermal conductivity induced by phase transitions

The dips and bumps of thermal conductivity as a function of magnetic field
shown in Fig.3(a) and (d) is a generic feature that exists in a wide range in
the parameter space. These phenomena can be understood from the transition
between different flux sectors and the gap closing accompanied by the
topological transitions that changes the Chern number. When the magnetic field
is small, it tends to open a gap in the Majorana fermion bands for $0$\- and
$\pi$-flux sectors. The low energy band structure $E(\mathbf{k})$ of $0$\- and
$\pi$-flux sectors for the parameters in Fig.3(a) at zero and small $h$ are
shown in Fig.5. The spectrum at $h=0$ without magnetic field is gapless for
both flux sectors. When a small magnetic field is added, the spectrum in both
flux sectors open a gap. From Eq.(11) the thermal conductivity is related to
the density of states near the Fermi energy, therefore this gap reduces the
thermal conductivity, leading to the decrease of $\kappa_{xx}$ at small $h$
for $0$\- and $\pi$-flux sectors in Fig.3(a).

Figure 5:  Dispersion of Majorana bands near the Fermi energy for $0$\- and
$\pi$\- flux sectors with $K_{1}=1,K_{3}=0.3,K_{3}^{\prime}=0.34$. For both
flux sectors, the spectrum are gapless at $h=0$, and a small $h$ induced by
magnetic field can open a gap. The gap in the $\pi$-flux sector is much
smaller than that for the $0$-flux sector.

When magnetic field increases, there is a transition from $\pi$\- to $0$-flux
sector, which leads to a dip in thermal conductivity at $h/K_{1}=0.165$. As
the magnetic field further increases, it can induce topological transitions
that change the Chern number for the Majorana bands. The Chern number for each
flux sector is shown in Fig.3(c), and a comparison with Fig.3(a) shows that
the local maxima of thermal conductivity coincide with the change in Chern
number. This is because the gap needs to close near transitions that change
the Chern number, which leads to an increase in the density of states near
Fermi energy. Because $\kappa_{xx}$ is related to the density of state, this
change in Chern number leads to an increase of $\kappa_{xx}$, which is the
origin of the bump of thermal conductivity near $h/K_{1}=0.3$.

Figure 6:  Chern number of the $0$-flux sector with finite $H_{2}$ and
vanishing $H_{I}$. This phase diagram is independent of the magnitude of
$H_{2}$ as long as $H_{2}$ is finite. Interaction $H_{I}$ induced by the
magnetic field can rescale the effective $K_{3}$ and $K_{3}^{\prime}$ along
the white arrow, leading to a phase transition that changes the Chern number.

The Majorana interaction term $H_{I}$ induced by the magnetic field also plays
an important role in the change of Chern number. Consider the $0$-flux sector
for simplicity in which all the $\Delta$’s have the same value in the mean
field solution. If we turn off interaction $H_{I}$ but keep a finite $H_{2}$
(see Eq.(4)), the phase diagram of Chern number for the $0$-flux sector with
different $K_{3}$ and $K_{3}^{\prime}$ is shown in Fig.6. From Eq.(7) the
auxiliary field $\Delta$ modifies the NN hopping terms in the mean field
Hamiltonian by changing $K_{1}$ to $K_{1}+2\Delta$. Because an overall scaling
of Hamiltonian does not modify the Chern number, a system with parameters
$\\{K_{1}+2\Delta,K_{3},K_{3}^{\prime}\\}$ has the same Chern number as the
one with parameters
$\\{K_{1},K_{3}\frac{K_{1}}{K_{1}+2\Delta},K_{3}^{\prime}\frac{K_{1}}{K_{1}+2\Delta}\\}$.
Therefore, the interaction term effectively rescales $K_{3}$ and
$K_{3}^{\prime}$ by a factor $\frac{K_{1}}{K_{1}+2\Delta}$. The mean field
solution shows $\Delta$ is a positive number that grows with $h$. Therefore,
if the system without magnetic field has $K_{3},K_{3}^{\prime}$ located at,
e.g., the red dot in Fig.6, then as the magnetic field increases, $h$ and
$\Delta$ will increase and the effective $K_{3}$ and $K_{3}^{\prime}$ will
follow the white arrow in Fig.6. When the effective $K_{3}$ and
$K_{3}^{\prime}$ cross a phase boundary between different Chern numbers, the
Chern number will change and the gap will close, leading to an increase in
$\kappa_{xx}$. This explains the bump of $\kappa_{xx}$ in the $0$-flux sector
in Fig.3(a) and (d) at $h/K_{1}=0.3$. As long as $K_{3}$ and $K_{3}^{\prime}$
are chosen to be at $C=4$ or $C=-2$ region close to the phase boundary, the
bump in $K_{xx}$ is expected to occur, which is not sensitive to the exact
choice of $K_{3}$ and $K_{3}^{\prime}$. Therefore, the change of Chern number
induced by magnetic field is a robust feature of this Majorana system which is
independent of the detailed parameter choice.

##  VI Discussions

We have shown that in quantum spin liquids an external magnetic field can
induce dips and bumps in thermal conductivity $\kappa_{xx}$ via phase
transitions between different flux sectors and different Chern numbers. It is
tempting to compare our results to the recent experiments where oscillations
of $\kappa_{xx}$ with magnetic field was observed. The dip-bump features that
we found in $\kappa_{xx}$ become weaker at lower temperature, as shown in
Fig.4, which could be consistent with the report in Ref Czajka  _et al._ ,
2021 of a weakening trend of the oscillatory features at lower temperatures.
Our model is a possible mechanism for the non-monotonic behavior in thermal
conductivity, but it is difficult for a direct comparison between our model
and realistic materials. Our model predicts the change in Chern number
accompanies the local maxima of thermal conductivity, which will lead to a
jump in thermal Hall conductivity $\kappa_{xy}$ because $\kappa_{xy}$ is
related to Chern number, and this change in $\kappa_{xy}$ will be sharper at
lower temperature. The change of $\kappa_{xy}$ near the bump of $\kappa_{xx}$
has not been reported in candidate materials. The relation between
$\kappa_{xy}$ and the non-monotonic behavior of $\kappa_{xx}$ are left as an
open question for future investigation.

##  VII Acknowledgment

This work is supported by the Natural Sciences and Engineering Research
Council of Canada (NSERC) and the Center for Quantum Materials at the
University of Toronto. H.Y.K acknowledges the support by the Canadian
Institute for Advanced Research (CIFAR) and the Canada Research Chairs
Program. Computations were performed on the Niagara supercomputer at the
SciNet HPC Consortium. SciNet is funded by: the Canada Foundation for
Innovation under the auspices of Compute Canada; the Government of Ontario;
Ontario Research Fund - Research Excellence; and the University of Toronto.

##  Appendix A Thermal conductivity in other flux sectors

We have computed the thermal conductivity for various flux sectors including
$0,\pi,\frac{1}{2}\pi,\frac{1}{3}\pi,\frac{2}{3}\pi,\frac{1}{4}\pi$, and
$\frac{3}{4}\pi$. The plots for $0$\- and $\pi$-flux sectors are shown in the
main text. Here we show the results for all these flux sectors in Fig.S1. The
free energy for each flux sector as a function of $h$ induced by magnetic
field is plotted in Fig.S1(a). It shows the ground state with the lowest free
energy is in $0$\- or $\pi$-flux sector. The thermal conductivity
$\kappa_{xx}$ is plotted in Fig.S1(b), which shows the dependence of
$\kappa_{xx}$ on magnetic field is non-monotonic in general within each flux
sector. This non-monotonic behaviour is related to the change in Chern number
as shown in Fig.S1(c). Within each flux sector, each local maximum of
$\kappa_{xx}$ coincides with a topological transition that changes the Chern
number. This is because the change in Chern number accompanies the closure of
Majorana band gap and an increase in density of states, which leads to
increased thermal conductivity. Since the transitions that change the Chern
number within each flux sector generically exist, and the transitions between
different flux sectors driven by free energy can also lead to jumps in thermal
conductivity, these transitions provide a generic mechanism to induce non-
monotonic dependence of thermal conductivity on magnetic field.

Figure S1:  Free energy (a), thermal conductivity (b) and Chern number (c) as
a function of $h$ induced by magnetic field for various flux sectors with
$K_{1}=1,K_{3}=0.3,K_{3}^{\prime}=0.34,T=\frac{1}{24}$. In (c) the overlapped
curves are slightly shifted for better visibility.

## References

  * Anderson (1973) P.W. Anderson, “Resonating valence bonds: A new kind of insulator?” Materials Research Bulletin 8, 153–160 (1973).
  * Wen (2002) Xiao-Gang Wen, “Quantum orders and symmetric spin liquids,” Phys. Rev. B 65, 165113 (2002).
  * Balents (2010) Leon Balents, “Spin liquids in frustrated magnets,” Nature 464, 199–208 (2010).
  * Savary and Balents (2016) Lucile Savary and Leon Balents, “Quantum spin liquids: a review,” Reports on Progress in Physics 80, 016502 (2016).
  * Zhou  _et al._ (2017) Yi Zhou, Kazushi Kanoda, and Tai-Kai Ng, “Quantum spin liquid states,” Rev. Mod. Phys. 89, 025003 (2017).
  * Hermanns  _et al._ (2018) M. Hermanns, I. Kimchi, and J. Knolle, “Physics of the kitaev model: Fractionalization, dynamic correlations, and material connections,” Annual Review of Condensed Matter Physics 9, 17–33 (2018).
  * Knolle and Moessner (2019) J. Knolle and R. Moessner, “A field guide to spin liquids,” Annual Review of Condensed Matter Physics 10, 451–472 (2019).
  * Broholm  _et al._ (2020) C. Broholm, R. J. Cava, S. A. Kivelson, D. G. Nocera, M. R. Norman, and T. Senthil, “Quantum spin liquids,” Science 367, eaay0668 (2020).
  * Kitaev (2006) Alexei Kitaev, “Anyons in an exactly solved model and beyond,” Annals of Physics 321, 2–111 (2006), january Special Issue.
  * Jackeli and Khaliullin (2009) G. Jackeli and G. Khaliullin, “Mott insulators in the strong spin-orbit coupling limit: From heisenberg to a quantum compass and kitaev models,” Phys. Rev. Lett. 102, 017205 (2009).
  * Rau  _et al._ (2014) Jeffrey G. Rau, Eric Kin-Ho Lee, and Hae-Young Kee, “Generic spin model for the honeycomb iridates beyond the kitaev limit,” Phys. Rev. Lett. 112, 077204 (2014).
  * Rau and Kee (2014) Jeffrey G. Rau and Hae-Young Kee, “Trigonal distortion in the honeycomb iridates: Proximity of zigzag and spiral phases in na2iro3,” (2014), arXiv:1408.4811 .
  * Rau  _et al._ (2016) Jeffrey G. Rau, Eric Kin-Ho Lee, and Hae-Young Kee, “Spin-orbit physics giving rise to novel phases in correlated systems: Iridates and related materials,” Annual Review of Condensed Matter Physics 7, 195–221 (2016).
  * Trebst (2017) Simon Trebst, “Kitaev materials,” (2017), arXiv:1701.07056 .
  * Kasahara  _et al._ (2018) Y. Kasahara, T. Ohnishi, Y. Mizukami, O. Tanaka, Sixiao Ma, K. Sugii, N. Kurita, H. Tanaka, J. Nasu, Y. Motome, T. Shibauchi, and Y. Matsuda, “Majorana quantization and half-integer thermal quantum hall effect in a kitaev spin liquid,” Nature 559, 227–231 (2018).
  * Gohlke  _et al._ (2018) Matthias Gohlke, Gideon Wachtel, Youhei Yamaji, Frank Pollmann, and Yong Baek Kim, “Quantum spin liquid signatures in kitaev-like frustrated magnets,” Phys. Rev. B 97, 075126 (2018).
  * Takagi  _et al._ (2019) Hidenori Takagi, Tomohiro Takayama, George Jackeli, Giniyat Khaliullin, and Stephen E. Nagler, “Concept and realization of kitaev quantum spin liquids,” Nature Reviews Physics 1, 264–280 (2019).
  * Motome and Nasu (2020) Yukitoshi Motome and Joji Nasu, “Hunting majorana fermions in kitaev magnets,” Journal of the Physical Society of Japan 89, 012002 (2020).
  * Ye  _et al._ (2020) Mengxing Ye, Rafael M. Fernandes, and Natalia B. Perkins, “Phonon dynamics in the kitaev spin liquid,” Phys. Rev. Research 2, 033180 (2020).
  * Fuchs  _et al._ (2020) Jean-Noël Fuchs, Sourabh Patil, and Julien Vidal, “Parity of chern numbers in the kitaev honeycomb model and the sixteenfold way,” Phys. Rev. B 102, 115130 (2020).
  * Li and Okamoto (2021) Shaozhi Li and Satoshi Okamoto, “Thermal hall effect in the kitaev-heisenberg system with spin-phonon coupling,” (2021), arXiv:2112.06107 .
  * Feng  _et al._ (2021) Kexin Feng, Mengxing Ye, and Natalia B. Perkins, “Temperature evolution of the phonon dynamics in the kitaev spin liquid,” Phys. Rev. B 103, 214416 (2021).
  * Singh and Gegenwart (2010) Yogesh Singh and P. Gegenwart, “Antiferromagnetic mott insulating state in single crystals of the honeycomb lattice material ${\text{na}}_{2}{\text{iro}}_{3}$,” Phys. Rev. B 82, 064412 (2010).
  * Chaloupka  _et al._ (2010) Ji ří Chaloupka, George Jackeli, and Giniyat Khaliullin, “Kitaev-heisenberg model on a honeycomb lattice: Possible exotic phases in iridium oxides ${A}_{2}{\mathrm{iro}}_{3}$,” Phys. Rev. Lett. 105, 027204 (2010).
  * Liu  _et al._ (2011) X. Liu, T. Berlijn, W.-G. Yin, W. Ku, A. Tsvelik, Young-June Kim, H. Gretarsson, Yogesh Singh, P. Gegenwart, and J. P. Hill, “Long-range magnetic ordering in na2iro3,” Phys. Rev. B 83, 220403 (2011).
  * Singh  _et al._ (2012) Yogesh Singh, S. Manni, J. Reuther, T. Berlijn, R. Thomale, W. Ku, S. Trebst, and P. Gegenwart, “Relevance of the heisenberg-kitaev model for the honeycomb lattice iridates ${A}_{2}{\mathrm{iro}}_{3}$,” Phys. Rev. Lett. 108, 127203 (2012).
  * Choi  _et al._ (2012) S. K. Choi, R. Coldea, A. N. Kolmogorov, T. Lancaster, I. I. Mazin, S. J. Blundell, P. G. Radaelli, Yogesh Singh, P. Gegenwart, K. R. Choi, S.-W. Cheong, P. J. Baker, C. Stock, and J. Taylor, “Spin waves and revised crystal structure of honeycomb iridate ${\mathrm{na}}_{2}{\mathrm{iro}}_{3}$,” Phys. Rev. Lett. 108, 127204 (2012).
  * Modic  _et al._ (2014) K. A. Modic, Tess E. Smidt, Itamar Kimchi, Nicholas P. Breznay, Alun Biffin, Sungkyun Choi, Roger D. Johnson, Radu Coldea, Pilanda Watkins-Curry, Gregory T. McCandless, Julia Y. Chan, Felipe Gandara, Z. Islam, Ashvin Vishwanath, Arkady Shekhter, Ross D. McDonald, and James G. Analytis, “Realization of a three-dimensional spin–anisotropic harmonic honeycomb iridate,” Nature Communications 5, 4203 (2014).
  * Kim  _et al._ (2016) Heung-Sik Kim, Yong Baek Kim, and Hae-Young Kee, “Revealing frustrated local moment model for pressurized hyperhoneycomb iridate: Paving the way toward a quantum spin liquid,” Phys. Rev. B 94, 245127 (2016).
  * Plumb  _et al._ (2014) K. W. Plumb, J. P. Clancy, L. J. Sandilands, V. Vijay Shankar, Y. F. Hu, K. S. Burch, Hae-Young Kee, and Young-June Kim, “$\alpha-{\mathrm{rucl}}_{3}$: A spin-orbit assisted mott insulator on a honeycomb lattice,” Phys. Rev. B 90, 041112 (2014).
  * Sandilands  _et al._ (2015) Luke J. Sandilands, Yao Tian, Kemp W. Plumb, Young-June Kim, and Kenneth S. Burch, “Scattering continuum and possible fractionalized excitations in $\alpha\text{$-$}{\mathrm{rucl}}_{3}$,” Phys. Rev. Lett. 114, 147201 (2015).
  * Kim  _et al._ (2015) Heung-Sik Kim, Vijay Shankar V., Andrei Catuneanu, and Hae-Young Kee, “Kitaev magnetism in honeycomb ${\text{rucl}}_{3}$ with intermediate spin-orbit coupling,” Phys. Rev. B 91, 241110 (2015).
  * Sears  _et al._ (2015) J. A. Sears, M. Songvilay, K. W. Plumb, J. P. Clancy, Y. Qiu, Y. Zhao, D. Parshall, and Young-June Kim, “Magnetic order in $\alpha-{\text{rucl}}_{3}$: A honeycomb-lattice quantum magnet with strong spin-orbit coupling,” Phys. Rev. B 91, 144420 (2015).
  * Majumder  _et al._ (2015) M. Majumder, M. Schmidt, H. Rosner, A. A. Tsirlin, H. Yasuoka, and M. Baenitz, “Anisotropic ${\mathrm{ru}}^{3+}4{d}^{5}$ magnetism in the $\alpha-{\mathrm{rucl}}_{3}$ honeycomb system: Susceptibility, specific heat, and zero-field nmr,” Phys. Rev. B 91, 180401 (2015).
  * Sandilands  _et al._ (2016) Luke J. Sandilands, Yao Tian, Anjan A. Reijnders, Heung-Sik Kim, K. W. Plumb, Young-June Kim, Hae-Young Kee, and Kenneth S. Burch, “Spin-orbit excitations and electronic structure of the putative kitaev magnet $\alpha-{\mathrm{rucl}}_{3}$,” Phys. Rev. B 93, 075144 (2016).
  * Banerjee  _et al._ (2016) A. Banerjee, C. A. Bridges, J.-Q. Yan, A. A. Aczel, L. Li, M. B. Stone, G. E. Granroth, M. D. Lumsden, Y. Yiu, J. Knolle, S. Bhattacharjee, D. L. Kovrizhin, R. Moessner, D. A. Tennant, D. G. Mandrus, and S. E. Nagler, “Proximate kitaev quantum spin liquid behaviour in a honeycomb magnet,” Nature Materials 15, 733–740 (2016).
  * Do  _et al._ (2017) Seung-Hwan Do, Sang-Youn Park, Junki Yoshitake, Joji Nasu, Yukitoshi Motome, Yong Seung Kwon, D. T. Adroja, D. J. Voneshen, Kyoo Kim, T.-H. Jang, J.-H. Park, Kwang-Yong Choi, and Sungdae Ji, “Majorana fermions in the kitaev quantum spin system $\alpha$-rucl3,” Nature Physics 13, 1079–1084 (2017).
  * Hirobe  _et al._ (2017) Daichi Hirobe, Masahiro Sato, Yuki Shiomi, Hidekazu Tanaka, and Eiji Saitoh, “Magnetic thermal conductivity far above the néel temperature in the kitaev-magnet candidate $\alpha\text{$-$}{\mathrm{rucl}}_{3}$,” Phys. Rev. B 95, 241112 (2017).
  * Czajka  _et al._ (2021) Peter Czajka, Tong Gao, Max Hirschberger, Paula Lampen-Kelley, Arnab Banerjee, Jiaqiang Yan, David G. Mandrus, Stephen E. Nagler, and N. P. Ong, “Oscillations of the thermal conductivity in the spin-liquid state of $\alpha$-rucl3,” Nature Physics 17, 915–919 (2021).
  * Bruin  _et al._ (2021) J. A. N. Bruin, R. R. Claus, Y. Matsumoto, N. Kurita, H. Tanaka, and H. Takagi, “Robustness of the thermal hall effect close to half-quantization in a field-induced spin liquid state,” (2021), arXiv:2104.12184 .
  * Suetsugu  _et al._ (2022) S. Suetsugu, Y. Ukai, M. Shimomura, M. Kamimura, T. Asaba, Y. Kasahara, N. Kurita, H. Tanaka, T. Shibauchi, J. Nasu, Y. Motome, and Y. Matsuda, “Evidence for the first-order topological phase transition in a kitaev spin liquid candidate $\alpha$-rucl3,” (2022), arXiv:2203.00275 .
  * Zhang  _et al._ (2019) Shang-Shun Zhang, Zhentao Wang, Gábor B. Halász, and Cristian D. Batista, “Vison crystals in an extended kitaev model on the honeycomb lattice,” Phys. Rev. Lett. 123, 057201 (2019).
  * Zhang  _et al._ (2020) Shang-Shun Zhang, Cristian D. Batista, and Gábor B. Halász, “Toward kitaev’s sixteenfold way in a honeycomb lattice model,” Phys. Rev. Research 2, 023334 (2020).
  * Lieb (1994) Elliott H. Lieb, “Flux phase of the half-filled band,” Phys. Rev. Lett. 73, 2158–2161 (1994).
  * Li and Franz (2018) Chengshu Li and Marcel Franz, “Majorana-hubbard model on the honeycomb lattice,” Phys. Rev. B 98, 115123 (2018).
  * Katsura  _et al._ (2010) Hosho Katsura, Naoto Nagaosa, and Patrick A. Lee, “Theory of the thermal hall effect in quantum magnets,” Phys. Rev. Lett. 104, 066403 (2010).
  * Mahan (2013) Gerald D Mahan,  _Many-particle physics_ (Springer Science & Business Media, 2013).
  * Durst and Lee (2000) Adam C. Durst and Patrick A. Lee, “Impurity-induced quasiparticle transport and universal-limit wiedemann-franz violation in d-wave superconductors,” Phys. Rev. B 62, 1270–1290 (2000).
  * Wang  _et al._ (2021) Wen O. Wang, Jixun K. Ding, Brian Moritz, Edwin W. Huang, and Thomas P. Devereaux, “Magnon heat transport in a 2d mott insulator,” (2021), arXiv:2109.01119 .