# Keldysh Theory of Thermal Transport in Multiband Hamiltonians

Luqman Saleem<EMAIL_ADDRESS>Physical Science and Engineering
Division, King Abdullah University of Science and Technology (KAUST), Thuwal
23955, Saudi Arabia  Udo Schwingenschlögl<EMAIL_ADDRESS>Physical Science and Engineering Division, King Abdullah University of Science
and Technology (KAUST), Thuwal 23955, Saudi Arabia  Aurélien Manchon
<EMAIL_ADDRESS>Aix-Marseille Université, CNRS, CINaM, Marseille,
France

###### Abstract

We establish a comprehensive theoretical framework for systems subjected to a
static uniform temperature gradient, employing the non-equilibrium Keldysh-
Dyson formalism. This framework interprets the statistical force due to the
temperature gradient as a mechanical force, utilizing both Luttinger’s scalar
and Moreno-Coleman-Tatara’s vector potentials, which collectively emulate the
gauge invariance stemming from the conservation of energy. Our approach has
the ability to treat heat current and heat magnetization on an equal footing,
thereby extending and generalizing previous formalisms. The derived result for
the thermal conductivity is applied to investigate the thermal characteristics
of Weyl magnons in a stacked honeycomb ferromagnet featuring a trivial
insulator phase, a magnon Chern insulator phase, and three Weyl magnon phases.
Against the expectation from the Berry curvature, the magnon Chern insulator
phase exhibits the highest transverse thermal conductivity.

††preprint: APS/123-QED

##  I Introduction

The theoretical description of quantum transport induced by an electromagnetic
stimulus is well established in condensed matter physics [1, 2, 3, 4, 5]. In
the minimal coupling scheme, the electromagnetic field is inserted into the
Hamiltonian via Peierl’s substitution [6] and the adoption of an appropriate
gauge (e.g., length, velocity, Landau) enables the study of a wealth of
phenomena [7, 8, 9, 10]. In contrast, the theoretical description of quantum
transport under a temperature gradient remains a challenging task due to the
absence of a clear fundamental theory for this type of non-equilibrium force
[11, 12]. This hinders the progress of emergent fields of research aiming at
harvesting thermal transport, such as spin caloritronics [13, 14, 15, 16, 17],
topological magnonics [18, 19, 20, 21, 22, 23], quantum heat transport [24,
25], and topological heat transport [26, 12].

Unlike electromagnetic fields, a temperature gradient gives rise to a
statistical force that cannot be readily incorporated into the Hamiltonian.
Luttinger [11] introduced a “gravitational” scalar potential which enables the
study of thermally induced transport of electrons [3, 27, 28, 29, 30, 31, 32],
magnons [33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45], and phonons [46,
47]. However, this approach results in diverging thermal transport
coefficients in the zero temperature limit for systems without time reversal
symmetry [48, 33, 36]. To address this issue, various attempts were made,
including Ref. [29]’s proposal to subtract the divergenceless equilibrium
magnetization currents from the thermal transport coefficients calculated
using Luttinger’s scalar potential. However, the scaling law employed by this
approach lacks a microscopic justification and cannot account for interaction
effects and disorder.

Shitade [49] presented a gauge-invariant theory of gravity by imposing local
space-time translation symmetry, based on scalar and vector gauge fields. The
scalar gauge field was identified to correspond to Luttinger’s scalar
potential, while the vector gauge field gives rise to tortional electric and
magnetic fields. The author utilized this formalism to compute the thermal
conductivity as a response to the tortional electric field and the heat
magnetization as a response to the tortional magnetic field. However, the
assumption of local space-time translation symmetry was not explicitly
addressed in the context of a temperature gradient. Relying on the Cartan and
Keldysh approaches in curved spacetime, the formalism also often exhibits
technical complexity and lacks intuitive clarity.

Moreno and Coleman [50] represented a temperature gradient by a vector
potential rather than a scalar potential. This approach was later adopted by
Tatara [51, 52] to address the issue of diverging thermal transport
coefficients. It was argued that the vector potential, also known as the
thermal vector potential, in conjunction with the magnetic vector potential
that represents electromagnetic fields, can provide accurate thermal transport
coefficients by considering diamagnetic heat currents. Nonetheless, the
approach proved insufficient to fully address the issue of diverging thermal
transport coefficients in systems without time reversal symmetry due to the
distinct nature of thermal transport. The heat current, which quantifies
entropy changes, must incorporate a component arising from the magnetization
because the thermodynamic expression for heat inherently includes this
component [53, 54, 55, 56]. The magnetization also influences electrically
insulating excitations like phonons (lattice vibrations) and magnons
(elementary magnetic excitations). Typically, such excitations possess non-
zero orbital magnetic moments. In the presence of a temperature gradient these
orbital magnetic moments induce magnetization currents, contributing to the
physical currents. Consequently, an accurate description of thermal transport
necessitates the consideration of both conventional heat currents,
proportional to the temperature gradient, and magnetization currents,
proportional to the magnetization gradient.

The experimental discovery of a magnon-mediated Hall effect [57] has spurred
extensive research into understanding the magnon transport in the presence of
a temperature gradient. Systems featuring Dzyaloshinskii-Moriya interaction
[58, 59] exhibit various intriguing transport phenomena, including the magnon
spin-momentum locking [60], magnon thermal Hall effect in both ferromagnetic
[61, 38, 33, 36] and antiferromagnetic [37, 62, 63, 64, 65, 66, 67, 41]
materials, magnon spin Nernst effect [68, 69, 70, 71, 72, 73, 42, 74],
magnonic magnetization torque [75, 68, 76], and magnonic Edelstein effect
[43]. Additionally, topological phenomena were proposed, such as magnon
topological insulators [77, 38, 39, 78, 79, 18, 80, 81, 82] and Dirac [83],
Floquet [84, 85], and Weyl [86, 87, 88] magnons. Developing a thorough
theoretical framework for describing the magnon transport in the presence of a
temperature gradient is of utmost importance to unlock the potential of
magnonics in next-generation spintronic devices [20].

Unfortunately, most theoretical formulations of magnon transport rely on
Luttinger’s scalar potential. The widely used thermal conductivity formula of
Matsumoto et al. [33] attempts to address the divergence of the thermal
transport coefficients by incorporating a correction factor resembling orbital
magnetization, proportional to $\bm{\mathbf{r}}\times\bm{\mathbf{v}}$.
However, this formula contains a position operator that is not well defined
for extended systems. Furthermore, it primarily focuses on the contribution of
the Berry curvature, neglecting other potential contributions, and is only
applicable to non-interacting magnons. Given the recent surge of research on
interacting magnons [89, 90], a comprehensive theory based on Green’s
functions to account for interaction effects is necessary.

The Baym-Kadanoff [91] or Keldysh-Dyson [92, 93] formalism is a powerful tool
for studying non-equilibrium systems. Its extension to temperature gradient
perturbations will enable the calculation of both heat current and heat
magnetization in both fermionic and bosonic systems. A Green’s function
formalism is also readily adaptable to multiband Hamiltonians and effects
related to the Berry curvature of multiband Bloch states. This will enable a
more accurate description of the thermal transport in interacting and
disordered systems, and pave the way for the design and optimization of
magnonic devices.

This study presents a theoretical framework to calculate linear and nonlinear
responses to a steady and uniform temperature gradient. The temperature
gradient is described by both Luttinger’s scalar and thermal vector
potentials, which together result in thermal electric and magnetic fields.
Just as the charge current is proportional to the applied electric field, the
heat current is proportional to the thermal electric field; and the heat
magnetization is related to the thermal magnetic field. The response function
is then computed within the Keldysh-Dyson formalism in the stationary case.

The present approach has major advantages over previous approaches [29, 33,
34, 35, 36, 37, 49, 40]. Firstly, the use of a vector potential to represent
the temperature gradient eliminates previously neglected diamagnetic heat
currents. Secondly, the heat currents can be expressed in terms of Fermi or
Bose distribution functions, enabling the separation of intrinsic and
extrinsic contributions to the thermal transport coefficients. Thirdly, the
formalism is readily applicable to multiband systems, providing a versatile
tool for investigating thermal transport in magnets with multiple magnetic
sublattices as well as in magnetic heterostructures. Fourthly, the formalism
can be extended to other transport quantities (e.g., particle transport,
nonequilibrium densities) and time-dependent non-uniform problems.

This paper is structured as follows: Section II introduces the thermal vector
potential, its associated Hamiltonian, and the thermal electromagnetic field.
Section III develops a Keldysh-Dyson formalism for non-equilibrium Green’s
functions, expanding them in terms of the thermal electromagnetic field tensor
and establishing a linear response theory. Section IV derives the thermal
conductivity, obtaining the heat current as a response to the thermal electric
field and the heat magnetization as a response to the thermal magnetic field.
It is shown that in the clean limit the result is equivalent to the formula of
Matsumoto et al. [33, 34, 35, 36, 37, 40]. Section V focuses on the thermal
properties of Weyl magnons in a honeycomb ferromagnet. Section VI provides a
summary and discussion of the findings.

##  II Heat Transport and Thermal Vector Potential

The thermal vector potential $\bm{\mathbf{A}}_{T}(t,\bm{\mathbf{r}})$ is
related to Luttinger’s scalar potential $\phi(t,\bm{\mathbf{r}})$ and the
temperature gradient as [51]

$\displaystyle\partial_{t}\bm{\mathbf{A}}_{T}(t,\bm{\mathbf{r}})=\bm{\mathbf{\nabla}}\phi(t,\bm{\mathbf{r}})=\frac{\bm{\mathbf{\nabla}}T}{T},$
(1)

with $\partial_{t}\equiv\partial/\partial t$. Unlike the electromagnetic
theory, where the magnetic vector potential possesses U(1) gauge symmetry,
$\bm{\mathbf{A}}_{T}(t,\bm{\mathbf{r}})$ does not possess gauge symmetry in
the formal sense because the energy conservation arises from the global
translation symmetry with respect to time. However, it can be shown that it
couples to the energy in minimal form, i.e.,
$\bm{\mathbf{p}}\rightarrow\bm{\mathbf{p}}-\bm{\mathbf{A}}_{T}\epsilon_{\bm{\mathbf{k}}}$,
where $\bm{\mathbf{p}}$ is the canonical momentum and
$\epsilon_{\bm{\mathbf{k}}}$ is the energy in momentum space [51]. In real
space-time operator form this relation becomes
$\bm{\mathbf{p}}\rightarrow\bm{\mathbf{p}}-i\hbar\bm{\mathbf{A}}_{T}\partial_{t}$.
The minimal form can be attributed to a “gauge invariance” resulting from the
energy conservation law. In fact, one can attribute a part of the temperature
gradient to Luttinger’s scalar potential and the other part to the thermal
vector potential as

$\frac{\bm{\mathbf{\nabla}}T}{T}=\bm{\mathbf{\nabla}}\phi(t,\bm{\mathbf{r}})+\partial_{t}\bm{\mathbf{A}}_{T}(t,\bm{\mathbf{r}}).$
(2)

Consequently, similar to the electromagnetic theory, we obtain gauge
invariance under the transformations $\phi\to\phi+d\chi/dt$ and
$\bm{\mathbf{A}}_{T}\to\bm{\mathbf{A}}_{T}-\bm{\mathbf{\nabla}}\chi$, where
$\chi$ is a scalar field [51]. To elevate the thermal theory to the level of
the electromagnetic theory, we define a thermal electric field

$\bm{\mathbf{E}}_{T}(t,\bm{\mathbf{r}})\equiv-\partial_{t}\bm{\mathbf{A}}_{T}(t,\bm{\mathbf{r}})-\bm{\mathbf{\nabla}}\phi(t,\bm{\mathbf{r}})$
(3)

and a thermal magnetic field

$\bm{\mathbf{B}}_{T}(t,\bm{\mathbf{r}})\equiv\bm{\mathbf{\nabla}}\times\bm{\mathbf{A}}_{T}(t,\bm{\mathbf{r}}).$
(4)

$\bm{\mathbf{E}}_{T}(t,\bm{\mathbf{r}})$ gives rise to the mechanical force
due to the temperature gradient. The heat current density can be defined as
quantum expectation value

$\bm{\mathbf{J}}_{Q}(t,\bm{\mathbf{r}})=\pm
i\text{Tr}\big{(}\hat{\bm{\mathbf{J}}}_{Q}(t,\bm{\mathbf{r}})\hat{G}_{\bm{\mathbf{E}}_{T}}^{<}(t,\bm{\mathbf{r}})\big{)}.$
(5)

Here and in the following, the upper sign is for Bosons and the lower sign is
for Fermions. $\hat{\bm{\mathbf{J}}}_{Q}(t,\bm{\mathbf{r}})$ is the heat
current density operator defined as product of the velocity operator and time
derivative [49] and $\hat{G}_{\bm{\mathbf{E}}_{T}}^{<}(t,\bm{\mathbf{r}})$ is
the lesser component of the Green’s function (proportional to
$\bm{\mathbf{E}}_{T}(t,\bm{\mathbf{r}})$). The hat sign indicates the matrix
of operators in the internal degrees of freedom. Moreover, a thermodynamical
definition can be constructed for the heat magnetization $\bm{\mathbf{M}}_{Q}$
by noticing that the thermal vector potential couples to the heat current
density [94],

$\displaystyle\bm{\mathbf{M}}_{Q}=-\left(\frac{\partial\Omega}{\partial\bm{\mathbf{B}}_{T}}\right)_{\mu,T}.$
(6)

$\Omega=E-TS-\mu N$ is the free energy with $E$, $S$, $\mu$, and $N$ denoting
the internal energy, entropy, chemical potential, and particle number,
respectively.

##  III Keldysh-Dyson formalism for the thermal vector potential

We start by defining a general unperturbed Hamiltonian for Bosons or Fermions
as

$\displaystyle{\hat{H}}_{0}=\int
d\bm{\mathbf{r}}{\Psi}^{\dagger}(\bm{\mathbf{r}})\left({\hat{h}}_{0}(\bm{\mathbf{p}})+\hat{V}(\bm{\mathbf{r}})\right)\Psi(\bm{\mathbf{r}}),$
(7)

where
${\hat{h}}_{0}(\bm{\mathbf{p}})={\left(-i\hbar\bm{\mathbf{\nabla}}\right)}^{2}/2m$
is the non-interacting kinetic part, $\hat{V}(\bm{\mathbf{r}})$ is a general
one-particle potential, and $\Psi^{\dagger}(\bm{\mathbf{r}})$ and
$\Psi(\bm{\mathbf{r}})$ are the creation and annihilation field operators,
respectively.

Moving forward, we adopt a four-vector notation: Time and space coordinates
are represented as $x\equiv x^{\mu}=(t,\bm{\mathbf{r}})$ and
$x_{\mu}=(-t,\bm{\mathbf{r}})$, derivatives are represented as
${\partial}_{x^{\mu}}=({\partial}_{t},\bm{\mathbf{\nabla}})$ and
${\partial}_{x_{\mu}}=(-{\partial}_{t},\bm{\mathbf{\nabla}})$, and potentials
are represented as $A^{\mu}_{T}(x)=(\phi(x),{\bm{\mathbf{A}}}_{T}(x))$ and
${{(A_{T})}_{\mu}}(x)=(-\phi(x),{\bm{\mathbf{A}}}_{T}(x))$. We proceed by
applying a temperature gradient via the thermal electromagnetic tensor

$F^{\mu\nu}_{T}(x)={\partial}_{x_{\mu}}A^{\nu}_{T}(x)-{\partial}_{x_{\nu}}A^{\mu}_{T}(x).$
(8)

A “gauge-invariant” Hamiltonian that incorporates both Luttinger’s scalar and
thermal vector potentials is formulated as

$\displaystyle\hat{H}(t)=$ $\displaystyle\int
d\bm{\mathbf{r}}\Psi^{\dagger}(x)\left(\hat{h}_{0}(\bm{\mathbf{p}}-i\hbar\bm{\mathbf{A}}_{T}\partial_{t})+\hat{V}(\bm{\mathbf{r}})\right)\Psi(x)$
$\displaystyle+i\hbar\int
d\bm{\mathbf{r}}\phi(x)\Psi^{\dagger}(x)\partial_{t}\Psi(x).$ (9)

The second term is Luttinger’s Hamiltonian [11]. While $\hat{H}(t)$ does not
include interactions such as the particle-particle interaction, the following
formalism still applies even when these interactions are present. However, a
self-energy correction according to the Feynman rule is necessary.

###  III.1 Green’s functions

The non-equilibrium Green’s functions are defined on the Keldysh contour $C$
as [4]

$\displaystyle\mathcal{G}(x_{1},x_{2})=-i\big{\langle}T_{C}{\left(\Psi(x_{1})\Psi^{\dagger}(x_{2})\right)}\big{\rangle}.$
(10)

$T_{C}$ is the time-ordering operator. The Dyson equations are

$\displaystyle\left(i\hbar
D_{t_{1}}-{\hat{h}}_{0}(-i\hbar\bm{\mathbf{D}}_{x_{1}})\right)\mathcal{G}(x_{1},x_{2})$
$\displaystyle-\int
dx_{3}\Sigma(x_{1},x_{3})\mathcal{G}(x_{3},x_{2})=\delta(x_{1}-x_{2}),$ (11a)
$\displaystyle\mathcal{G}(x_{1},x_{2})\left(-i\hbar
D_{t_{2}}-{\hat{h}}_{0}(i\hbar\bm{\mathbf{D}}_{x_{2}})\right)$
$\displaystyle-\int
dx_{3}\mathcal{G}(x_{1},x_{3})\Sigma(x_{3},x_{2})=\delta(x_{1}-x_{2}),$ (11b)

where we define a gauge covariant derivative analogous to the electromagnetic
theory as

$D_{{(x_{n})}_{\mu}}=\partial_{{(x_{n})}_{\mu}}+A_{T}^{\mu}(x_{n})\partial_{t_{n}}.$
(12)

In real-time representation both the Green’s function $\mathcal{G}$ and the
self energy $\Sigma$ have three independent components,

$\displaystyle\mathcal{G}(x_{1},x_{2})$
$\displaystyle=\left(\begin{array}[]{cc}\hat{G}^{R}(x_{1},x_{2})&2\hat{G}^{<}(x_{1},x_{2})\\\
0&\hat{G}^{A}(x_{1},x_{2})\end{array}\right),$ (13c)
$\displaystyle\Sigma(x_{1},x_{2})$
$\displaystyle=\left(\begin{array}[]{cc}\hat{\Sigma}^{R}(x_{1},x_{2})&2\hat{\Sigma}^{<}(x_{1},x_{2})\\\
0&\hat{\Sigma}^{A}(x_{1},x_{2})\end{array}\right).$ (13f)

The lesser, advanced, and retarded Green’s functions are defined as

$\displaystyle\begin{aligned} \hat{G}^{<}(x_{1},x_{2})&=\mp
i\left\langle{\Psi}^{\dagger}(x_{2}){\Psi}(x_{1})\right\rangle,\\\
\hat{G}^{A}(x_{1},x_{2})&=+i\theta(t_{2}-t_{1})\left\langle{\left[{\Psi}(x_{1}),{\Psi}^{\dagger}(x_{2})\right]}_{\mp}\right\rangle,\\\
\hat{G}^{R}(x_{1},x_{2})&=-i\theta(t_{1}-t_{2})\left\langle{\left[{\Psi}(x_{1}),{\Psi}^{\dagger}(x_{2})\right]}_{\mp}\right\rangle.\end{aligned}$
(14)

The lesser component is directly connected to the density matrix and can be
used to calculate expectation values of operators.

In a compact form, the Dyson equations (11a) and (11b) read

$\displaystyle\left(\left((\mathcal{G}^{(0)})^{-1}-\Sigma\right)\star\mathcal{G}\right)(x_{1},x_{2})$
$\displaystyle=\delta(x_{1}-x_{2}),$ (15a)
$\displaystyle\left(\mathcal{G}\star\left((\mathcal{G}^{(0)})^{-1}-\Sigma\right)\right)(x_{1},x_{2})$
$\displaystyle=\delta(x_{1}-x_{2})$ (15b)

with the convolution star product defined as

$\displaystyle(\mathcal{A}\star\mathcal{B})(x_{1},x_{2})=\int
dx_{3}\mathcal{A}(x_{1},x_{3})\mathcal{B}(x_{3},x_{2})$ (16)

and

$\displaystyle\begin{aligned}
(\mathcal{G}^{(0)})^{-1}(x_{1},x_{2})=\left(i\hbar
D_{t_{1}}-{\hat{h}}_{0}(-i\hbar\bm{\mathbf{D}}_{x_{1}})\right)\delta(x_{1}-x_{2}),\\\
=\left(-i\hbar
D_{t_{2}}-{\hat{h}}_{0}(i\hbar\bm{\mathbf{D}}_{x_{2}})\right)\delta(x_{1}-x_{2}).\end{aligned}$
(17)

We transform the $(x_{1},x_{2})$ coordinates into the relative microscopic $x$
and center of mass macroscopic $X\equiv X^{\mu}$ coordinates using the Wigner
representation $(X,x)=\left(1/2(x_{1}+x_{2}),x_{1}-x_{2}\right)$. We continue
to the $(X,p)$ coordinates by Fourier transforming $x^{\mu}$ and $x_{\mu}$ to
$p^{\mu}\equiv(\epsilon,\bm{\mathbf{p}})$ and
$p_{\mu}\equiv(-\epsilon,\bm{\mathbf{p}})$, respectively. It is known that in
the Wigner representation the convolution star product becomes the Moyal
product [95]

$\displaystyle(\mathcal{A}\star\mathcal{B})(X,p)=\mathcal{A}(X,p)\text{e}^{i\hbar/2\left(\partial_{X^{\mu}}\partial_{p_{\mu}}-\partial_{p_{\mu}}\partial_{X^{\mu}}\right)}\mathcal{B}(X,p).$
(18)

To write the Dyson equations in the $(X,p)$ coordinates, we calculate
$(\mathcal{G}^{(0)})^{-1}(X,p)$ by defining

$\displaystyle\pi_{\mu}(x_{1},x_{2})=-i\hbar
D_{x_{1}^{\mu}}\delta(x_{1}-x_{2})=i\hbar D_{x_{2}^{\mu}}\delta(x_{1}-x_{2})$
(19)

and Fourier transforming

$\displaystyle\pi_{\mu}(X,p)$ $\displaystyle=\int
dxe^{-ip_{\mu}x^{\mu}/\hbar}\pi_{\mu}(X,x)$
$\displaystyle=p_{\mu}-\epsilon(A_{T})_{\mu}(X,0)+i\hbar\left(\partial_{t}(A_{T})_{\mu}(X,t)\right)_{t=0}.$
(20)

The last term is neglected, as we are interested in the behavior of the system
for $t>0$. $\pi_{\mu}(X,p)$ is the mechanical momentum (which is different
from the canonical momentum $p_{\mu}$).

Finally, we have in $(X,p)$ coordinates

$\displaystyle(\mathcal{G}^{(0)})^{-1}(\pi(X,p))={\pi}^{0}-{\hat{h}^{*}}_{0}(\bm{\mathbf{\pi}}),$
(21)

where ${\hat{h}^{*}}_{0}(\bm{\mathbf{\pi}})$ is the generalization of
${\hat{h}}_{0}(\bm{\mathbf{p}})$ in the presence of Luttinger’s scalar and
thermal vector potentials. It is obtained by replacing all direct products by
star products [5]. The final form of the Dyson equations in terms of the Moyal
product is

$\displaystyle{\left(\left((\mathcal{G}^{(0)})^{-1}-\Sigma\right)\star\mathcal{G}\right)}(X,p)$
$\displaystyle=1,$ (22a)
$\displaystyle{\left(\mathcal{G}\star\left((\mathcal{G}^{(0)})^{-1}-\Sigma\right)\right)}(X,p)$
$\displaystyle=1.$ (22b)

###  III.2 Static uniform thermal field

The off-diagonal elements of $F_{T}^{\mu\nu}$ are the components of the
thermal electric and magnetic fields. We decompose $A_{T}^{\mu}(X,0)$ into
${\tilde{A}}^{\mu}_{T}(X,0)$, which covers the static uniform parts of the
thermal electric and magnetic fields
(${\tilde{F}}^{\mu\nu}_{T}={\partial}_{X_{\mu}}{\tilde{A}}^{\nu}_{T}-{\partial}_{X_{\nu}}{\tilde{A}}^{\mu}_{T}=\text{constant}$)
and $\delta A^{\mu}_{T}(X,0)=A^{\mu}_{T}(X,0)-{\tilde{A}}^{\mu}_{T}(X,0)$,
which covers the dynamic and the non-uniform parts. Similarly, the static
uniform mechanical momentum is
${\tilde{\pi}}^{\mu}(X,p)=p^{\mu}-\epsilon\tilde{A}_{T}^{\mu}(X,0)$. Changing
the variables of the Moyal product (18) and Dyson equations (22a) and (22b)
from the canonical momentum $p$ to the static uniform kinetic momentum
$\tilde{\pi}$ yields

$\displaystyle{\left(\mathcal{A}\star\mathcal{B}\right)}(X,\tilde{\pi})=\mathcal{A}(X,\tilde{\pi}){\text{e}^{i\hbar/2\left(\partial_{X^{\mu}}\partial_{{\tilde{\pi}}_{\mu}}-\partial_{{\tilde{\pi}}_{\mu}}\partial_{X^{\mu}}+{\tilde{\pi}}^{0}{\tilde{F}}^{\mu\nu}_{T}\partial_{{\tilde{\pi}}^{\mu}}\partial_{{\tilde{\pi}}^{\nu}}\right)}}\mathcal{B}(X,\tilde{\pi}),$
(23a)
$\displaystyle\left({\tilde{\pi}}^{0}-\epsilon\delta\phi(X,0)-{\hat{h}}_{0^{*}}(\tilde{\bm{\mathbf{\pi}}}-\epsilon\delta{\bm{\mathbf{A}}}_{T}(X,0))-\Sigma(X,\tilde{\pi})\right)\star\mathcal{G}(X,\tilde{\pi})=1,$
(23b)
$\displaystyle\mathcal{G}(X,\tilde{\pi})\star\left({\tilde{\pi}}^{0}-\epsilon\delta\phi(X,0)-{\hat{h}}_{0^{*}}(\tilde{\bm{\mathbf{\pi}}}-\epsilon\delta{\bm{\mathbf{A}}}_{T}(X,0))-\Sigma(X,\tilde{\pi})\right)=1.$
(23c)

We address in the following only static uniform thermal electromagnetic
fields, i.e., $\delta A_{T}(X,0)=0$, meaning that we do not need to
distinguish between $\pi$ and $\tilde{\pi}$. Furthermore, we restrict
ourselves to steady-state thermal electromagnetic fields, which simplifies
(23a), (23b), and (23c) to

$\displaystyle(\mathcal{A}\star\mathcal{B})(\pi)=\mathcal{A}(\pi){\text{e}^{(i\hbar/2){\pi}^{0}F^{\mu\nu}_{T}\partial_{{\pi}^{\mu}}\partial_{{\pi}^{\nu}}}}$
$\displaystyle\mathcal{B}(\pi),$ (24a)
$\displaystyle\left({\pi}^{0}-{\hat{h}}_{0}(\bm{\mathbf{\pi}})-\Sigma(\pi)\right)\star\mathcal{G}(\pi)$
$\displaystyle=1,$ (24b)
$\displaystyle\mathcal{G}(\pi)\star\left({\pi}^{0}-{\hat{h}}_{0}(\bm{\mathbf{\pi}})-\Sigma(\pi)\right)$
$\displaystyle=1.$ (24c)

From now on, all quantities are functions of $\pi$ unless otherwise stated.

###  III.3 Expansion of the Dyson equations

While often the Dyson equations (24b) and (24c) cannot be solved exactly, we
can perturbatively expand them using the Moyal product (24a), applying the
Green’s function
${\mathcal{G}}_{0}={({\pi}^{0}-{\hat{h}}_{0}(\bm{\mathbf{\pi}})-{\Sigma}_{0})}^{-1}$
from the left (right) to the first (second) equation, and using the identity
${\mathcal{G}}^{-1}_{0}{\mathcal{G}}_{0}={\mathcal{G}}_{0}{\mathcal{G}}^{-1}_{0}=1$,

$\displaystyle\mathcal{G}$
$\displaystyle={\mathcal{G}}_{0}\left(1+(\Sigma-{\Sigma}_{0})\mathcal{G}-\sum^{\infty}_{n=1}\frac{1}{n!}\left({\pi}^{0}-{\hat{h}}_{0}(\bm{\mathbf{\pi}})-\Sigma\right)\left(\prod^{n}_{i=1}\frac{i\hbar{\pi}^{0}}{2}F^{{\mu}_{i}{\nu}_{i}}_{T}\partial_{{\pi}^{{\mu}_{i}}}\partial_{{\pi}^{{\nu}_{i}}}\right)\mathcal{G}\right),$
(25a) $\displaystyle\mathcal{G}$
$\displaystyle=\left(1+\mathcal{G}(\Sigma-{\Sigma}_{0})-\sum^{\infty}_{n=1}\frac{1}{n!}\mathcal{G}\left(\prod^{n}_{i=1}\frac{i\hbar{\pi}^{0}}{2}F^{{\mu}_{i}{\nu}_{i}}_{T}\partial_{{\pi}^{{\mu}_{i}}}\partial_{{\pi}^{{\nu}_{i}}}\right)\left({\pi}^{0}-{\hat{h}}_{0}(\bm{\mathbf{\pi}})-\Sigma\right)\right){\mathcal{G}}_{0}.$
(25b)

Here and in the following, the subscript $0$ in the Green’s functions and
self-energies indicates that the quantity is taken in the absence of external
fields. To expand Eqs. (25a) and (25b) in powers of $\hbar F^{\mu\nu}_{T}/2$
we aim for solutions of the form

$\displaystyle\mathcal{G}$
$\displaystyle={\mathcal{G}}_{0}+\sum^{\infty}_{n=1}\frac{1}{n!}\left(\prod^{n}_{i=1}\frac{\hbar}{2}F^{{\mu}_{i}{\nu}_{i}}_{T}\right){\mathcal{G}}_{{\mu}_{1}{\nu}_{1},\cdots,{\mu}_{n}{\nu}_{n}},$
(26a) $\displaystyle\Sigma$
$\displaystyle={\Sigma}_{0}+\sum^{\infty}_{n=1}\frac{1}{n!}\left(\prod^{n}_{i=1}\frac{\hbar}{2}F^{{\mu}_{i}{\nu}_{i}}_{T}\right){\Sigma}_{{\mu}_{1}{\nu}_{1},\cdots,{\mu}_{n}{\nu}_{n}}.$
(26b)

We find ${\mathcal{G}}_{{\mu}_{1}{\nu}_{1},\cdots,{\mu}_{n}{\nu}_{n}}$ by
inserting these equations into Eq. (25a) or (25b) and comparing the left and
right hand sides. We can calculate by Eq. (13) the retarded, advanced, and
lesser Green’s functions up to any power of $\hbar F^{\mu\nu}_{T}/2$. Below,
we explicitly give the zeroth order (equilibrium Green’s function theory) and
first order (linear response theory) results.

The equilibrium Green’s function theory is derived from the equilibrium Dyson
equation

${\left({\epsilon}-{\hat{h}}_{0}(\bm{\mathbf{p}})-{\Sigma}_{0}(p)\right)}{\mathcal{G}}_{0}(p)=1$
(27)

as

$\displaystyle G^{R(A)}_{0}(p)$
$\displaystyle={\left(\epsilon-{\hat{h}}_{0}(\bm{\mathbf{p}})-{\Sigma}^{R(A)}_{0}(p)\right)}^{-1},$
(28a) $\displaystyle G^{<}_{0}(p)$
$\displaystyle=\mp\left(G^{A}_{0}(p)-G^{R}_{0}(p)\right)f_{\mp},$ (28b)
$\displaystyle{\Sigma}^{<}_{0}(p)$
$\displaystyle=\mp\left({\Sigma}^{A}_{0}(p)-{\Sigma}^{R}_{0}(p)\right)f_{\mp},$
(28c)

where $f_{\mp}(\epsilon)={(e^{\beta\epsilon}\mp 1)}^{-1}$ denotes the Bose and
Fermi distribution functions.

The linear response theory is derived by linearizing Eq. (25a) or (25b) in
$\hbar F^{\mu\nu}_{T}/2$,

$\displaystyle{\mathcal{G}}_{\mu\nu}={\mathcal{G}}_{0}{\Sigma}_{\mu\nu}{\mathcal{G}}_{0}+i{\pi}^{0}{\mathcal{G}}_{0}\left({\partial}_{{\pi}^{\mu}}{\mathcal{G}}^{-1}_{0}\right){\mathcal{G}}_{0}\left(\partial_{{\pi}^{\nu}}{\mathcal{G}}^{-1}_{0}\right){\mathcal{G}}_{0}.$
(29)

${\mathcal{G}}_{\mu\nu}$ is antisymmetric under the exchange of $\mu$ and
$\nu$, because $F^{\mu\nu}_{T}$ is antisymmetric under this exchange. By
rewriting Eq. (29) as

$\displaystyle{\mathcal{G}}_{\mu\nu}$
$\displaystyle={\mathcal{G}}_{0}{\Sigma}_{\mu\nu}{\mathcal{G}}_{0}+\frac{i{\pi}^{0}}{2}\mathcal{G}_{0}\left({\partial}_{{\pi}^{\mu}}{\mathcal{G}}^{-1}_{0}\right){\mathcal{G}}_{0}\left({\partial}_{{\pi}^{\nu}}{\mathcal{G}}^{-1}_{0}\right)\mathcal{G}_{0}$
$\displaystyle-\frac{i{\pi}^{0}}{2}\mathcal{G}_{0}\left({\partial}_{{\pi}^{\nu}}{\mathcal{G}}^{-1}_{0}\right){\mathcal{G}}_{0}\left({\partial}_{{\pi}^{\mu}}{\mathcal{G}}^{-1}_{0}\right)\mathcal{G}_{0}$
(30)

we obtain the thermal electric and magnetic Green’s functions

$\displaystyle{\mathcal{G}}_{{\bm{\mathbf{E}}}_{T}}$
$\displaystyle={\mathcal{G}}_{0}{\Sigma}_{{\bm{\mathbf{E}}}_{T}}{\mathcal{G}}_{0}+\frac{i{\pi}^{0}}{2}\mathcal{G}_{0}\left({\partial}_{{\pi}^{0}}{\mathcal{G}}^{-1}_{0}\right){\mathcal{G}}_{0}\left(\bm{\mathbf{\nabla_{\pi}}}{\mathcal{G}}^{-1}_{0}\right)\mathcal{G}_{0}$
$\displaystyle-\frac{i{\pi}^{0}}{2}\mathcal{G}_{0}\left(\bm{\mathbf{\nabla_{\pi}}}{\mathcal{G}}^{-1}_{0}\right){\mathcal{G}}_{0}\left({\partial}_{{\pi}^{0}}{\mathcal{G}}^{-1}_{0}\right)\mathcal{G}_{0},$
(31a) $\displaystyle{\mathcal{G}}_{{\bm{\mathbf{B}}}_{T}}$
$\displaystyle={\mathcal{G}}_{0}\left({\Sigma}_{{\bm{\mathbf{B}}}_{T}}+\frac{i{\pi}^{0}}{2}\left(\bm{\mathbf{\nabla_{\pi}}}{\mathcal{G}}^{-1}_{0}\times{\mathcal{G}}_{0}\left(\bm{\mathbf{\nabla_{\pi}}}{\mathcal{G}}^{-1}_{0}\right)\right)\right){\mathcal{G}}_{0}.$
(31b)

We can decomposed these matrix Green’s functions into their components using
the fact that for a product
$\mathcal{G}={\mathcal{G}}_{1}{\mathcal{G}}_{2}\cdots{\mathcal{G}}_{n}$ of
matrix Green’s function the retarded, advanced, and lesser components are

$\displaystyle\hat{G}^{R(A)}$
$\displaystyle=\hat{G}^{R(A)}_{1}\hat{G}^{R(A)}_{2}\cdots\hat{G}^{R(A)}_{n},$
(32) $\displaystyle\hat{G}^{<}$
$\displaystyle=\hat{G}^{<}_{1}\hat{G}^{A}_{2}\cdots\hat{G}^{A}_{n}+\hat{G}^{R}_{1}\hat{G}^{<}_{2}\cdots\hat{G}^{A}_{n}+\cdots+\hat{G}^{R}_{1}\hat{G}^{R}_{2}\cdots\hat{G}^{<}_{n}.$

Employing the relations

$\displaystyle\begin{split}\bm{\mathbf{\nabla_{\pi}}}\hat{G}^{R}_{0}=&-\hat{G}^{R}_{0}\bm{\mathbf{\nabla_{\pi}}}({\hat{G}_{0}}^{R})^{-1}\hat{G}^{R}_{0},\\\
\partial_{\pi^{0}}\hat{G}^{R}_{0}=&-\hat{G}^{R}_{0}\partial_{\pi^{0}}({\hat{G}_{0}}^{R})^{-1}\hat{G}^{R}_{0},\end{split}$
(33)

we obtain

$\displaystyle\hat{G}^{R(A)}_{{\bm{\mathbf{E}}}_{T}}=$ $\displaystyle\
\hat{G}^{R(A)}_{0}\left({\hat{\Sigma}}^{R(A)}_{{\bm{\mathbf{E}}}_{T}}+\frac{i{\pi}^{0}}{2}\left({\partial}_{{\pi}^{0}}\left({\hat{G}_{0}}^{R(A)}\right)^{-1}\hat{G}^{R(A)}_{0}\bm{\mathbf{\nabla_{\pi}}}\left({\hat{G}_{0}}^{R(A)}\right)^{-1}-\bm{\mathbf{\nabla_{\pi}}}\left({\hat{G}_{0}}^{R(A)}\right)^{-1}\hat{G}^{R(A)}_{0}{\partial}_{{\pi}^{0}}\left({\hat{G}_{0}}^{R(A)}\right)^{-1}\right)\right)\hat{G}^{R(A)}_{0},$
(34a) $\displaystyle\hat{G}^{R(A)}_{{\bm{\mathbf{B}}}_{T}}=$ $\displaystyle\
\hat{G}^{R(A)}_{0}\left({\hat{\Sigma}}^{R(A)}_{{\bm{\mathbf{B}}}_{T}}+\frac{i{\pi}^{0}}{2}\left(\bm{\mathbf{\nabla_{\pi}}}\left({\hat{G}_{0}}^{R(A)}\right)^{-1}\times\hat{G}^{R(A)}_{0}\bm{\mathbf{\nabla_{\pi}}}\left({\hat{G}_{0}}^{R(A)}\right)^{-1}\right)\right)\hat{G}^{R(A)}_{0},$
(34b) $\displaystyle\hat{G}^{<}_{{\bm{\mathbf{E}}}_{T}}=$
$\displaystyle\mp\left(\hat{G}^{A}_{{\bm{\mathbf{E}}}_{T}}-\hat{G}^{R}_{{\bm{\mathbf{E}}}_{T}}\right)f_{\mp}+\hat{G}^{R}_{0}\left({\hat{\Sigma}}^{<}_{\bm{\mathbf{E}}_{T}}\pm\left({\hat{\Sigma}}^{A}_{\bm{\mathbf{E}}_{T}}-{\hat{\Sigma}}^{R}_{\bm{\mathbf{E}}_{T}}\right)f_{\mp}\right)\hat{G}^{A}_{0}$
$\displaystyle\pm\frac{i{\pi}^{0}}{2}\hat{G}^{R}_{0}\left(\left({\hat{\Sigma}}^{A}_{0}-{\hat{\Sigma}}^{R}_{0}\right)\hat{G}^{A}_{0}\bm{\mathbf{\nabla_{\pi}}}\left({\hat{G}_{0}}^{A}\right)^{-1}-\bm{\mathbf{\nabla_{\pi}}}\left({\hat{G}_{0}}^{R}\right)^{-1}\hat{G}^{R}_{0}\left({\hat{\Sigma}}^{A}_{0}-{\hat{\Sigma}}^{R}_{0}\right)\right)\hat{G}^{A}_{0}f_{\mp}^{\prime},$
(34c) $\displaystyle\hat{G}^{<}_{{\bm{\mathbf{B}}}_{T}}=$
$\displaystyle\mp\left(\hat{G}^{A}_{{\bm{\mathbf{B}}}_{T}}-\hat{G}^{R}_{{\bm{\mathbf{B}}}_{T}}\right)f_{\mp}+\hat{G}^{R}_{0}\left({\hat{\Sigma}}^{<}_{{\bm{\mathbf{B}}}_{T}}\pm\left({\hat{\Sigma}}^{A}_{{\bm{\mathbf{B}}}_{T}}-{\hat{\Sigma}}^{R}_{{\bm{\mathbf{B}}}_{T}}\right)f_{\mp}\right)\hat{G}^{A}_{0}$
(34d)

with the abbreviations $f_{\mp}\equiv f_{\mp}(\pi^{0})$ and
$f_{\mp}^{\prime}\equiv{\partial}_{{\pi}^{0}}f_{\mp}$. These equations have
solutions of the form

$\displaystyle\hat{G}^{<}_{{\bm{\mathbf{E}}}_{T}}$
$\displaystyle=\hat{G}^{<}_{{\bm{\mathbf{E}}}_{T},\text{I}}f_{\mp}+\hat{G}^{<}_{{\bm{\mathbf{E}}}_{T},\text{II}}f_{\mp}^{\prime},$
(35a) $\displaystyle{\hat{\Sigma}}^{<}_{{\bm{\mathbf{E}}}_{T}}$
$\displaystyle={\hat{\Sigma}}^{<}_{{\bm{\mathbf{E}}}_{T},\text{I}}f_{\mp}+{\hat{\Sigma}}^{<}_{{\bm{\mathbf{E}}}_{T},\text{II}}f_{\mp}^{\prime},$
(35b) $\displaystyle\hat{G}^{<}_{{\bm{\mathbf{B}}}_{T}}$
$\displaystyle=\hat{G}^{<}_{{\bm{\mathbf{B}}}_{T},\text{I}}f_{\mp}+\hat{G}^{<}_{{\bm{\mathbf{B}}}_{T},\text{II}}f_{\mp}^{\prime},$
(35c) $\displaystyle{\hat{\Sigma}}^{<}_{{\bm{\mathbf{B}}}_{T}}$
$\displaystyle={\hat{\Sigma}}^{<}_{{\bm{\mathbf{B}}}_{T},\text{I}}f_{\mp}+{\hat{\Sigma}}^{<}_{{\bm{\mathbf{B}}}_{T},\text{II}}f_{\mp}^{\prime}.$
(35d)

The first and second terms always represent the interband and intraband
contributions, respectively [5]. By inserting the solutions into Eqs.
(34a)-(34d) and comparing the left and right hand sides, we find

$\displaystyle\hat{G}^{<}_{\bm{\mathbf{E}},\text{I}}=$
$\displaystyle\mp\left(\hat{G}^{A}_{{\bm{\mathbf{E}}}_{T}}-\hat{G}^{R}_{{\bm{\mathbf{E}}}_{T}}\right),$
(36a) $\displaystyle\hat{G}^{<}_{{\bm{\mathbf{E}}}_{T},\text{II}}=$
$\displaystyle\
\hat{G}^{R}_{0}{\hat{\Sigma}}^{<}_{{\bm{\mathbf{E}}}_{T},\text{II}}\hat{G}^{A}_{0}\mp\frac{i{\pi}^{0}}{2}\bm{\mathbf{\nabla_{\pi}}}(\hat{G}^{R}_{0}+\hat{G}^{A}_{0})$
$\displaystyle\pm
i{\pi}^{0}\hat{G}^{R}_{0}\bm{\mathbf{\nabla_{\pi}}}({\hat{h}}_{0}(\bm{\mathbf{\pi}})+\frac{1}{2}\left({\hat{\Sigma}}^{R}_{0}+{\hat{\Sigma}}^{A}_{0}\right))\hat{G}^{A}_{0},$
(36b) $\displaystyle{\hat{\Sigma}}^{<}_{\bm{\mathbf{E}},\text{I}}=$
$\displaystyle\mp\left({\hat{\Sigma}}^{A}_{{\bm{\mathbf{E}}}_{T}}-{\hat{\Sigma}}^{R}_{{\bm{\mathbf{E}}}_{T}}\right),$
(36c) $\displaystyle\hat{G}^{<}_{{\bm{\mathbf{B}}}_{T},\text{I}}=$
$\displaystyle\mp\left(\hat{G}^{A}_{{\bm{\mathbf{B}}}_{T}}-\hat{G}^{R}_{{\bm{\mathbf{B}}}_{T}}\right),$
(36d) $\displaystyle\hat{G}^{<}_{{\bm{\mathbf{B}}}_{T},\text{II}}=$
$\displaystyle\ \hat{\Sigma}^{<}_{{\bm{\mathbf{B}}}_{T},\text{II}}=0,$ (36e)
$\displaystyle{\hat{\Sigma}}^{<}_{{\bm{\mathbf{B}}}_{T},\text{I}}=$
$\displaystyle\mp\left({\hat{\Sigma}}^{A}_{{\bm{\mathbf{B}}}_{T}}-{\hat{\Sigma}}^{R}_{{\bm{\mathbf{B}}}_{T}}\right).$
(36f)

Eqs. (34a), (34b), (35a)-(35d), and (36a)-(36f) provide a complete
description.

##  IV Heat Transport

The thermal conductivity tensor is defined as [29]

$\kappa_{ij}=\tilde{\kappa}_{ij}+\frac{2}{T}\epsilon_{ijk}M_{Q}^{k},$ (37)

where $\tilde{\kappa}_{ij}$ is the contribution from heat current density,
$M_{Q}^{k}$ is the $k$th component of $\bm{\mathbf{M}}_{Q}$, and
$\epsilon_{ijk}$ is the Levi-Civita tensor. In the Wigner representation Eq.
(5) becomes [49]

$\displaystyle\bm{\mathbf{J}}_{Q}$ $\displaystyle=\pm
i\text{Tr}\int\frac{d\pi^{0}}{2\pi}(\pi^{0}\bm{\mathbf{\hat{v}}}\mathcal{G})^{<},$
$\displaystyle=\pm
i\hbar\text{Tr}\int\frac{d\pi^{0}}{2\pi}\pi^{0}\bm{\mathbf{\hat{v}}}\left(\hat{G}_{\bm{\mathbf{E}}_{T},\text{I}}^{<}f_{\mp}+\hat{G}_{\bm{\mathbf{E}}_{T},\text{II}}^{<}f^{\prime}_{\mp}\right)\bm{\mathbf{E}}_{T}.$
(38)

Its contribution to ${\kappa}_{ij}$ is

$\displaystyle\tilde{\kappa}_{ij}$
$\displaystyle=\pm\frac{i\hbar}{T}\text{Tr}\int\frac{d\pi^{0}}{2\pi}\pi^{0}\hat{v}_{i}\left(\hat{G}_{\bm{\mathbf{E}}_{T},\text{I}}^{<}f_{\mp}+\hat{G}_{\bm{\mathbf{E}}_{T},\text{II}}^{<}f^{\prime}_{\mp}\right)_{j},$
(39)

where $\hat{v}_{i}$ is the velocity operator in direction $i$, the temperature
gradient points in direction $j$, and $i,j=x,y,z$. Before calculating the
contribution of the heat magnetization to $\kappa_{ij}$ by Eq. (6), it is
common practice [29, 49] to calculate the auxiliary heat magnetization
$\tilde{\bm{\mathbf{M}}}_{Q}=-\left(\partial
K/\partial\bm{\mathbf{B}}_{T}\right)_{T,\mu}$, where $K=E-\mu N$ is the grand
canonical energy. We obtain
${\tilde{\bm{\mathbf{M}}}}_{Q}=\beta^{-1}\partial\left({\beta}^{2}{\bm{\mathbf{M}}}_{Q}\right)/\partial\beta$
and have in the Wigner representation

$\displaystyle K=\pm
i\text{Tr}\int\frac{d\pi^{0}}{2\pi}(\pi^{0}\mathcal{G})^{<},$ (40)

resulting in

$\displaystyle\tilde{\bm{\mathbf{M}}}_{Q}$ $\displaystyle=\pm
i\hbar\text{Tr}\int\frac{d\pi^{0}}{2\pi}\pi^{0}\hat{G}_{\bm{\mathbf{B}}_{T},\text{I}}^{<}f_{\mp}.$
(41)

By integration with respect to $\beta$, we arrive at

$\displaystyle\bm{\mathbf{M}}_{Q}$
$\displaystyle=\pm\frac{i\hbar}{\beta^{2}}\text{Tr}\int_{\infty}^{\beta}d\beta\int\frac{d\pi^{0}}{2\pi}\pi^{0}\hat{G}_{\bm{\mathbf{B}}_{T},\text{I}}^{<}\beta
f_{\mp}$ (42)

and

$\displaystyle\kappa_{ij}=$
$\displaystyle\pm\frac{i\hbar}{T}\text{Tr}\int\frac{d\pi^{0}}{2\pi}\pi^{0}\bigg{(}\hat{v}_{i}\left(\hat{G}_{\bm{\mathbf{E}}_{T},\text{I}}^{<}\right)_{j}f_{\mp}+\hat{v}_{i}\left(\hat{G}_{\bm{\mathbf{E}}_{T},\text{II}}^{<}\right)_{j}f^{\prime}_{\mp}$
$\displaystyle+\epsilon_{ijk}\frac{2}{\beta^{2}}\int_{\infty}^{\beta}d\beta\left(\hat{G}_{\bm{\mathbf{B}}_{T},\text{I}}^{<}\right)_{k}\beta
f_{\mp}\bigg{)}.$ (43)

Together with Eqs. (34) and (36) this equation provides an accurate
description of the thermal transport properties in multiband systems and is
the central result of this work.

Figure 1: (a) Lattice representation of the stacked honeycomb ferromagnet with
the spins oriented along the $z$ direction. Within a layer the nearest and
next-nearest neighbors are separated by the vectors
$\bm{\mathbf{d}}_{\beta}^{(1)}$ and $\bm{\mathbf{d}}_{\beta}^{(2)}$,
respectively. The distance between two honeycomb lattice sites A and B is $a$
and the distance between adjacent layers is $c$. We set $a=c=1$. The
intralayer Heisenberg interaction between nearest neighbors is $J$. The
interlayer Heisenberg interaction between sites A and B in adjacent layers is
$J_{A}$ and $J_{B}$, respectively. (b) The Dzyaloshinskii-Moriya interaction
$D$, represented by black dots, acts perpendicularly to the layers and
involves next-nearest neighbors.

To prove that our approach is consistent with the literature we consider the
clean limit $\hat{\Sigma}=0$. We have
$\hat{G}_{0}^{A(R)}=\left(\pi^{0}-\hat{h}_{0}\mp i\eta\right)^{-1}$, where
$\eta$ is a small positive real number and
$\bm{\mathbf{\nabla_{\bm{\mathbf{p}}}}}\left(\hat{G}_{0}^{{R(A)}}\right)^{-1}=-\hat{\bm{\mathbf{v}}}$.
Using the eigenbasis of $\hat{h}_{0}$,
$\hat{h}_{0}|n\rangle=\epsilon_{n}|n\rangle$, to evaluate the trace and the
residue theorem to perform the $\pi^{0}$ integration we find

$\displaystyle\tilde{\kappa}_{ij}$ $\displaystyle=\frac{1}{\hbar
T}\sum_{n}\left(m_{n}^{k}\epsilon_{n}f_{\mp}(\epsilon_{n})-\Omega_{n}^{k}\epsilon_{n}^{2}f_{\mp}(\epsilon_{n})\right),$
(44) $\displaystyle M_{Q}^{k}$
$\displaystyle=\frac{1}{2\hbar}\sum_{n}\left(2\Omega_{n}^{k}\int_{\infty}^{\epsilon_{n}}d\epsilon\epsilon
f_{\mp}(\epsilon)-m_{n}^{k}\epsilon_{n}f_{\mp}(\epsilon_{n})\right).$ (45)

The magnetic moment reads

$\displaystyle m_{n}^{k}$
$\displaystyle=-2\epsilon_{ijk}\text{Im}\sum_{m}\frac{\langle
n|\hat{v}_{i}|m\rangle\langle
m|\hat{v}_{j}|n\rangle}{(\epsilon_{n}-\epsilon_{m})}$ (46)

and the Berry curvature reads

$\displaystyle\Omega_{n}^{k}$
$\displaystyle=-2\epsilon_{ijk}\text{Im}\sum_{m}\frac{\langle
n|\hat{v}_{i}|m\rangle\langle
m|\hat{v}_{j}|n\rangle}{(\epsilon_{n}-\epsilon_{m})^{2}}.$ (47)

By combining Eqs. (44) and (45) we arrive at

$\displaystyle{\kappa}_{ij}$ $\displaystyle=-\frac{1}{\hbar
T}\sum_{n}\Omega_{n}^{k}\left(\epsilon_{n}^{2}f_{\mp}(\epsilon_{n})-2\int_{\infty}^{\epsilon_{n}}d\epsilon\epsilon
f_{\mp}(\epsilon)\right)$
$\displaystyle=-\frac{1}{\hbar\beta^{2}T}\sum_{n}\Omega_{n}^{k}c_{2}(f_{\mp}),$
(48)

where

$\displaystyle
c_{2}(f_{\mp})=-\int^{\infty}_{\epsilon_{n}}d\epsilon(\beta\epsilon)^{2}f_{\mp}^{\prime}(\epsilon).$
(49)

This is the expression derived in Refs. [33, 34, 35, 36, 37, 29, 49, 40].

Figure 2: Phase diagram of the stacked honeycomb ferromagnet in the
$(K_{-},J_{-})$ plane. The black dots mark the points chosen to analyze the
behavior of the Chern number and thermal conductivity in phases 1 to 5. Figure
3: (a) Band structure and (b) Berry curvature of the lower magnon band in the
$k_{x}k_{z}$ plane at $k_{y}=0$. WM nodes are evident in phases 3 to 5. The
Berry curvature vector is indicated by arrows. The color represents the
divergence of the Berry curvature with blue indicating negative values and red
indicating positive values. Figure 4: Berry curvature of the lower magnon
band: (a) $\Omega_{1}^{x}$ at $k_{x}=0$, (b) $\Omega_{1}^{y}$ at
$k_{y}=\pi/3$, and (c) $\Omega_{1}^{z}$ at $k_{z}=\pi$. Blue color indicates
negative values and red color indicates positive values.

##  V Thermal Hall effect of Weyl Magnons in a Stacked Honeycomb Ferromagnet

###  V.1 Model

The thermal Hall effects in two-dimensional ferromagnetic and
antiferromagnetic magnonic systems (induced by spin-orbit coupling) were
studied extensively [38, 61, 33, 36, 37, 41, 65, 67, 63, 64, 62, 66], while
less attention was given to three-dimensional magnonic systems. In analogy to
electrons, magnons can exhibit non-trivial topological properties. One example
are the Weyl magnons (WMs) in three-dimensional magnets [87, 88, 86]. To study
the impact of WMs on $\kappa_{ij}$, we address the toy model of the stacked
honeycomb ferromagnet shown in Fig. 1, whose topological properties were
determined in Ref. [96]. The Hamiltonian reads

$\displaystyle H=$ $\displaystyle-J\sum_{\langle
i,j\rangle,n}\bm{\mathbf{S}}_{i,n}\cdot\bm{\mathbf{S}}_{j,n}+\sum_{\llangle
i,j\rrangle}\bm{\mathbf{D}}\cdot\left(\bm{\mathbf{S}}_{i,n}\times\bm{\mathbf{S}}_{j,n}$
$\displaystyle-\sum_{i,n}K_{i}\left(S_{i,n}^{z}\right)^{2}-\sum_{i,\langle
n,m\rangle}J_{i}\bm{\mathbf{S}}_{i,n}\cdot\bm{\mathbf{S}}_{i,m},$

where $\langle\cdots\rangle$ and $\llangle\cdots\rrangle$ represent the
nearest neighbor and next-nearest neighbor sites, respectively, and
$\bm{\mathbf{S}}_{i,n}$ represents the spin at site $i$ (A or B) of layer $n$.
The first term is the symmetric Heisenberg interaction with $J>0$ and the
second term is the antisymmetric Dzyaloshinskii-Moriya interaction with
Dzyaloshinskii-Moriya vector $\bm{\mathbf{D}}$. Only the perpendicular
component $D\nu_{ij}$ (with $\nu_{ij}$ alternating between $-1$ and $1$ on
consecutive intralayer bonds) of $\bm{\mathbf{D}}$ survives [58, 59]. The
third term is the anisotropy energy and the last term is the interlayer
interaction. A slightly different model with third-nearest neighbors was
proposed in Ref. [97]. However, this model has only either WM phases or magnon
Chern insulator phases with bands of Chern number either $-2$ or $1$, while
our aim is to compare $\kappa_{ij}$ of different phases, including a trivial
insulator phase with bands of Chern number $0$.

We start our analysis of $H$ with the standard linear spin-wave method [98,
99]. Focusing solely on the lowest order in the magnon creation
($\gamma_{i,n}^{\dagger}$) and annihilation ($\gamma_{i,n}$) operators, we
utilize the Holstein-Primakoff transformation [100]

$\displaystyle\begin{aligned}
S_{i,n}^{+}=&S_{i,n}^{x}+iS_{i,n}^{y}=\sqrt{2S_{i,n}}\sqrt{1-{\frac{\gamma_{i,n}^{\dagger}\gamma_{i,n}}{2S_{i,n}}}}\gamma_{i,n},\\\
S_{i,n}^{-}=&S_{i,n}^{x}-iS_{i,n}^{y}=\sqrt{2S_{i,n}}\gamma_{i,n}^{\dagger}\;\sqrt{1-{\frac{\gamma_{i,n}^{\dagger}\gamma_{i,n}}{2S_{i,n}}}}\;,\\\
S_{i,n}^{z}=&(S_{i,n}-\gamma_{i,n}^{\dagger}\gamma_{i,n}).\end{aligned}$ (50)

Keeping only the quadratic terms, we then perform a Fourier transformation
using the relation

$\gamma_{i,n}^{(\dagger)}=\frac{1}{\sqrt{L}}\sum_{\bm{\mathbf{k}}}\text{e}^{(-)i\bm{\mathbf{r}}_{i,n}\cdot\bm{\mathbf{k}}}\gamma_{\bm{\mathbf{k}}}^{(\dagger)},$
(51)

where $L$ denotes the number of unit cells,
$\gamma_{\bm{\mathbf{k}},n}^{\dagger}$ and $\gamma_{\bm{\mathbf{k}},n}$ denote
the magnon operators in momentum space, and $\bm{\mathbf{r}}_{i,n}$ denotes
the position vector. We obtain

$H=\sum_{\bm{\mathbf{k}}}\Gamma_{\bm{\mathbf{k}}}^{\dagger}h(\bm{\mathbf{k}})\Gamma_{\bm{\mathbf{k}}}$
(52)

with the basis vectors $\Gamma_{\bm{\mathbf{k}}}$ and

$h(\bm{\mathbf{k}})=\sum_{\alpha\;=\;0,x,y,z}h_{\alpha}(\bm{\mathbf{k}})\sigma_{\alpha},$
(53)

where $\sigma_{0}$ is the identity matrix, $\sigma_{x}$, $\sigma_{y}$, and
$\sigma_{z}$ are the Pauli matrices, and

$\displaystyle\begin{split}h_{0}(\bm{\mathbf{k}})&=3JS+K_{+}S+J_{+}S\left(1-\cos(k_{z}c)\right),\\\
h_{x}(\bm{\mathbf{k}})&=-JS\sum_{\beta=1}^{3}\cos(\bm{\mathbf{k}}\cdot\bm{\mathbf{d}}_{\beta}^{(1)}),\\\
h_{y}(\bm{\mathbf{k}})&=-JS\sum_{\beta=1}^{3}\sin(\bm{\mathbf{k}}\cdot\bm{\mathbf{d}}_{\beta}^{(1)}),\\\
h_{x}(\bm{\mathbf{k}})&=2DS\sum_{\beta=1}^{3}\sin(\bm{\mathbf{k}}\cdot\bm{\mathbf{d}}_{\beta}^{(2)})+K_{-}S+J_{-}S\left(1-\cos(k_{z}c)\right)\end{split}$
(54)

with $K_{\pm}=K_{\text{A}}\pm K_{\text{B}}$ and $J_{\pm}=J_{\text{A}}\pm
J_{\text{B}}$. The vectors $\bm{\mathbf{d}}_{\beta}^{(1)}$ between nearest
neighbor sites and $\bm{\mathbf{d}}_{\beta}^{(2)}$ between next-nearest
neighbor sites are shown in Fig. 1.

Figure 5: Berry curvature ($z$ component) of the lower magnon band for (a)
phase 2 and (b) phase 3 at different values of $k_{z}$. The maxima transfer
from the K point (white square) to the K’ point (black square) as $k_{z}$
changes from $0$ to $\pi$. In phase 3 the maxima become minima abruptly as
$k_{z}$ passes through the WM node at
$k_{z}=\cos^{-1}(K_{-}/J_{-}+1-3\sqrt{3}D/J_{-})\approx 0.2677\pi$. Blue color
indicates negative values and red color indicates positive values.

The phase diagram in the $(K_{-},J_{-})$ plane is constructed for
$D=0.2J,K_{+}=12D$, and $J_{+}=2D$ (following Ref. [96]). The phase boundaries
are obtained by evaluating
$h_{x}(\bm{\mathbf{k}})=h_{y}(\bm{\mathbf{k}})=h_{z}(\bm{\mathbf{k}})=0$ to
determine the parameters at which bands cross each other. This gives rise to
four boundaries, $K_{-}=\pm 3\sqrt{3}D$ and $K_{-}=-2J_{-}\pm 3\sqrt{3}D$, and
generates one gapped trivial insulator phase (phase 1), one gapped magnon
Chern insulator phase (phase 2), and three gapless WM phases (phases 3 to 5),
as shown in Fig. 2. Phase 3 has a pair of WM nodes at
$\left(-4\pi/3\sqrt{3},0,\pm\cos^{-1}(K_{-}/J_{-}+1-3\sqrt{3}D/J_{-})\right)$,
phase 4 has a pair of WM nodes at
$\left(-4\pi/3\sqrt{3},0,\pm\cos^{-1}(K_{-}/J_{-}+1+3\sqrt{3}D/J_{-})\right)$,
and phase 5 combines the two pairs.

To analyze the phases, we consider the points $(K_{-},J_{-})$ marked in Fig.
2: $(4\sqrt{3},3\sqrt{3})D$ for phase 1, $(\sqrt{3},-\sqrt{3})D$ for phase 2,
$(2\sqrt{3},3\sqrt{3})D$ for phase 3, $(-2\sqrt{3},-3\sqrt{3})D$ for phase 4,
and $(-4\sqrt{3},4\sqrt{3})D$ for phase 5. The band structure and Berry
curvature of the lower magnon band in the $k_{x}k_{z}$ plane at $k_{y}=0$ are
addressed in Fig. 3. Monopoles of the Berry curvature are present at the WM
nodes in phases 3 to 5.

Figure 6: Two-dimensional thermal conductivity ($xy$ component) for (a) phase
1 (trivial insulator), (b) phase 2 (magnon Chern insulator), (c) WM phases 3
and 4, and (d) WM phase 5. The Chern number of the lower magnon band is
indicated by white ($\text{C}_{1}$ = 0) and gray ($\text{C}_{1}$ = 1) colors.

###  V.2 Chern number and thermal conductivity

We decompose the three-dimensional Hamiltonian into two-dimensional
Hamiltonians indexed by $k_{z}$,

$h(\bm{\mathbf{k}})=\sum_{k_{z}}h_{k_{x},k_{y}}(k_{z}),$ (55)

and calculate the Chern number of the $n$th magnon band ($n=1$ or $2$) in the
traditional way as

$\text{C}_{n}(k_{z})=\frac{1}{2\pi}\int dk_{y}\int
dk_{x}\Omega_{n}^{z}(k_{z}).$ (56)

Similarly, Eq. (IV) implies for the three-dimensional thermal conductivity

$\kappa_{ij}^{3D}=\int\frac{dk_{z}}{2\pi}\kappa_{ij}^{2D}(k_{z}),$ (57)

where

$\kappa_{ij}^{2D}(k_{z})=-\frac{1}{\hbar\beta^{2}T}\sum_{n}\int dk_{y}\int
dk_{x}\Omega_{n}^{k}(k_{z})c_{2}(f_{+})$ (58)

is the two-dimensional thermal conductivity.

Figure 7: (a) Chern number of the lower magnon band and (b) two-dimensional
thermal conductivity ($xy$ component) at $T$ = 100 K in the ($K_{-},J_{-}$)
plane for different values of $k_{z}$.

WM nodes appear in pairs with opposite chirality for planes not orthogonal to
the $k_{z}$ axis, resulting in zero net Berry curvature, as illustrated in
Fig. 4(a,b) for the $k_{y}k_{z}$ and $k_{x}k_{z}$ planes. As a consequence,
$\kappa_{yz}^{3D}$ and $\kappa_{zx}^{3D}$ are zero. On the other hand, for
planes perpendicular to the $k_{z}$ axis, see Fig. 4(c) for the $k_{x}k_{y}$
plane, the net Berry curvature is non-zero, leading to non-zero
$\kappa_{xy}^{3D}$. Focusing on the $k_{x}k_{y}$ plane, we examine how the $z$
component of the Berry curvature, $\Omega_{n}^{z}$, behaves at different
values of $k_{z}$ to understand the contribution of the WM nodes to
$\kappa_{xy}^{3D}$. Note that in the WM phases most contributions are due to
the lower magnon band ($n=1$) because of the Bose distribution. It turns out
that $\Omega_{1}^{z}$ is minimal at $k_{z}=\pm\pi$ and maximal at $k_{z}=0$ in
phase 1. It is an even function of $k_{z}$ in phase 2 with maxima that are
roughly ten times higher than in phase 1. Moreover, in phase 2 the maxima
transfer from the K point to the K’ point as $k_{z}$ changes from $0$ to
$\pi$, see Fig. 5(a). In all the WM phases the maxima of $\Omega_{1}^{z}$
become minima abruptly as $k_{z}$ passes through a WM node, see Fig. 5(b) for
the WM node at $k_{z}=\cos^{-1}(K_{-}/J_{-}+1-3\sqrt{3}D/J_{-})\approx
0.2677\pi$ in phase 3. This results in a smaller $\Omega_{1}^{z}$ than that of
phase 2, leading to a smaller $\kappa_{xy}^{3D}$.

Figure 6 presents $\text{C}_{\text{1}}$ and $\kappa_{xy}^{2D}$ as functions of
$k_{z}$ for all five phases. Phase 1 (trivial insulator) has
$\text{C}_{\text{1}}=0$ and phase 2 (magnon Chern insulator) has
$\text{C}_{\text{1}}=1$. $\text{C}_{\text{1}}$ of the WM phases changes
between 0 and 1 as $k_{z}$ passes through WM nodes. $\kappa_{xy}^{2D}$ remains
almost constant in phases 1 and 2, being significantly larger for the magnon
Chern insulator. While it is influenced significantly by $\text{C}_{\text{1}}$
in the WM phases. Phases 3 and 4 show the same values because of identical
$|K_{-}|$ and $|J_{-}|$. Figure 7 presents $\text{C}_{\text{1}}$ and
$\kappa_{xy}^{2D}$ at $T$ = 100 K in the ($K_{-},J_{-}$) plane.
$\kappa_{xy}^{2D}$ transitions a gradually near the topological boundaries,
which can be attributed to the emergence of chiral edge states, a
characteristic feature of two-dimensional magnon Chern insulators. As $T$
increases, resulting in a higher magnon density, the transition at the
topological boundaries becomes more pronounced.

Figure 8 presents $\kappa_{xy}^{3D}$ as a function of $T$ and at $T$ = 100 K
in the ($K_{-},J_{-}$) plane. The absolute value increases for increasing $T$
in all the phases. The highest absolute values are observed for the magnon
Chern insulator, which is not surprising, because in this phase, unlike the WM
phases, $\Omega_{1}^{z}$ maintains its sign throughout the $k_{x}k_{y}$ plane,
as shown in Fig. 5(a).

##  VI Summary and Discussion

We presented a theoretical framework for investigating the response of heat
currents to a temperature gradient using the Keldysh-Dyson formalism. The core
of this approach is an expansion of the Dyson equations through the Moyal
product within the Wigner representation. The temperature gradient is
described using both Luttinger’s scalar potential ($\phi$) and Moreno-Coleman-
Tatara’s vector potential ($\bm{\mathbf{A}}_{T}$). While the thermal field
does not arise from a formal gauge symmetry, we can achieve “gauge invariance”
under the transformations $\phi\to\phi+\frac{d\chi}{dt}$ and
$\bm{\mathbf{A}}_{T}\to\bm{\mathbf{A}}_{T}-\bm{\mathbf{\nabla}}\chi$ by
assigning a portion of the temperature gradient to $\phi$ and the remainder to
$\bm{\mathbf{A}}_{T}$. This establishes a one-to-one correspondence between
the responses induced by an electromagnetic field and by a temperature
gradient.

We derived the thermal conductivity tensor encompassing contributions from the
heat current and heat magnetization, applicable to general multiband
Hamiltonians. Our formalism treats the heat current and heat magnetization on
equal footing, mirroring the usage of the method of Ref. [5] for calculating
the electric current density and subsequently calculating the orbital
magnetization [101]. Our result for clean limit reduces to the well
established expression [33, 34, 35, 36, 37, 40, 49, 29].

Figure 8: Three-dimensional thermal conductivity ($xy$ component) (a) as a
function of $T$ and (b) at $T$ = 100 K in the ($K_{-},J_{-}$) plane.

The introduced approach has several advantages over previous formalisms.
Firstly, it readily accounts for diamagnetic currents by replacing the
canonical momentum with the kinetic momentum. Secondly, it can be applied to
interacting and/or disordered systems. Thirdly, its ability to address non-
linear responses extends beyond the scope of formalisms based on the Kubo
formula. Fourthly, it can be extended to other transport quantities (e.g.,
particle transport, nonequilibrium densities) and time-dependent non-uniform
problems.

Using the derived formula for the thermal conductivity in the clean limit, the
thermal properties of a stacked honeycomb ferromagnet were analyzed for all
its phases (trivial insulator, magnon Chern insulator, and three WM phases
with WM nodes at different points in the momentum space). Decomposing the
three-dimensional Hamiltonian into two-dimensional Hamiltonians indexed by
$k_{z}$, we observed that the Chern number remains constant for all values of
$k_{z}$ for the trivial and magnon Chern insulator phases, while it changes
abruptly from 0 to 1 as $k_{z}$ passes through WM nodes. Accordingly, in the
trivial and magnon Chern insulator phases the $xy$-component of the two-
dimensional thermal conductivity exhibits only minor variations as a function
of $k_{z}$, while it changes significantly as $k_{z}$ passes through WM nodes
in the WM phases. The primary contributions to the three-dimensional thermal
conductivity in the WM phases are due to the WM nodes (monopoles of the Berry
curvature). While one might expect the absolute value of the $xy$-component of
the three-dimensional thermal conductivity to be higher in the WM phases than
in the other two phases, interestingly, the magnon Chern insulator exhibits
the highest absolute values, which can be understood by noting that in the WM
phases the Berry curvature changes sign as $k_{z}$ passes through WM nodes.

###### Acknowledgements.

The authors acknowledge insightful discussions with Gen Tatara, Konstantinos
Sourounis, and Diego Garcia Ovalle. This work was supported by King Abdullah
University of Science and Technology.

## References

  * Kubo [1957] R. Kubo, Statistical-mechanical theory of irreversible processes. I. General theory and simple applications to magnetic and conduction problems, J. Phys. Soc. Jpn. 12, 570 (1957).
  * Kubo  _et al._ [1957] R. Kubo, M. Yokota, and S. Nakajima, Statistical-mechanical theory of irreversible processes. II. Response to thermal disturbance, J. Phys. Soc. Jpn. 12, 1203 (1957).
  * Smrcka and Středa [1977] L. Smrcka and P. Středa, Transport coefficients in strong magnetic fields, J. Phys. C: Solid State Phys. 10, 2153 (1977).
  * Rammer and Smith [1986] J. Rammer and H. Smith, Quantum field-theoretical methods in transport theory of metals, Rev. Mod. Phys. 58, 323 (1986).
  * Onoda  _et al._ [2006] S. Onoda, N. Sugimoto, and N. Nagaosa, Theory of non-equilibirum states driven by constant electromagnetic fields, Prog. Theor. Phys. 116, 61 (2006).
  * Peierls [1933] R. Peierls, On the theory of diamagnetism of conduction electrons, Z. Phys. 80, 763 (1933).
  * Xiao  _et al._ [2010] D. Xiao, M.-C. Chang, and Q. Niu, Berry phase effects on electronic properties, Rev. Mod. Phys. 82, 1959 (2010).
  * Nagaosa  _et al._ [2010] N. Nagaosa, J. Sinova, S. Onoda, A. H. MacDonald, and N. P. Ong, Anomalous Hall effect, Rev. Mod. Phys. 82, 1539 (2010).
  * Sinova  _et al._ [2015] J. Sinova, S. O. Valenzuela, J. Wunderlich, C. H. Back, and T. Jungwirth, Spin Hall effects, Rev. Mod. Phys. 87, 1213 (2015).
  * Chang  _et al._ [2023] C.-Z. Chang, C.-X. Liu, and A. H. MacDonald, Colloquium: Quantum anomalous Hall effect, Rev. Mod. Phys. 95, 011002 (2023).
  * Luttinger [1964] J. M. Luttinger, Theory of thermal transport coefficients, Phys. Rev. 135, A1505 (1964).
  * Chernodub  _et al._ [2022] M. N. Chernodub, Y. Ferreiros, A. G. Grushin, K. Landsteiner, and M. A. Vozmediano, Thermal transport, geometry, and anomalies, Phys. Rep. 977, 1 (2022).
  * Bauer  _et al._ [2012] G. E. W. Bauer, E. Saitoh, and B. J. van Wees, Spin caloritronics, Nat. Mater. 11, 391 (2012).
  * Boona  _et al._ [2014] S. R. Boona, R. C. Myers, and J. P. Heremans, Spin caloritronics, Energy Environ. Sci. 7, 885 (2014).
  * Yokoyama and Murakami [2014] T. Yokoyama and S. Murakami, Spintronics and spincaloritronics in topological insulators, Physica E Low Dimens. Syst. Nanostruct. 55, 1 (2014).
  * Yu  _et al._ [2017] H. Yu, S. D. Brechet, and J.-P. Ansermet, Spin caloritronics, origin and outlook, Phys. Lett. A 381, 825 (2017).
  * Uchida [2021] K.-I. Uchida, Transport phenomena in spin caloritronics, Proc. Jpn. Acad., Ser. B, Phys. Biol. Sci. 97, 69 (2021).
  * Rückriegel  _et al._ [2018] A. Rückriegel, A. Brataas, and R. A. Duine, Bulk and edge spin transport in topological magnon insulators, Phys. Rev. B 97, 081106 (2018).
  * Aguilera  _et al._ [2020] E. Aguilera, R. Jaeschke-Ubiergo, N. Vidal-Silva, L. E. F. F. Torres, and A. S. Nunez, Topological magnonics in the two-dimensional van der Waals magnet Cr$\text{I}_{3}$, Phys. Rev. B 102, 024409 (2020).
  * Barman  _et al._ [2021] A. Barman, G. Gubbiotti, S. Ladak, A. O. Adeyeye, M. Krawczyk, J. Gräfe, C. Adelmann, S. Cotofana, A. Naeemi, V. I. Vasyuchka, B. Hillebrands, S. A. Nikitov, H. Yu, D. Grundler, A. V. Sadovnikov, A. A. Grachev, S. E. Sheshukova, J.-Y. Duquesne, M. Marangolo, G. Csaba, W. Porod, V. E. Demidov, S. Urazhdin, S. O. Demokritov, E. Albisetti, D. Petti, R. Bertacco, H. Schultheiss, V. V. Kruglyak, V. D. Poimanov, S. Sahoo, J. Sinha, H. Yang, M. Münzenberg, T. Moriyama, S. Mizukami, P. Landeros, R. A. Gallardo, G. Carlotti, J.-V. Kim, R. L. Stamps, R. E. Camley, B. Rana, Y. Otani, W. Yu, T. Yu, G. E. W. Bauer, C. Back, G. S. Uhrig, O. V. Dobrovolskiy, B. Budinska, H. Qin, S. van Dijken, A. V. Chumak, A. Khitun, D. E. Nikonov, I. A. Young, B. W. Zingsem, and M. Winklhofer, The 2021 magnonics roadmap, J. Phys.: Condens. Matter 33, 413001 (2021).
  * Wang and Wang [2021] X. S. Wang and X. R. Wang, Topological magnonics, J. Appl. Phys. 129, 151101 (2021).
  * McClarty [2022] P. A. McClarty, Topological magnons: A review, Annu. Rev. Condens. Matter Phys. 13, 171 (2022).
  * Zhuo  _et al._ [2023] F. Zhuo, J. Kang, A. Manchon, and Z. Cheng, Topological phases in magnonics, Adv. Phys. Res. in press, 2300054 (2023).
  * Pekola and Karimi [2021] J. P. Pekola and B. Karimi, Colloquium: Quantum heat transport in condensed matter systems, Rev. Mod. Phys. 93, 041001 (2021).
  * Ridley  _et al._ [2022] M. Ridley, N. W. Talarico, D. Karlsson, N. L. Gullo, and R. Tuovinen, A many-body approach to transport in quantum systems: From the transient regime to the stationary state, J. Phys. A Math. Theor. 55, 273001 (2022).
  * Rivas and Martin-Delgado [2017] Á. Rivas and M. A. Martin-Delgado, Topological heat transport and symmetry-protected boson currents, Sci. Rep. 7, 6350 (2017).
  * Oji and Streda [1985] H. Oji and P. Streda, Theory of electronic thermal transport: Magnetoquantum corrections to the thermal transport coefficients, Phys. Rev. B 31, 7291 (1985).
  * Michaeli and Finkel’stein [2009a] K. Michaeli and A. M. Finkel’stein, Quantum kinetic approach for studying thermal transport in the presence of electron-electron interactions and disorder, Phys. Rev. B 80, 115111 (2009a).
  * Qin  _et al._ [2011] T. Qin, Q. Niu, and J. Shi, Energy magnetization and the thermal Hall effect, Phys. Rev. Lett. 107, 236601 (2011).
  * Schwiete and Finkel’stein [2014] G. Schwiete and A. M. Finkel’stein, Thermal transport and Wiedemann-Franz law in the disordered Fermi liquid, Phys. Rev. B 90, 060201 (2014).
  * Eich  _et al._ [2014] F. G. Eich, A. Principi, M. Di Ventra, and G. Vignale, Luttinger-field approach to thermoelectric transport in nanoscale conductors, Phys. Rev. B 90, 115116 (2014).
  * Kapustin and Spodyneiko [2021] A. Kapustin and L. Spodyneiko, Microscopic formulas for thermoelectric transport coefficients in lattice systems, Phys. Rev. B 104, 035150 (2021).
  * Matsumoto and Murakami [2011a] R. Matsumoto and S. Murakami, Theoretical prediction of a rotating magnon wave packet in ferromagnets, Phys. Rev. Lett. 106, 197202 (2011a).
  * Matsumoto and Murakami [2011b] R. Matsumoto and S. Murakami, Berry phase and thermal transport coefficients in magnon systems, J. Phys.: Conf. Ser. 302, 012025 (2011b).
  * Murakami [2011] S. Murakami, Two-dimensional topological insulators and their edge states, J. Phys.: Conf. Ser. 302, 012019 (2011).
  * Matsumoto and Murakami [2011c] R. Matsumoto and S. Murakami, Rotational motion of magnons and the thermal Hall effect, Phys. Rev. B 84, 184406 (2011c).
  * Matsumoto  _et al._ [2014] R. Matsumoto, R. Shindou, and S. Murakami, Thermal Hall effect of magnons in magnets with dipolar interaction, Phys. Rev. B 89, 054420 (2014).
  * Mook  _et al._ [2014a] A. Mook, J. Henk, and I. Mertig, Magnon Hall effect and topology in kagome lattices: A theoretical investigation, Phys. Rev. B 89, 134409 (2014a).
  * Mook  _et al._ [2014b] A. Mook, J. Henk, and I. Mertig, Edge states in topological magnon insulators, Phys. Rev. B 90, 024412 (2014b).
  * Murakami and Okamoto [2017] S. Murakami and A. Okamoto, Thermal Hall effect of magnons, J. Phys. Soc. Jpn. 86, 011010 (2017).
  * Mook  _et al._ [2019a] A. Mook, J. Henk, and I. Mertig, Thermal Hall effect in noncollinear coplanar insulating antiferromagnets, Phys. Rev. B 99, 014427 (2019a).
  * Mook  _et al._ [2019b] A. Mook, R. R. Neumann, J. Henk, and I. Mertig, Spin Seebeck and spin Nernst effects of magnons in noncollinear antiferromagnetic insulators, Phys. Rev. B 100, 100401 (2019b).
  * Li  _et al._ [2020a] B. Li, A. Mook, A. Raeliarijaona, and A. A. Kovalev, Magnonic analog of the Edelstein effect in antiferromagnetic insulators, Phys. Rev. B 101, 024427 (2020a).
  * Mook  _et al._ [2020] A. Mook, R. R. Neumann, A. Johansson, J. Henk, and I. Mertig, Origin of the magnetic spin Hall effect: Spin current vorticity in the Fermi sea, Phys. Rev. Res. 2, 023065 (2020).
  * Neumann  _et al._ [2022] R. R. Neumann, A. Mook, J. Henk, and I. Mertig, Thermal Hall effect of magnons in collinear antiferromagnetic insulators: Signatures of magnetic and topological phase transitions, Phys. Rev. Lett. 128, 117201 (2022).
  * Qin  _et al._ [2012] T. Qin, J. Zhou, and J. Shi, Berry curvature and the phonon Hall effect, Phys. Rev. B 86, 104305 (2012).
  * Bhalla and Das [2021] P. Bhalla and N. Das, Optical phonon contribution to the thermal conductivity of a quantum paraelectric, J. Phys.: Condens. Matter 33, 345401 (2021).
  * Agarwalla  _et al._ [2011] B. K. Agarwalla, L. Zhang, J.-S. Wang, and B. Li, Phonon Hall effect in ionic crystals in the presence of static magnetic field, Eur. Phys. J. B 81, 197 (2011).
  * Shitade [2014] A. Shitade, Heat transport as torsional responses and Keldysh formalism in a curved spacetime, Prog. Theor. Exp. Phys. 2014, 123I01 (2014).
  * Moreno and Coleman [1996] J. Moreno and P. Coleman, Thermal currents in highly correlated systems, arXiv:cond-mat/9603079v2 (1996).
  * Tatara [2015a] G. Tatara, Thermal vector potential theory of transport induced by a temperature gradient, Phys. Rev. Lett. 114, 196601 (2015a).
  * Tatara [2015b] G. Tatara, Thermal vector potential theory of magnon-driven magnetization dynamics, Phys. Rev. B 92, 064405 (2015b).
  * Gusynin  _et al._ [2014] V. P. Gusynin, S. G. Sharapov, and A. A. Varlamov, Anomalous thermospin effect in the low-buckled Dirac materials, Phys. Rev. B 90, 155107 (2014).
  * Gusynin  _et al._ [2015] V. P. Gusynin, S. G. Sharapov, and A. A. Varlamov, Spin Nernst effect and intrinsic magnetization in two-dimensional Dirac materials, Low Temp. Phys. 41, 342 (2015).
  * Cooper  _et al._ [1997] N. R. Cooper, B. I. Halperin, and I. M. Ruzin, Thermoelectric response of an interacting two-dimensional electron gas in a quantizing magnetic field, Phys. Rev. B 55, 2344 (1997).
  * Michaeli and Finkel’stein [2009b] K. Michaeli and A. M. Finkel’stein, Quantum kinetic approach to the calculation of the Nernst effect, Phys. Rev. B 80, 214516 (2009b).
  * Onose  _et al._ [2010] Y. Onose, T. Ideue, H. Katsura, Y. Shiomi, N. Nagaosa, and Y. Tokura, Observation of the magnon Hall effect, Science 329, 297 (2010).
  * Dzyaloshinsky [1958] I. Dzyaloshinsky, A thermodynamic theory of “weak” ferromagnetism of antiferromagnetics, J. Phys. Chem. Solid. 4, 241 (1958).
  * Moriya [1960] T. Moriya, Anisotropic superexchange interaction and weak ferromagnetism, Phys. Rev. 120, 91 (1960).
  * Okuma [2017] N. Okuma, Magnon spin-momentum locking: Various spin vortices and Dirac magnons in noncollinear antiferromagnets, Phys. Rev. Lett. 119, 107205 (2017).
  * Katsura  _et al._ [2010] H. Katsura, N. Nagaosa, and P. A. Lee, Theory of the thermal Hall effect in quantum magnets, Phys. Rev. Lett. 104, 066403 (2010).
  * Mook  _et al._ [2016a] A. Mook, J. Henk, and I. Mertig, Spin dynamics simulations of topological magnon insulators: From transverse current correlation functions to the family of magnon Hall effects, Phys. Rev. B 94, 174444 (2016a).
  * Owerre [2017a] S. A. Owerre, Topological thermal Hall effect in frustrated kagome antiferromagnets, Phys. Rev. B 95, 014422 (2017a).
  * Doki  _et al._ [2018] H. Doki, M. Akazawa, H.-Y. Lee, J. H. Han, K. Sugii, M. Shimozawa, N. Kawashima, M. Oda, H. Yoshida, and M. Yamashita, Spin thermal Hall conductivity of a kagome antiferromagnet, Phys. Rev. Lett. 121, 097203 (2018).
  * Laurell and Fiete [2018] P. Laurell and G. A. Fiete, Magnon thermal Hall effect in kagome antiferromagnets with Dzyaloshinskii-Moriya interactions, Phys. Rev. B 98, 094419 (2018).
  * Kim  _et al._ [2019] S. K. Kim, K. Nakata, D. Loss, and Y. Tserkovnyak, Tunable magnonic thermal Hall effect in skyrmion crystal phases of ferrimagnets, Phys. Rev. Lett. 122, 057204 (2019).
  * Lu  _et al._ [2019] Y. Lu, X. Guo, V. Koval, and C. Jia, Topological thermal Hall effect driven by spin-chirality fluctuations in frustrated antiferromagnets, Phys. Rev. B 99, 054409 (2019).
  * Kovalev and Zyuzin [2016] A. A. Kovalev and V. Zyuzin, Spin torque and Nernst effects in Dzyaloshinskii-Moriya ferromagnets, Phys. Rev. B 93, 161106 (2016).
  * Zyuzin and Kovalev [2016] V. A. Zyuzin and A. A. Kovalev, Magnon spin Nernst effect in antiferromagnets, Phys. Rev. Lett. 117, 217203 (2016).
  * Cheng  _et al._ [2016] R. Cheng, S. Okamoto, and D. Xiao, Spin Nernst effect of magnons in collinear antiferromagnets, Phys. Rev. Lett. 117, 217202 (2016).
  * Shiomi  _et al._ [2017] Y. Shiomi, R. Takashima, and E. Saitoh, Experimental evidence consistent with a magnon Nernst effect in the antiferromagnetic insulator MnP$\text{S}_{3}$, Phys. Rev. B 96, 134425 (2017).
  * Zhang  _et al._ [2018] Y. Zhang, S. Okamoto, and D. Xiao, Spin-Nernst effect in the paramagnetic regime of an antiferromagnetic insulator, Phys. Rev. B 98, 035424 (2018).
  * Zyuzin and Kovalev [2018] V. A. Zyuzin and A. A. Kovalev, Spin Hall and Nernst effects of Weyl magnons, Phys. Rev. B 97, 174407 (2018).
  * Li  _et al._ [2020b] B. Li, S. Sandhoefner, and A. A. Kovalev, Intrinsic spin Nernst effect of magnons in a noncollinear antiferromagnet, Phys. Rev. Res. 2, 013079 (2020b).
  * Manchon  _et al._ [2014] A. Manchon, P. B. Ndiaye, J.-H. Moon, H.-W. Lee, and K.-J. Lee, Magnon-mediated Dzyaloshinskii-Moriya torque in homogeneous ferromagnets, Phys. Rev. B 90, 224403 (2014).
  * Kovalev  _et al._ [2017] A. A. Kovalev, V. A. Zyuzin, and B. Li, Pumping of magnons in a Dzyaloshinskii-Moriya ferromagnet, Phys. Rev. B 95, 165106 (2017).
  * Zhang  _et al._ [2013] L. Zhang, J. Ren, J.-S. Wang, and B. Li, Topological magnon insulator in insulating ferromagnet, Phys. Rev. B 87, 144101 (2013).
  * Nakata  _et al._ [2017] K. Nakata, S. K. Kim, J. Klinovaja, and D. Loss, Magnonic topological insulators in antiferromagnets, Phys. Rev. B 96, 224414 (2017).
  * McClarty  _et al._ [2018] P. A. McClarty, X.-Y. Dong, M. Gohlke, J. G. Rau, F. Pollmann, R. Moessner, and K. Penc, Topological magnons in Kitaev magnets at high fields, Phys. Rev. B 98, 060404 (2018).
  * Daniels  _et al._ [2019] M. W. Daniels, W. Yu, R. Cheng, J. Xiao, and D. Xiao, Topological spin Hall effects and tunable skyrmion Hall effects in uniaxial antiferromagnetic insulators, Phys. Rev. B 99, 224433 (2019).
  * Díaz  _et al._ [2019] S. A. Díaz, J. Klinovaja, and D. Loss, Topological magnons and edge states in antiferromagnetic skyrmion crystals, Phys. Rev. Lett. 122, 187203 (2019).
  * Wang  _et al._ [2020] X. S. Wang, A. Brataas, and R. E. Troncoso, Bosonic bott index and disorder-induced topological transitions of magnons, Phys. Rev. Lett. 125, 217202 (2020).
  * Lu  _et al._ [2021] Y.-S. Lu, J.-L. Li, and C.-T. Wu, Topological phase transitions of Dirac magnons in honeycomb ferromagnets, Phys. Rev. Lett. 127, 217202 (2021).
  * Owerre [2017b] S. A. Owerre, Floquet topological magnons, J. Phys. Commun. 1, 021002 (2017b).
  * Owerre [2019] S. A. Owerre, Photo-induced Floquet Weyl magnons in noncollinear antiferromagnets, Ann. Phys. 406, 14 (2019).
  * Mook  _et al._ [2016b] A. Mook, J. Henk, and I. Mertig, Tunable magnon Weyl points in ferromagnetic pyrochlores, Phys. Rev. Lett. 117, 157204 (2016b).
  * Owerre [2018a] S. A. Owerre, Weyl magnons in noncoplanar stacked kagome antiferromagnets, Phys. Rev. B 97, 094412 (2018a).
  * Owerre [2018b] S. A. Owerre, Topological thermal Hall effect due to Weyl magnons, Can. J. Phys. 96, 1216 (2018b).
  * Mook  _et al._ [2021] A. Mook, K. Plekhanov, J. Klinovaja, and D. Loss, Interaction-stabilized topological magnon insulator in ferromagnets, Phys. Rev. X 11, 021061 (2021).
  * Gohlke  _et al._ [2023] M. Gohlke, A. Corticelli, R. Moessner, P. A. McClarty, and A. Mook, Spurious symmetry enhancement in linear spin wave theory and interaction-induced topology in magnons, Phys. Rev. Lett. 131, 186702 (2023).
  * Baym and Kadanoff [1961] G. Baym and L. P. Kadanoff, Conservation laws and correlation functions, Phys. Rev. 124, 287 (1961).
  * Dyson [1949] F. J. Dyson, The $S$ matrix in quantum electrodynamics, Phys. Rev. 75, 1736 (1949).
  * Keldysh [1965] L. V. Keldysh, Diagram technique for nonequilibrium processes, Sov. Phys. JETP 20, 1018 (1965).
  * Zhang  _et al._ [2020] Y. Zhang, Y. Gao, and D. Xiao, Thermodynamics of energy magnetization, Phys. Rev. B 102, 235161 (2020).
  * [95] H. Haug and A.-P. Jauho, Boltzmann limit, in: Quantum Kinetics in Transport and Optics of Semiconductors, Springer, Berlin (2008).
  * Su and Wang [2017] Y. Su and X. R. Wang, Chiral anomaly of Weyl magnons in stacked honeycomb ferromagnets, Phys. Rev. B 96, 104437 (2017).
  * Li and Nevidomskyy [2021] S. Li and A. H. Nevidomskyy, Topological Weyl magnons and thermal Hall effect in layered honeycomb ferromagnets, Phys. Rev. B 104, 104419 (2021).
  * Petit [2011] S. Petit, Numerical simulations and magnetism, Collection SFN 12, 105 (2011).
  * Toth and Lake [2015] S. Toth and B. Lake, Linear spin wave theory for single-Q incommensurate magnetic structures, J. Phys.: Condens. Matter 27, 166002 (2015).
  * Holstein and Primakoff [1940] T. Holstein and H. Primakoff, Field dependence of the intrinsic domain magnetization of a ferromagnet, Phys. Rev. 58, 1098 (1940).
  * Zhu  _et al._ [2012] G. Zhu, S. A. Yang, C. Fang, W. M. Liu, and Y. Yao, Theory of orbital magnetization in disordered systems, Phys. Rev. B 86, 214415 (2012).