# Energy magnetization and transport in systems with a non-zero Berry
curvature in a magnetic field

Archisman Panigrahi  Dept. of Physics, Indian Institute of Science, Bangalore
560012, India  Subroto Mukerjee  Dept. of Physics, Indian Institute of
Science, Bangalore 560012, India

###### Abstract

We demonstrate that the well-known expression for the charge magnetization of
a sample with a non-zero Berry curvature can be obtained by demanding that the
Einstein relation holds for the electric transport current. We extend this
formalism to the transport energy current and show that the energy
magnetization must satisfy a particular condition. We provide a physical
interpretation of this condition and relate the energy magnetization to
circulating energy currents due to chiral edge states. We further obtain an
expression for the energy magnetization analogous to the one previously
obtained for the charge magnetization. We also solve the Boltzmann Transport
Equation for the non-equilibrium distribution function in 2D for systems with
a non-zero Berry curvature in a magnetic field. This distribution function can
be used to obtain the regular Hall response in time-reversal invariant samples
with a non-zero Berry curvature, for which there is no anomalous Hall
response.

##  I Introduction

The effect of the Berry curvature on thermoelectric transport in crystals has
received a lot of attention in recent years. Following the seminal work of
Berry on the adiabatic evolution of quantum states [1], it was realized that
the geometric phase identified by him has important consequences for the
semiclassical dynamics of electron wavepackets in crystals. In particular, the
Berry curvature associated with this phase contributes an anomalous term to
the velocity of the wavepacket [2, 3, 4, 5, 6, 7]. This anomalous velocity can
give rise to a Hall response without an external magnetic field, resulting in
the anomalous Hall effect [8], the Magnus Hall Effect [9, 10], as well as many
other interesting electronic transport phenomena [11]. Similarly, in response
to a temperature gradient, without any external magnetic field, it generate a
transverse Hall voltage, known as the anomalous Nernst effect [12, 9, 13, 14].
Coupled with the Boltzmann transport theory, the modified semiclassical
equations have been employed to study transport in topological insulators
[15], Chern Insulators [16], Weyl Semi-Metals [17, 18, 19, 20, 21, 13, 14,
22], Kondo Insulators [23], Rashba systems [24, 25], optical lattices and
quasicrystals [26, 27], superconductors [28], non-Hermitian systems [29, 30],
as well as in various other systems [31, 32, 33, 34, 35, 36, 37]. Non-linear
effects in transport have also been studied within this formalism [38, 39, 40,
41, 42, 43].

In the linear response regime of thermoelectric transport, the electric
current and heat current are given by the expressions,

$~{}\begin{aligned}
\bm{j}^{e}&=\stackrel{{\scriptstyle\leftrightarrow}}{{L}}_{11}\cdot\left(\bm{E}+\frac{\gradient{\mu}}{e}\right)+\stackrel{{\scriptstyle\leftrightarrow}}{{L}}_{12}\cdot(-\gradient{T})\\\
\bm{j}^{Q}&=\stackrel{{\scriptstyle\leftrightarrow}}{{L}}_{21}\cdot\left(\bm{E}+\frac{\gradient{\mu}}{e}\right)+\stackrel{{\scriptstyle\leftrightarrow}}{{L}}_{22}\cdot(-\gradient{T}),\end{aligned}$
(1)

where the heat current $\bm{j}^{Q}$ is related to the energy current
$\bm{j}^{E}$ and the number density current $\bm{j}^{N}$ by the relation
$\bm{j}^{Q}=\bm{j}^{E}-\mu\bm{j}^{N}$ [44], $\mu$ being the chemical
potential, and the electric current $\bm{j}^{e}$ is related to the number
density current $\bm{j}^{N}$ by the relation111Charge of an electron is taken
to be $-e$ in this paper. $\bm{j}^{e}=-e\bm{j}^{N}$. Here
$\stackrel{{\scriptstyle\leftrightarrow}}{{L}}_{11}$ is the electric
conductivity tensor, $\stackrel{{\scriptstyle\leftrightarrow}}{{L}}_{22}$, the
thermal conductivity tensor222See Appendix A for more discussion about this
quantity., and $\stackrel{{\scriptstyle\leftrightarrow}}{{L}}_{12}$ and
$\stackrel{{\scriptstyle\leftrightarrow}}{{L}}_{21}$ are the Peltier
conductivity coefficients. The coefficients
$\stackrel{{\scriptstyle\leftrightarrow}}{{L}}_{12}$ and
$\stackrel{{\scriptstyle\leftrightarrow}}{{L}}_{21}$ are related via the
Onsager relation [45, 46, 47],

$\stackrel{{\scriptstyle\leftrightarrow}}{{L}}_{21}=T\stackrel{{\scriptstyle\leftrightarrow}}{{L}}_{12}.$
(2)

It is important to note that the Onsager relation holds only for the
coefficients that relate transport currents to the applied gradients and not
the total currents. The total currents for systems with broken time reversal
symmetry typically have diamagnetic contributions, which need to be subtracted
in order to obtain the transport currents [47]. This has been explicitly
demonstrated for systems with a non-zero Berry curvature [48, 11].

An additional set of relations are the Einstein relations which require the
transport currents generated by an electric field $\bm{E}$ to be equal to
those due to a chemical potential gradient $\gradient{\mu}$ of strength
$e\bm{E}$ [48, 47]. It is well known that the bound electric current in a
magnetic sample can be expressed as the curl of the charge magnetization [49].
Similarly, there can be circulating energy currents in a sample, which can be
expressed as the curl of a quantity called the energy magnetization [47]. It
has been demonstrated that the Einstein relation holds for the electric
current in systems with a non-zero Berry curvature [48]. In this paper, we
show that for the Einstein relation to hold for the energy current or the heat
current in the above type of system, the energy magnetization has to satisfy a
certain specific condition (Eq. (18)). We provide an interpretation of this
condition for a Chern insulator. While the charge magnetization due to Berry
curvature can be calculated [48] from first principles, we demonstrate that
the same expression can be obtained by demanding that the Einstein relation
holds for the electric transport current. Employing a similar procedure for
the transport energy current, we show that an expression for the energy
magnetization can be found once the Bloch wavefunctions and the band structure
of the system are known.

As will be shown in section VI, the transport currents contain pieces which
depend on the equilibrium and the non-equilibrium distribution functions. The
equilibrium distribution function (Fermi function) contributes to the
intrinsic anomalous Hall and Nernst responses, which occur without any
external magnetic field, for a system with non-zero Berry curvature. The non-
equilibrium distribution function, which is calculated in section VI up to
leading order in the external potential and temperature gradients for a two-
dimensional system, is responsible for the regular Hall and Nernst responses
and captures the effects of Berry curvature on these. Similar expressions for
the non-equilibrium part of the distribution were obtained in several papers
[18, 50, 51], but only in the context of chiral magnetic effects, which are
absent in two dimensions. When the Berry curvature $\bm{\Omega}=0$, the non-
equilibrium parts of the distributions obtained in the above-mentioned papers
only lead to regular Ohmic transport, but not the regular Hall effect.
Neglecting parts of the non-equilibrium distribution function is justified for
systems like Weyl semimetals, where the anomalous Hall effect is much stronger
than the regular Hall effect. However, there are systems (e.g., bilayer
graphene) with non-zero Berry curvature, which do not show the anomalous Hall
effect due to intrinsic time reversal symmetry. In such samples, the Berry
curvature modifies the regular Hall effect, which can be calculated with the
solution (Eq. (27)) of the Boltzmann Transport Equation obtained in this
paper.

There are two main results of this paper. The first, obtained through the
calculations of sections II-V, is a derivation of the charge and energy
magnetization assuming the validity of the Einstein relation. The second
result, obtained in section VI, is a complete solution of the Boltzmann
Transport Equation in two-dimensions in the linear response regime, which can
be used to calculate transport currents for a two-dimensional system with a
non-zero Berry curvature and in the presence of a magnetic field.

The paper is organized as follows: We describe the overall formalism of
calculating currents from the semiclassical equations in section II and then
in section III, we discuss how the orbital magnetization affects the electric
and energy currents. In section IV, we demonstrate how we can obtain the known
expression of the charge magnetization by demanding that the Einstein relation
holds for the electric transport current. Section V contains one of the two
central results of this paper. Here we employ the same method as for the
charge magnetization to find a condition that the energy magnetization has to
obey. We also provide a physical interpretation of this condition and then use
it to find the expression of the energy magnetization. In section VI, we
derive the complete expressions for the transport heat current density, and
solve the Boltzmann transport equation for the electron distribution function
up to linear order in the potential and temperature gradient in the presence
of a magnetic field and non-zero Berry curvature. As mentioned above, the
distribution function can be used to obtain the regular Hall and Nernst
responses.

##  II Formalism

The semiclassical equations of motion for the position and crystal momentum of
a Bloch wavepacket are333In this paper, the charge of an electron is $-e$. [3,
4, 7],

$\dot{\bm{r}}=\frac{1}{\hbar}\frac{\partial\varepsilon_{{}_{\bm{k}}}}{\partial\bm{k}}-\dot{\bm{k}}\times\bm{\Omega}(\bm{k})$
(3) $\hbar\dot{\bm{k}}=-e\left(\bm{E}+\dot{\bm{r}}\times\bm{B}\right)$ (4)

Here,
$\bm{\Omega}(\bm{k})=i\bra{\bm{\nabla}_{{}_{\bm{k}}}u_{n,\bm{k}}}\times\ket{\bm{\nabla}_{{}_{\bm{k}}}u_{n,\bm{k}}}$
is the Berry curvature in the reciprocal space, and $u_{n,\bm{k}}$ is the
periodic part of the Bloch wavefunction. The energy eigenvalues are modified
due to the orbital magnetic moment [48] $\bm{m}_{{}_{\bm{k}}}$ of an
wavepacket444See section III for more discussion about the orbital magnetic
moment.,
$\varepsilon_{{}_{\bm{k}}}=\varepsilon_{0}(\bm{k})-\bm{m}_{{}_{\bm{k}}}\cdot\bm{B}$,
where $\varepsilon_{0}(\bm{k})$ is the original band structure energy at zero
magnetic field.

The coupled differential equations of motion for $\dot{\bm{r}}$ and
$\dot{\bm{k}}$ can be decoupled with some simple vector algebra [50], and the
resulting decoupled equations are,

$~{}\dot{\bm{r}}=\frac{\frac{1}{\hbar}\frac{\partial\varepsilon_{\bm{k}}}{\partial\bm{k}}+\frac{e}{\hbar}(\bm{E}\times\bm{\Omega})+\frac{e}{\hbar^{2}}(\bm{\Omega}\cdot\frac{\partial\varepsilon_{\bm{k}}}{\partial\bm{k}})\bm{B}}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}$
(5)

$~{}\dot{\bm{k}}=-\frac{\frac{e}{\hbar}\bm{E}+\frac{e}{\hbar^{2}}\frac{\partial\varepsilon_{\bm{k}}}{\partial\bm{k}}\times\bm{B}+\frac{e^{2}}{\hbar^{2}}(\bm{E}\cdot\bm{B})\bm{\Omega}}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}$
(6)

In a two-dimension sample, the Berry curvature $\bm{\Omega}$ is along the $z$
axis, while $\frac{\partial\varepsilon_{\bm{k}}}{\partial\bm{k}}$ is in the
$xy$ plane, then the quantity
$\bm{\Omega}\cdot\frac{\partial\varepsilon_{\bm{k}}}{\partial\bm{k}}=0$. Also,
in two-dimensions, $\bm{\Omega}$ is along the $z$ axis, whereas $\bm{k}$ is
restricted to the $xy$ plane. Then we can neglect the term
$(\bm{E}\cdot\bm{B})\bm{\Omega}$ in the expression of $\dot{\bm{k}}$. Then the
equations (5) and (6) simplify to,

$~{}\dot{\bm{r}}=\frac{\frac{1}{\hbar}\frac{\partial\varepsilon_{\bm{k}}}{\partial\bm{k}}+\frac{e}{\hbar}(\bm{E}\times\bm{\Omega})}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}$
(7)

$~{}\dot{\bm{k}}=-\frac{\frac{e}{\hbar}\bm{E}+\frac{e}{\hbar^{2}}\frac{\partial\varepsilon_{\bm{k}}}{\partial\bm{k}}\times\bm{B}}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}$
(8)

It has been shown that the equations of motion (Eq. (3) and Eq. (4)) violate
Liouville’s theorem, and for a phase space volume element $\Delta V$, the
quantity $\Delta
V\left(1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}(\bm{k})\right)$ remains a
constant of motion [52]. As a result, the calculation of the expectation value
of any operator $\hat{\mathcal{O}}$ over all states requires the introduction
of an additional factor of
$\left(1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}\left(\bm{k}\right)\right)$ in
the integrand.

$~{}\expectationvalue{\hat{\mathcal{O}}}=2\int\frac{d\bm{k}}{(2\pi)^{d}}\left(1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}\left(\bm{k}\right)\right)\expectationvalue{\hat{\mathcal{O}}}_{{}_{\bm{k}}}\tilde{f}_{{}_{\bm{k}}}$
(9)

The factor of 2 is for the spin degeneracy555In principle, the energy
eigenvalues will also be modified due to Zeeman coupling of the spin magnetic
moment with the magnetic field. But that will not change any of the main
results, and we can neglect it.. Here, $\tilde{f}_{{}_{\bm{k}}}$ is the
distribution function and in thermodynamic equilibrium, it reduces to the
Fermi distribution function,
$f_{{}_{\bm{k}}}=\frac{1}{e^{\beta\left(\varepsilon_{{}_{\bm{k}}}-\mu\right)}+1}$.
When an external electromagnetic field, chemical potential gradient or
temperature gradient is applied, the distribution function is modified to
$\tilde{f}_{{}_{\bm{k}}}=f_{{}_{\bm{k}}}+g_{{}_{\bm{k}}}$, where
$g_{{}_{\bm{k}}}$ is the non-equilibrium contribution. We calculate
$g_{{}_{\bm{k}}}$ for a two dimensional system within the Boltzmann transport
formalism in Section VI. Up to linear order in the external fields,
$g_{{}_{\bm{k}}}$ contributes to regular Ohmic conduction, the regular Hall
effect, and the regular Nernst effect. In the presence of a non-zero Berry
curvature, the equilibrium part of the distribution contributes to the
anomalous Hall effect and the anomalous Nernst effect, up to linear order in
the external electric field, temperature gradient, and chemical potential
gradient.

##  III Orbital Magnetic Moment

We first briefly review the calculation of the orbital magnetic moment of a
Bloch wavepacket, which is responsible for the circulating magnetization
energy currents and electric currents. While a Bloch wavepacket [7, 4] is
localized at a point (say, $\bm{r}_{0}$), the electron is not necessarily
localized there, and can have an angular momentum due to the motion about the
center of the wavepacket, giving rise to an orbital magnetic moment given by
the expression [48, 53, 54, 55],

$~{}\bm{m}_{\bm{k}}=-\frac{e}{2m}\bra{\psi_{\bm{k},\bm{r}_{0}}}(\hat{\bm{r}}-\bm{r}_{0})\times\hat{\bm{p}}\ket{\psi_{\bm{k},\bm{r}_{0}}}.$
(10)

Note that $m$ is the bare electron mass, not the effective mass of the Bloch
state.

Figure 1: In addition to its velocity, a Bloch wavepacket can have an angular
momentum about its center, which gives rise to an orbital magnetic moment.

In the absence of an external magnetic field, the total electric current
density is given by666Also see Appendix B and Appendix C for a self-contained
derivation of how the orbital magnetic moment appears in Eq. (11). [11, 32],

$~{}\begin{aligned} \left.\bm{j}^{e}_{\text{total
}}\right|_{\bm{B}=0}&=-e\int\frac{2d\bm{k}}{(2\pi)^{d}}\left[{g}_{\bm{k}}\frac{1}{\hbar}\frac{\partial\varepsilon_{0}({{\bm{k}}})}{\partial\bm{k}}+{f}_{\bm{k}}\frac{e}{\hbar}(\bm{E}\times\bm{\Omega}(\bm{k}))\right]\\\
&+\bm{\nabla}\times\int\frac{2d\bm{k}}{(2\pi)^{d}}f_{\bm{k}}\bm{m}_{\bm{k}}\end{aligned}$
(11)

Similarly, in the absence of a magnetic field777If there is a non-zero
magnetic field, we only need to make a handful of changes to Eq. (11) and Eq.
(12). The energy eigenvalues of the Bloch eigenstates will be modified as
$\varepsilon_{\bm{k}}=\varepsilon_{0}(\bm{k})-\bm{m}_{\bm{k}}\cdot\bm{B}$, and
the distribution functions $f_{\bm{k}}$ and $g_{\bm{k}}$ will also be modified
(see section VI). Also, the Jacobian factor
$(1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega})$ would not be canceled in the
expressions containing $\bm{m}_{\bm{k}}$, and we would have terms like
$\int\frac{2d\bm{k}}{(2\pi)^{d}}(1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega})f_{\bm{k}}\bm{m}_{\bm{k}}$
in Eq. (11) and Eq. (12), and the rest of the calculation would be identical
up to leading order in the magnetic field., the total energy current can be
shown to be [48],

$\displaystyle\left.\bm{j}^{E}_{\text{total }}\right|_{\bm{B}=0}$
$\displaystyle=\int\frac{2d\bm{k}}{(2\pi)^{d}}\varepsilon_{0}(\bm{k})\left[{g}_{\bm{k}}\frac{1}{\hbar}\frac{\partial\varepsilon_{0}(\bm{k})}{\partial\bm{k}}+\frac{e}{\hbar}{f}_{\bm{k}}(\bm{E}\times\bm{\Omega}(\bm{k}))\right]$
(12)
$\displaystyle-\bm{E}\times\int\frac{2d\bm{k}}{(2\pi)^{d}}f_{\bm{k}}\bm{m}_{\bm{k}}$

In both Eq. (11) and Eq. (12), the first term that does not involve the energy
magnetization $\bm{m_{k}}$ can be considered as the contribution due to the
movement of the center of the Bloch wavepacket, while the second term is the
contribution of the movement of the electron about the center of the
wavepacket.

##  IV Obtaining the expression for charge magnetization from the Einstein
relation of the electric current

In this section we demonstrate that the known expression of the charge
magnetization can be recovered by demanding that the Einstein relation holds
for the electric transport current. The electric transport current can be
obtained by subtracting the charge magnetization current from the total
electric current. The magnetization electric current is the curl of the charge
magnetization, whose expression has been obtained in Ref.[48]. Here we show
that the same expression of the charge magnetization can be obtained by
demanding that the Einstein relation holds for the electric transport current.
In the next section, we will employ the same strategy to find the expression
for the energy magnetization (whose curl is the bound energy current). The
expression for the electric magnetization current is,

$\bm{j}^{e}_{M}=\curl{\bm{M}^{e}}=\gradient{\mu}\times\frac{\partial\bm{M}^{e}}{\partial\mu}+\gradient{T}\times\frac{\partial\bm{M}^{e}}{\partial
T}.$

Consequently, we can calculate the transport current
$\bm{j}^{e}_{\text{transport}}=\bm{j}^{e}_{\text{total}}-\bm{j}^{e}_{M}$, and
demanding that the terms can only depend on the combination
$(e\bm{E}+\gradient{\mu})$, we get the condition (see Appendix D),

$\frac{\partial\bm{M}^{e}}{\partial\mu}=\int\frac{2d\bm{k}}{(2\pi)^{d}}\left[\frac{\partial
f_{\bm{k}}}{\partial\mu}\bm{m}_{\bm{k}}+{f}_{\bm{k}}\frac{e}{\hbar}\bm{\Omega}(\bm{k})\right].$
(13)

After integrating and applying appropriate boundary conditions888The charge
magnetization should be zero for an empty band., we get,

$\bm{M}^{e}=\int\frac{2d\bm{k}}{(2\pi)^{d}}\left[\bm{m}(\bm{k})f_{\bm{k}}+k_{B}T\frac{e\bm{\Omega}}{\hbar}\log(1+e^{-\beta(\varepsilon_{0}(\bm{k})-\mu)})\right]$
(14)

which is the same as the expression obtained in Ref.[48].

##  V Obtaining the expression for energy magnetization from the Einstein
relation of the energy current

The expression in Eq. (12) is the total energy current density, which is due
to the sum of the transport and magnetization parts of the energy current. The
magnetization energy current is the curl of energy magnetization, which is
[47], $\bm{M}^{E}=\bm{M}_{0}^{E}-(\bm{E}\cdot\bm{r})\bm{M}^{e}$. Here
$\bm{M}_{0}^{E}$ is the energy magnetization at zero external electric field,
and $\bm{M}^{e}$ is the usual charge magnetization999In Ref.[47], the authors
use the notation $\bm{M}^{E}=\bm{M}_{0}^{E}+\phi(\bm{r})\bm{M}^{N}$, where
$\bm{M}^{N}$ (whose curl is bound number density current) is
$\frac{\bm{M}^{e}}{-e}$, and $\phi(\bm{r})=e\bm{E}\cdot\bm{r}$ is the
potential energy due to the (constant) electric field. (whose curl is the
bound electric current density).

In the presence of an electric field, we need to add this additional term
because the charges carry the potential whose gradient gives rise to the
field. Both the bare energy magnetization $\bm{M}_{0}^{E}$ and the charge
magnetization $\bm{M}^{e}$ are functions of the chemical potential $\mu$ and
temperature $T$. Then, the circulating, bound energy current density is,

$~{}\begin{aligned}
\bm{j}^{E}_{M}=\curl{\bm{M}_{0}^{E}}-\bm{E}\times\bm{M}^{e}-\underbrace{(\bm{E}\cdot\bm{r})\curl{\bm{M}^{e}}}_{\text{2nd
order quantity}}\end{aligned}$ (15)

Here $\curl{\bm{M}^{e}}$ is already a function of $\gradient{\mu}$ and
$\gradient{T}$, and hence, the term $(\bm{E}\cdot\bm{r})\curl{\bm{M}^{e}}$ is
of second order, and we drop it. It follows that the transport energy current
density is (see Appendix E),

$\displaystyle\bm{j}_{\text{transport }}^{E}$
$\displaystyle=\bm{j}_{\text{total }}^{E}-\bm{j}^{E}_{M}$ (16)
$\displaystyle=\int\frac{2d\bm{k}}{(2\pi)^{d}}\varepsilon_{0}(\bm{k})\frac{1}{\hbar}\frac{\partial\varepsilon_{0}(\bm{k})}{\partial\bm{k}}{g}_{\bm{k}}$
$\displaystyle+e\bm{E}\times\bm{V}_{1}-\gradient{\mu}\times\frac{\partial\bm{M}_{0}^{E}}{\partial\mu}-\gradient{T}\times\frac{\partial\bm{M}_{0}^{E}}{\partial
T}$

where

$\bm{V}_{1}=\int\frac{2d\bm{k}}{(2\pi)^{d}}\frac{\bm{\Omega}(\bm{k})}{\hbar}\left[\varepsilon_{0}f_{\bm{k}}+k_{B}T\log(1+e^{-\beta(\varepsilon_{0}(\bm{k})-\mu)})\right].$
(17)

In Eq. (16), there is an asymmetry between $e\bm{E}$ and $\gradient{\mu}$. For
the Einstein relation to hold for the energy current, the allowed terms can
only depend on the combination $(e\bm{E}+\gradient{\mu})$. As shown in Sec.VI,
$g_{\bm{k}}$ is proportional to this particular combination. Thus, for the
Einstein relation to hold,

$\displaystyle-\frac{\partial\bm{M}_{0}^{E}}{\partial\mu}=\int$
$\displaystyle\frac{2d\bm{k}}{(2\pi)^{d}}\frac{\bm{\Omega}(\bm{k})}{\hbar}\bigg{[}\varepsilon_{0}(\bm{k})f_{\bm{k}}$
(18)
$\displaystyle+k_{B}T\log(1+e^{-\beta(\varepsilon_{0}(\bm{k})-\mu)})\bigg{]}.$

When the temperature scale ($k_{B}T$) is negligible compared to the chemical
potential $\mu$ (more precisely, negligible compared to
$(\mu-\varepsilon_{\text{min}})$, where $\varepsilon_{\text{min}}$ is the band
minima), which happens for ordinary metals in room temperature, it can be
shown that (See Appendix F) up to leading order (zeroth order) in temperature,

$-\frac{\partial\bm{M}_{0}^{E}}{\partial\mu}=\frac{\mu}{\hbar}\int\frac{2d\bm{k}}{(2\pi)^{d}}\bm{\Omega}(\bm{k})f_{\bm{k}}$

When there is more than one band, the result would be,

$-\frac{\partial\bm{M}_{0}^{E}}{\partial\mu}=\sum_{n}\frac{\mu}{\hbar}\int\frac{2d\bm{k}}{(2\pi)^{d}}\bm{\Omega}_{n}(\bm{k})f_{n\bm{k}},$
(19)

where $n$ is the band index. This has a particularly simple interpretation in
a two-dimensional sample. The magnitude of the circulating, bound energy
current along an edge (not to be confused with the energy current density
$\bm{j}^{E}$, which resides in the bulk) is given by $I^{E}_{0}=\Delta
z|\bm{M}^{E}_{0}\times\hat{n}|$ (see Fig. 2).

Figure 2: The circulating current and the magnetization in an effectively two-
dimensional sample of thickness $\Delta z$ are related as $I=\Delta
z|\bm{M}\times\hat{n}|$.

Then, the result can be rewritten as,

$\frac{\Delta I_{0}^{E}}{\Delta\mu}=\Delta
z\sum_{n}\frac{\mu}{\hbar}\int\frac{2d\bm{k}}{(2\pi)^{2}}\bm{\Omega}_{n}(\bm{k})f_{n\bm{k}}.$

A similar result for the charge magnetization and the circulating electric
current was shown in Refs.[54, 56], using the known expression [48, 53, 54]
for the charge magnetization. Note that we obtained this condition on the
energy magnetization current from a necessary condition (Eq. (18)) so that the
Einstein relation would hold. Let us discuss several cases, similar to the
discussion on the circulating electric current by R. Resta [56].

For a Chern insulator, the integral of the Berry curvature over the filled
bands (with all the factors of $2\pi$ included) is the sum of the Chern
numbers of those bands, which can be interpreted as the number of chiral edge
states [57, 58]. We get,

$\frac{\Delta I_{0}^{E}}{\Delta\mu}=\Delta
z\left(\frac{\mu}{h}\right)\left|N_{\circlearrowleft}-N_{\circlearrowright}\right|,$
(20)

where $N_{\circlearrowleft}$ ($N_{\circlearrowright}$) is the number of chiral
edge states circulating along the anti-clockwise (clockwise) direction. As the
chemical potential is varied, per unit change in chemical potential, each of
the chiral edge states contributes to the the circulating energy current by an
amount $\frac{\mu}{h}$ per unit thickness of the sample101010In Eq. (19) and
Eq. (20), the absolute value of the chemical potential appears. Its
implications are discussed in Appendix G..

When there are no chiral edge states (e.g. in a trivial insulator), this
quantity is zero. Also, when there are equal number of chiral edge states
circulating in the clockwise and counter-clockwise direction, their
contributions cancel.

In a metallic sample, the result would be similar to Eq. (20) for the filled
bands, but if the conduction band (which is partially filled) is of
topological nature, it would also partially contribute.

It has been shown [55, 56] that the derivative of the charge magnetization
with respect to the chemical potential has a term,
$\frac{e}{\hbar}\sum_{n}\int\frac{d\bm{k}}{(2\pi)^{d}}f_{n\bm{k}}\bm{\Omega}_{n}(\bm{k})$.
The sign difference with Eq. (19) is because, when an electron moves, the
local energy current density and the electric current density would be in
opposite directions, as the electron bears a negative charge.

While we do not have a closed expression for the energy magnetization (unlike
that of charge magnetization, whose closed expressions may be found in
Refs.[55, 48]), it can be calculated by numerically integrating Eq. (19) (or
Eq. (18), for finite temperature corrections),

$\bm{M}_{0}^{E}=-\int_{\tilde{\mu}=\varepsilon_{min}}^{\mu}d{\tilde{\mu}}\sum_{n}\frac{\tilde{\mu}}{\hbar}\int\frac{2d\bm{k}}{(2\pi)^{d}}\bm{\Omega}_{n}(\bm{k})f_{n\bm{k}}(\tilde{\mu}),$
(21)

with the initial condition that the energy magnetization would be zero when
the chemical potential is at the band minima $\varepsilon_{min}$, as then
there would be no electrons to generate the circulating currents.

##  VI Complete expressions for transport currents in a two-dimensional sample
with non-zero Berry curvature

The electric transport current density can be obtained by subtracting the curl
of the charge magnetization from the total electric current density [48],

$~{}\begin{aligned}
{{\bm{j}}^{e}}_{\text{transport}}&=-e\int\frac{2d\bm{k}}{(2\pi)^{d}}{g}_{\bm{k}}\frac{1}{\hbar}\frac{\partial\varepsilon_{0}({{\bm{k}}})}{\partial\bm{k}}\\\
&-\int\frac{2d\bm{k}}{(2\pi)^{d}}{f}_{\bm{k}}\frac{e^{2}}{\hbar}\left(\left[\bm{E}+\frac{\gradient{\mu}}{e}\right]\times\bm{\Omega}(\bm{k})\right)\\\
&-\frac{\gradient{T}}{T}\times\int\frac{2d\bm{k}}{(2\pi)^{d}}\frac{e}{\hbar}\bm{\Omega}(\bm{k})\bigg{[}f_{\bm{k}}(\varepsilon_{0}(\bm{k})-\mu)\\\
&+k_{B}T\log(1+e^{-\beta(\varepsilon_{0}(\bm{k})-\mu)})\bigg{]}\end{aligned}$
(22)

Also, the heat current density111111There is no “heat magnetization” analogous
to the charge magnetization or the energy magnetization [59]. is (See Appendix
E),

$\displaystyle\bm{j}_{\text{transport }}^{Q}$
$\displaystyle=\bm{j}_{\text{transport }}^{E}-\mu\bm{j}_{\text{transport
}}^{N}$ (23)
$\displaystyle=\int\frac{2d\bm{k}}{(2\pi)^{d}}(\varepsilon_{0}(\bm{k})-\mu)\frac{1}{\hbar}\frac{\partial\varepsilon_{0}(\bm{k})}{\partial\bm{k}}g_{k}$
$\displaystyle+e\bm{E}\times\bm{V}_{2}+\gradient{\mu}\times\bm{V}_{3}-\frac{\gradient{T}}{T}\times\bm{V}_{4}$
$\displaystyle-\gradient{\mu}\times\frac{\partial\bm{M}_{0}^{E}}{\partial\mu}-\gradient{T}\times\frac{\partial\bm{M}_{0}^{E}}{\partial
T},$

where

$\displaystyle\bm{V}_{2}=\int$
$\displaystyle\frac{2d\bm{k}}{(2\pi)^{d}}\frac{\bm{\Omega}(\bm{k})}{\hbar}\bigg{[}(\varepsilon_{0}(\bm{k})-\mu)f_{\bm{k}}$
(24)
$\displaystyle+k_{B}T\log(1+e^{-\beta(\varepsilon_{0}(\bm{k})-\mu)})\bigg{]},$

$\bm{V}_{3}=-\int\frac{2d\bm{k}}{(2\pi)^{d}}\frac{\mu}{\hbar}\bm{\Omega}(\bm{k})f_{\bm{k}}$,
and $\bm{V}_{4}=\mu\bm{V}_{2}$. Note that if we demand Einstein relation to
hold for the heat current, we recover Eq. (18). Also note that
$\bm{V}_{2}\rightarrow 0$ as $T\rightarrow 0$, implying
$\stackrel{{\scriptstyle\leftrightarrow}}{{L}}_{21}=0$ at $T=0$, which is
consistent with the Onsager relation, Eq. (2).

Once the band structure is known, the value of $\bm{M}^{E}_{0}$ can be found
by numerically integrating Eq. (18), and $\bm{\Omega}(\bm{k})$ can be obtained
from the Bloch wavefunctions. Then, if we can find the expression of the out
of equilibrium part of the distribution function, $g_{\bm{k}}$, the transport
parts of the electric and the heat currents can be calculated up to linear
order. We can calculate $g_{\bm{k}}$ with the Boltzmann transport equation
[60]. The total distribution function
$\tilde{f}_{{}_{\bm{k}}}=(f_{\bm{k}}+g_{\bm{k}})$ satisfies the equation,

$\left[\frac{\partial}{\partial
t}+\dot{\bm{k}}\cdot\frac{\partial}{\partial\bm{k}}+\dot{\bm{r}}\cdot\frac{\partial}{\partial\bm{r}}\right](f_{{}_{\bm{k}}}+g_{{}_{\bm{k}}})=-\frac{g_{{}_{\bm{k}}}}{\tau_{{}_{\bm{k}}}}.$

In general, the scattering timescale $\tau_{{}_{\bm{k}}}$ may depend on the
crystal momentum of the electron. However, only the electrons near the Fermi
surface can be effectively scattered, and if the scattering timescale only
depends on the energy of the Bloch wavefunction, then we can consider it to be
a constant. The following calculation remains valid even when
$\tau_{{}_{\bm{k}}}$ varies with the crystal momentum. When the external
fields are time independent, we can drop the explicit partial derivatives with
respect to time, and the equation takes the form,

$\frac{g_{{}_{\bm{k}}}}{\tau_{{}_{\bm{k}}}}+\dot{\bm{k}}\cdot\frac{\partial}{\partial\bm{k}}g_{{}_{\bm{k}}}+\dot{\bm{r}}\cdot\frac{\partial}{\partial\bm{r}}g_{{}_{\bm{k}}}=-\dot{\bm{k}}\cdot\frac{\partial}{\partial\bm{k}}f_{{}_{\bm{k}}}-\dot{\bm{r}}\cdot\frac{\partial}{\partial\bm{r}}f_{{}_{\bm{k}}}$
(25)

Using the decoupled equations of motion in two-dimensions (Eq. (7) and Eq.
(8)), the right hand side of the above equation simplifies to (see Appendix
H.1),

$\frac{\partial
f_{{}_{\bm{k}}}}{\partial\varepsilon_{{}_{\bm{k}}}}\frac{1}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}(\bm{k})}\frac{1}{\hbar}\frac{\partial\varepsilon_{{}_{\bm{k}}}}{\partial\bm{k}}\cdot\left[e\bm{E}+\bm{\nabla}{\mu}+\bm{\nabla}T\frac{\varepsilon_{{}_{\bm{k}}}-\mu}{T}\right].$

Here we can exchange $e\bm{E}\leftrightarrow\bm{\nabla}\mu$, i.e., the
Einstein relation holds for the currents generated by the non-equilibrium
distribution. We will show that the Onsager relation holds as well. We assume
that the fields are constant in space, and drop the term containing
$\frac{\partial{g_{{}_{\bm{k}}}}}{\partial\bm{r}}$ 121212When we solve the
equation after discarding the term
$\dot{\bm{r}}\cdot\frac{\partial{g_{{}_{\bm{k}}}}}{\partial\bm{r}}$, and
substitute the solution (see Eq. (27)) in this term, we would find that this
term is second order in the applied fields. Then up to linear order, we can
discard it.. Then, the BTE takes the form,

$\displaystyle\frac{g_{{}_{\bm{k}}}}{\tau_{{}_{\bm{k}}}}-\frac{\frac{e}{\hbar}\bm{E}+\frac{e}{\hbar^{2}}\frac{\partial\varepsilon_{{}_{\bm{k}}}}{\partial\bm{k}}\times\bm{B}}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}$
$\displaystyle\cdot\frac{\partial}{\partial\bm{k}}g_{{}_{\bm{k}}}=\frac{\partial
f}{\partial\varepsilon_{\bm{k}}}\frac{1}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}$
$\displaystyle\frac{1}{\hbar}\frac{\partial\varepsilon_{{}_{\bm{k}}}}{\partial\bm{k}}\cdot\left[e\bm{E}+\bm{\nabla}{\mu}+\bm{\nabla}T\frac{\varepsilon_{{}_{\bm{k}}}-\mu}{T}\right].$

We can further simplify this equation. If we neglect the term
$\bm{E}\cdot\frac{\partial g_{{}_{\bm{k}}}}{\partial\bm{k}}$ in the LHS, we
would find (in the next section) that $g_{{}_{\bm{k}}}$ is a linear function
of the electric field, the chemical potential gradient, and the temperature
gradient. Then, $\bm{E}\cdot\frac{\partial g_{{}_{\bm{k}}}}{\partial\bm{k}}$
would be quadratic in the applied fields, and the solution would remain self-
consistent up to the first order if we neglect it. Finally, the equation takes
the form,

$~{}\begin{aligned}
\frac{g_{{}_{\bm{k}}}}{\tau_{{}_{\bm{k}}}}-\underbrace{\frac{\frac{e}{\hbar^{2}}\frac{\partial\varepsilon_{{}_{\bm{k}}}}{\partial\bm{k}}\times\bm{B}}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\cdot\frac{\partial}{\partial\bm{k}}g_{{}_{\bm{k}}}}_{\text{treated
as a perturbation}}&=\frac{\partial
f}{\partial\varepsilon_{\bm{k}}}\frac{1}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\\\
&\frac{1}{\hbar}\frac{\partial\varepsilon_{{}_{\bm{k}}}}{\partial\bm{k}}\cdot\left[e\bm{E}+\gradient{\mu}+\gradient{T}\frac{\varepsilon_{{}_{\bm{k}}}-\mu}{T}\right]\end{aligned}$
(26)

We treat the second term in the LHS as a perturbation. An order of the
magnitude estimate in Appendix I shows that the magnetic field needs to be of
the order of hundreds of Tesla for this term to be of the same magnitude as
the other term, $\frac{g_{{}_{\bm{k}}}}{\tau_{{}_{\bm{k}}}}$. Since the
magnetic field in the laboratories are much less than that, we can treat it as
a perturbation. We would see that this term gives rise to the (ordinary) Hall
effect in samples with a non-zero Berry curvature, and that is why we would
not completely neglect it. Up to linear order in external fields, the solution
of Eq. (26) is (see Appendix H),

$\displaystyle g_{{}_{\bm{k}}}$ $\displaystyle=\frac{\partial
f}{\partial\varepsilon_{{}_{\bm{k}}}}\frac{1}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\frac{\tau_{{}_{\bm{k}}}}{\hbar}\frac{\partial\varepsilon_{{}_{\bm{k}}}}{\partial\bm{k}}\cdot\bm{S}$
(27)
$\displaystyle+\frac{\frac{e\tau_{{}_{\bm{k}}}}{\hbar^{2}}}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\frac{\partial\varepsilon_{{}_{\bm{k}}}}{\partial\bm{k}}\cdot\bm{B}\times\left[\frac{\partial}{\partial\bm{k}}\left[\frac{\frac{\partial
f}{\partial\varepsilon_{{}_{\bm{k}}}}\tau_{{}_{\bm{k}}}}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\right]\left(\frac{\bm{S}}{\hbar}\cdot\frac{\partial\varepsilon_{{}_{\bm{k}}}}{\partial\bm{k}}\right)\right]$
$\displaystyle+\frac{\frac{e\tau_{{}_{\bm{k}}}}{\hbar^{2}}}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\frac{\partial\varepsilon_{{}_{\bm{k}}}}{\partial\bm{k}}\cdot\bm{B}\times\left[\frac{\frac{\partial
f}{\partial\varepsilon_{{}_{\bm{k}}}}\tau_{{}_{\bm{k}}}}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\left[\left(\frac{1}{\hbar}\bm{S}\cdot\frac{\partial}{\partial\bm{k}}\right)\frac{\partial\varepsilon_{{}_{\bm{k}}}}{\partial\bm{k}}\right]\right]$

where
$\bm{S}=e\bm{E}+\gradient{\mu}+\frac{\varepsilon_{{}_{\bm{k}}}-\mu}{T}\gradient{T}$.131313It
can be easily verified that when there is no Berry curvature, and the band is
parabolic with an effective mass $m^{*}$, and the scattering time $\tau$ is
independent of $\bm{k}$, this solution produces the same electric conductivity
tensor which is also obtained from the Drude model,
$\stackrel{{\scriptstyle\leftrightarrow}}{{L}}_{11}=\stackrel{{\scriptstyle\leftrightarrow}}{{\sigma}}=\frac{ne^{2}\tau}{m^{*}}\begin{pmatrix}1&-\omega_{c}\tau\\\
\omega_{c}\tau&1\end{pmatrix}$, where $\omega_{c}=\frac{eB}{m^{*}}$ is the
cyclotron frequency. Having obtained Eq. (27), $g_{\bm{k}}$ has to be
substituted in the expressions of transport currents, Eq. (22) and Eq. (23),
and the transport coefficients
$\stackrel{{\scriptstyle\leftrightarrow}}{{L}}_{11}$,
$\stackrel{{\scriptstyle\leftrightarrow}}{{L}}_{12}$,
$\stackrel{{\scriptstyle\leftrightarrow}}{{L}}_{21}$ and
$\stackrel{{\scriptstyle\leftrightarrow}}{{L}}_{22}$ can be read off. However,
the actual computation will require the band structure and the Bloch
eigenfunctions (to calculate the Berry curvature).

In this solution, $\gradient{T}$ is always multiplied with a factor of
$\frac{\varepsilon_{{}_{\bm{k}}}-\mu}{T}$, and the expression of heat current
has a factor of $(\varepsilon_{{}_{\bm{k}}}-\mu)$ (Eq. (23)), from which it
follows that the Onsager relation holds for the contribution due to the non-
equilibrium part of the distribution function in presence of Berry curvature.
It has previously been demonstrated in Ref.[48] that the Onsager relation
holds for the contribution due to the equilibrium part of the distribution.

##  VII Conclusions

The Einstein relation has been shown to hold from certain microscopic theories
for the electric transport current in systems with non-zero Berry curvature
[48]. In this paper, we employ a complementary approach to first demonstrate
that assuming that the Einstein relation holds, an expression for the charge
magnetization can be obtained in a relatively straightforward manner which
agrees with the expression obtained for this quantity previously. We then
extend this argument to the transport energy current and the heat current, to
obtain a condition that the energy magnetization has to obey. We have
interpreted this condition physically and used it to obtain an expression for
the energy magnetization, analogous to the one that exists for the charge
magnetization. We have also solved the Boltzmann transport equation up to
linear order in potential and temperature gradients for a two-dimensional
system, which can be used to obtain the regular Hall response in systems like
bilayer graphene, which possess a non-zero Berry curvature, but display no
anomalous Hall Effect due to time reversal invariance.

## Acknowledgments

A. P. acknowledges support from the KVPY programme and S. M. thanks the
Department of Science and Technology, Government of India for support. We
thank Vijay Shenoy for discussions regarding the absolute value of the
chemical potential appearing in several expressions.

##  Appendix A About the heat conductivity tensor

Here we denoted thermal conductivity as the heat current per unit temperature
gradient, at zero electric field or zero chemical potential gradient.
Sometimes, thermal conductivity is denoted as heat current per unit
temperature gradient, when the net electric current is zero (an internal
electrochemical potential gradient is generated to maintain zero electric
current, but that internal electrochemical potential gradient in turn
contributes to the heat current). In that case, it can be shown that (see Eq.
(13.56) of Ref.[44]) the heat current (which is the heat conductivity, because
we set $\gradient{T}=1$) is
$\stackrel{{\scriptstyle\leftrightarrow}}{{L}}_{22}-\stackrel{{\scriptstyle\leftrightarrow}}{{L}}_{21}(\stackrel{{\scriptstyle\leftrightarrow}}{{L}}_{11})^{-1}\stackrel{{\scriptstyle\leftrightarrow}}{{L}}_{12}$.
But this additional term is a small correction, it is of the order
$\left(\frac{k_{B}T}{\varepsilon_{F}}\right)^{2}$.

##  Appendix B Self contained derivation of Eq. (11)

The wavefunction of a Bloch wavepacket located at $r_{0}$, and peaked at
crystal momentum $k_{0}$ is [7, 4] (here we follow the notation from [7]),

$\psi(\bm{r})_{n,\bm{k}_{0}}=\int
d\bm{k}w(\bm{k}-\bm{k}_{0})e^{i\bm{A}(\bm{k_{0}})\cdot(\bm{k}-\bm{k}_{0})}e^{-i\bm{k}\cdot\bm{r}_{0}}\left(e^{i\bm{k}\cdot\bm{r}}u_{n,\bm{k}}(\bm{r})\right).$
(28)

Here $w(\bm{k}-\bm{k}_{0})$ is an (arbitrary) weight function sharply peaked
at $\bm{k}=\bm{k}_{0}$, $u_{n,\bm{k}}(\bm{r})$ is the periodic part of the
Bloch wavefunction with crystal momentum $\bm{k}$ and band index $n$, and
$\bm{A}(\bm{k})=i\bra{u_{n,\bm{k}}}\bm{\nabla}_{{}_{\bm{k}}}\ket{u_{n,\bm{k}}}$
is the Berry connection. $u_{n,\bm{k}}(\bm{r})$ is normalized such that
$\bra{u_{n,\bm{k}}}\ket{u_{n,\bm{k}}}=\int_{\text{unit
cell}}d\bm{r}|u_{n,\bm{k}}(\bm{r})|^{2}=1$. It is to be noted that the
calculations in Appendix C does not depend on the actual form of the weight
function $w(\bm{k}-\bm{k}_{0})$.

An electron in a Bloch wavepacket centered at $\bm{r}_{0}$ may be found at
another point $\bm{r}_{1}$ ($\neq\bm{r}_{0}$), and it can contribute to the
local current density at $\bm{r}_{1}$. As an analogy, the total electron
density near an atom in an insulator is not just the electron density of that
atom, but the sum of electron densities belonging to all the atoms, evaluated
at that point. Of course, atoms far away from the point would contribute very
less, but nearby atoms can contribute significantly. Let us try write down the
expectation value of a general operator $\mathcal{O}$ at a point $\bm{r}_{1}$.
We would consider the sum of the contributions from all the wavepackets
centered at a point $\bm{r}_{0}$, and sum over all $\bm{r}_{0}$ [11]. That is,

$\expectationvalue{\hat{\mathcal{O}}}(\bm{r}_{1})=\int
d\bm{r}_{0}\int\frac{2d\bm{k}}{(2\pi)^{d}}(f+g)\left(1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}\right)\bra{\psi_{\bm{k},\bm{r}_{0}}}\frac{1}{2}\\{\hat{\mathcal{O}},\delta(\bm{\hat{r}}-\bm{r}_{1})\\}\ket{\psi_{\bm{k},\bm{r}_{0}}}$

Here $\bm{\hat{r}}$ is the position operator that acts on the states.
$\bm{r}_{1}$ is just a parameter, and here $\bm{r}_{1}$ acting on a state has
to be understood as the operator $\bm{r}_{1}\hat{\mathcal{I}}$, with
$\hat{\mathcal{I}}$ being the identity operator. And,
$\frac{1}{2}\\{\hat{\mathcal{O}},\delta(\bm{\hat{r}}-\bm{r}_{0})\\}$ denotes
the Hermitized operator
$\frac{\hat{\mathcal{O}}\delta(\bm{\hat{r}}-\bm{r}_{1})+\delta(\bm{\hat{r}}-\bm{r}_{1})\hat{\mathcal{O}}}{2}$,
in case $\hat{\bm{r}}$ and $\hat{\mathcal{O}}$ do not commute (for example, to
calculate the electric current, we need to find the expectation value of the
velocity operator, which does not commute with the position operator).

Now, since the state is centered around $\bm{r}_{0}$, we can take that into
account by expanding the delta function (up to leading order)[61],

$\displaystyle\delta(\bm{r}_{1}-\bm{\hat{r}})=\delta((\bm{r}_{1}-\bm{r}_{0})-(\bm{\hat{r}}-\bm{r}_{0}))$
$\displaystyle\approx\delta(\bm{r}_{1}-\bm{r}_{0})-(\bm{\hat{r}}-\bm{r}_{0})\cdot\bm{\nabla}_{\bm{r}_{1}}\delta(\bm{r}_{1}-\bm{r}_{0})$
$\displaystyle=\delta(\bm{r}_{1}-\bm{r}_{0})-\bm{\nabla}_{\bm{r}_{1}}\cdot[(\bm{\hat{r}}-\bm{r}_{0})\delta(\bm{r}_{1}-\bm{r}_{0})]$

Using this, we get,

$\displaystyle\expectationvalue{\hat{\mathcal{O}}}(\bm{r}_{1})$
$\displaystyle=\int\frac{2d\bm{k}}{(2\pi)^{d}}(f+g)\left(1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}\right)\bra{\psi_{\bm{k},\bm{r}_{1}}}\hat{\mathcal{O}}\ket{\psi_{\bm{k},\bm{r}_{1}}}$
(29)
$\displaystyle-\bm{\nabla}_{\bm{r}_{1}}\cdot\int\frac{2d\bm{k}}{(2\pi)^{d}}(f+g)\left(1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}\right)\bra{\psi_{\bm{k},\bm{r}_{1}}}\frac{1}{2}\\{\hat{\mathcal{O}},(\bm{\hat{r}}-\bm{r}_{1})\\}\ket{\psi_{\bm{k},\bm{r}_{1}}}$

We would work with this leading order expansion of the Delta function. For the
current density operator, we need to separately take the three components of
$\bm{\hat{j}}=-e\frac{\hat{\bm{p}}}{m}$ to be $\hat{\mathcal{O}}$.

Let us define the tensor,

$\mathcal{M}_{ij}=\bra{\psi_{\bm{k},\bm{r}_{1}}}\frac{1}{2}\\{-e\frac{\hat{p}_{j}}{m},(\bm{\hat{r}}-\bm{r}_{1})_{i}\\}\ket{\psi_{\bm{k},\bm{r}_{1}}}$

where $i,j$ runs from $1,2,\dots d$. It can be shown that (see Appendix C.1)
$\mathcal{M}_{ij}$ is completely anti-symmetric [32]. Then, in 3D universe
(the sample can still be 2 dimensional), we can construct a vector, whose
components are, $\tilde{m}_{i}=\frac{1}{2}\epsilon_{ijk}\mathcal{M}_{jk}$
[32]. It can be shown (see Appendix C.2) that this is identical to the $i$-th
component of the orbital magnetic moment $\bm{m}_{\bm{k}}$, defined in Eq.
(10). We can also invert this relation,
$\mathcal{M}_{jk}=\epsilon_{ijk}m_{i}$.

Note that we can rewrite the second term as,

$\displaystyle-\bm{\nabla}_{\bm{r}_{1}}\cdot\int\frac{2d\bm{k}}{(2\pi)^{d}}(f+g)\left(1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}\right)\bra{\psi_{\bm{k},\bm{r}_{1}}}\frac{1}{2}\\{\hat{\mathcal{O}},(\bm{\hat{r}}-\bm{r}_{1})\\}\ket{\psi_{\bm{k},\bm{r}_{1}}}$
$\displaystyle=-\partial_{i}\int\frac{2d\bm{k}}{(2\pi)^{d}}(f+g)\left(1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}\right)\bra{\psi_{\bm{k},\bm{r}_{1}}}\frac{1}{2}\\{\hat{\mathcal{O}},(\bm{\hat{r}}-\bm{r}_{1})_{i}\\}\ket{\psi_{\bm{k},\bm{r}_{1}}},$

with sum over $i$ implied. Then to get the $\mu$-th component of
$\expectationvalue{\bm{\hat{j}}^{e}}$, we take
$\mathcal{\hat{O}}=-e\frac{\hat{{p}_{\mu}}}{m}$. We denote
$[dk]=\frac{2d\bm{k}}{(2\pi)^{d}}(1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega})$.
We get,

$\displaystyle-\partial_{i}\int[dk](f+g)\mathcal{M}_{i\mu}$
$\displaystyle=-\partial_{i}\int[dk](f+g)m_{\nu}\epsilon_{\nu i\mu}$
$\displaystyle=+\partial_{i}\int[dk](f+g)m_{\nu}\epsilon_{i\nu\mu}$

This is the $\mu$-th component of
$\bm{\nabla}\times\int\frac{2d\bm{k}}{(2\pi)^{d}}(f+g)(1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega})\bm{m}_{\bm{k}}$.
When the electric field, the magnetic field, the chemical potential gradient
and the temperature gradients are constant, we can drop the term
$\bm{\nabla}\times\int
2d\bm{k}g(1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega})\bm{m}_{\bm{k}}$, because
$g$ is already a linear function in the constant fields $\bm{E}$,
$\gradient{\mu}$ and $\gradient{T}$.

The first term,
$\int\frac{2d\bm{k}}{(2\pi)^{d}}(f+g)(1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega})\bra{\psi_{\bm{k},\bm{r}_{1}}}\hat{\bm{j}^{e}}\ket{\psi_{\bm{k},\bm{r}_{1}}}$
is the contribution of the center of the wavepacket, which is (up to linear
order),
$-e\int\frac{2d\bm{k}}{(2\pi)^{d}}\left[{g}_{\bm{k}}\frac{1}{\hbar}\frac{\partial\varepsilon_{{}_{\bm{k}}}}{\partial\bm{k}}+{f}_{\bm{k}}\frac{e}{\hbar}(\bm{E}\times\bm{\Omega}(\bm{k}))\right]$
(note that the phase space correction factor
$(1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega})$ gets canceled). Finally we get, up
to linear order,

$~{}\bm{j}^{e}_{\text{total
}}=\expectationvalue{\hat{\bm{j}}^{e}}=-e\int\frac{2d\bm{k}}{(2\pi)^{d}}\left[{g}_{\bm{k}}\frac{1}{\hbar}\frac{\partial\varepsilon_{{}_{\bm{k}}}}{\partial\bm{k}}+{f}_{\bm{k}}\frac{e}{\hbar}(\bm{E}\times\bm{\Omega}(\bm{k}))\right]+\bm{\nabla}\times\int\frac{2d\bm{k}}{(2\pi)^{d}}f_{\bm{k}}\left(1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}\right)\bm{m}_{\bm{k}}.$
(30)

This equation has a simple interpretation. The net current is the contribution
of the movement of the center the wavepackets, as well as due to the rotation
about their individual centers. We want to calculate the properties of the
anomalous Hall effect, where $\bm{B}=0$. Then the expression simplifies to

$~{}\left.\bm{j}^{e}_{\text{total
}}\right|_{\bm{B}=0}=-e\int\frac{2d\bm{k}}{(2\pi)^{d}}\left[{g}_{\bm{k}}\frac{1}{\hbar}\frac{\partial\varepsilon_{0}({{\bm{k}}})}{\partial\bm{k}}+{f}_{\bm{k}}\frac{e}{\hbar}(\bm{E}\times\bm{\Omega}(\bm{k}))\right]+\bm{\nabla}\times\int\frac{2d\bm{k}}{(2\pi)^{d}}f_{\bm{k}}\bm{m}_{\bm{k}}.$
(31)

##  Appendix C Certain calculations related to orbital magnetic moment

###  C.1
$\mathcal{M}_{ij}=\bra{\psi_{\bm{k},\bm{r}_{1}}}\frac{1}{2}\\{-e\frac{\hat{p}_{j}}{m},(\bm{\hat{r}}-\bm{r}_{1})_{i}\\}\ket{\psi_{\bm{k},\bm{r}_{1}}}$
is a totally antisymmetric tensor

The effective Hamiltonian[44] acting on $u_{n,\bm{k}}$ in its effective
Schrodinger equation is,
$H(\bm{k})=\left[\frac{\hbar^{2}}{2m}(\bm{k}-i\bm{\nabla})^{2}+V(\bm{r})\right]=e^{-i\bm{k}\cdot\bm{r}}\hat{H}e^{i\bm{k}\cdot\bm{r}}$,
with $\hat{H}$ being the original Hamiltonian.

Let us define $\tilde{\bm{P}}(\bm{k})=\frac{m}{\hbar}\frac{\partial
H(\bm{k})}{\partial\bm{k}}=\hbar(\bm{k}-i\bm{\nabla})=e^{-i\bm{k}\cdot\bm{r}}\hat{\bm{P}}e^{i\bm{k}\cdot\bm{r}}$.
We would first prove some preliminary results which we would need later on.

####  C.1.1 Calculation of
$\langle{u_{\bm{k}}}|\tilde{\bm{P}}(\bm{k})|{u_{\bm{k}}}\rangle$

The idea of this calculation is based on [4].

$\displaystyle\expectationvalue{\tilde{\bm{P}}(\bm{k})}=\bra{u_{\bm{k}}}\tilde{\bm{P}}(\bm{k})\ket{u_{\bm{k}}}$
$\displaystyle=\frac{m}{\hbar}\bra{u_{\bm{k}}}\frac{\partial
H(\bm{k})}{\partial\bm{k}}\ket{u_{\bm{k}}}$
$\displaystyle=\frac{m}{\hbar}\left[\bra{u}\frac{\partial}{\partial\bm{k}}\left(H(\bm{k})\ket{u}\right)-\bra{u}H(\bm{k})\ket{\frac{\partial}{\partial\bm{k}}u}\right]$
$\displaystyle=\frac{m}{\hbar}\left[\bra{u}\frac{\partial}{\partial\bm{k}}\left(\varepsilon_{{}_{\bm{k}}}\ket{u}\right)-\varepsilon_{{}_{\bm{k}}}\bra{u}\ket{\frac{\partial}{\partial\bm{k}}u}\right]$
$\displaystyle=\frac{m}{\hbar}\left[\bra{u}\left(\varepsilon_{{}_{\bm{k}}}\ket{\frac{\partial}{\partial\bm{k}}u}\right)+\bra{u}\left(\frac{\partial\varepsilon_{{}_{\bm{k}}}}{\partial\bm{k}}\ket{u}\right)-\varepsilon_{{}_{\bm{k}}}\bra{u}\ket{\frac{\partial}{\partial\bm{k}}u}\right]$
$\displaystyle=\frac{m}{\hbar}\frac{\partial\varepsilon_{{}_{\bm{k}}}}{\partial\bm{k}}$

which is just the (bare) electron mass times its velocity. Note that the Berry
curvature term does not appear here because we took a single Bloch eigenstate,
not an wavepacket.

####  C.1.2 Calculation of
$\bra{\frac{\partial}{\partial\bm{k}}u_{\bm{k}}}\tilde{\bm{P}}(\bm{k})\ket{u_{\bm{k}}}$

This calculation is based on [4].

$\displaystyle\bra{\frac{\partial}{\partial
k_{\mu}}u_{\bm{k}}}\tilde{P}_{\nu}(\bm{k})\ket{u_{\bm{k}}}$
$\displaystyle=\frac{m}{\hbar}\bra{\frac{\partial}{\partial
k_{\mu}}u_{\bm{k}}}\frac{\partial H(\bm{k})}{\partial
k_{\nu}}\ket{u_{\bm{k}}}$
$\displaystyle=\frac{m}{\hbar}\left[\bra{\frac{\partial}{\partial
k_{\mu}}u}\frac{\partial}{\partial
k_{\nu}}\left(H(\bm{k})\ket{u}\right)-\bra{\frac{\partial}{\partial
k_{\mu}}u}H(\bm{k})\ket{\frac{\partial}{\partial k_{\nu}}u}\right]$
$\displaystyle=\frac{m}{\hbar}\left[\bra{\frac{\partial}{\partial
k_{\mu}}u}\frac{\partial}{\partial
k_{\nu}}\left(\varepsilon_{{}_{\bm{k}}}\ket{u}\right)-\bra{\frac{\partial}{\partial
k_{\mu}}u}H(\bm{k})\ket{\frac{\partial}{\partial k_{\nu}}u}\right]$
$\displaystyle=\frac{m}{\hbar}\left[\bra{\frac{\partial}{\partial
k_{\mu}}u}\left(\varepsilon_{{}_{\bm{k}}}\ket{\frac{\partial}{\partial
k_{\nu}}u}\right)+\bra{\frac{\partial}{\partial
k_{\mu}}u}\left(\frac{\partial\varepsilon_{{}_{\bm{k}}}}{\partial
k_{\nu}}\ket{u}\right)-\varepsilon_{{}_{\bm{k}}}\bra{\frac{\partial}{\partial
k_{\mu}}u}\ket{\frac{\partial}{\partial k_{\nu}}u}\right]$
$\displaystyle=\frac{m}{\hbar}\left[\bra{\frac{\partial}{\partial
k_{\mu}}u}\varepsilon_{{}_{\bm{k}}}-H(\bm{k})\ket{\frac{\partial}{\partial
k_{\nu}}u}\right]+\frac{m}{\hbar}\frac{\partial\varepsilon_{{}_{\bm{k}}}}{\partial
k_{\nu}}\bra{\frac{\partial}{\partial k_{\mu}}u}\ket{u}$
$\displaystyle=\frac{m}{\hbar}\left[\bra{\frac{\partial}{\partial
k_{\mu}}u}\varepsilon_{{}_{\bm{k}}}-H(\bm{k})\ket{\frac{\partial}{\partial
k_{\nu}}u}\right]+i{A}_{\mu}(\bm{k})\expectationvalue{{\tilde{P}}_{\nu}(\bm{k})}$

####  C.1.3 Calculation of $\mathcal{S}={}_{\text{all
space}}\bra{u_{\bm{k}_{1}}}e^{-i\bm{k}\cdot(\hat{\bm{r}}-\bm{r}_{0})}(\bm{r}-\bm{r}_{0})_{\mu}\hat{\bm{P}}_{\nu}e^{i\bm{k}\cdot(\bm{r}-\bm{r}_{0})}\ket{u_{\bm{k}_{2}}}_{\text{all
space}}$

Notation: ${}_{\text{all space}}\bra{\psi_{1}}\ket{\psi_{2}}_{\text{all
space}}=\int_{\text{all space}}d\bm{r}\psi_{1}^{*}\psi_{2}$ and
${}_{\text{cell}}\bra{\psi_{1}}\ket{\psi_{2}}_{\text{cell}}=\int_{\text{unit
cell}}d\bm{r}\psi_{1}^{*}\psi_{2}$. If $\psi_{1}$, $\psi_{2}$ are periodic
over unit cells, then ${}_{\text{all
space}}\bra{\psi_{1}}\ket{\psi_{2}}_{\text{all
space}}=N{}_{\text{cell}}\bra{\psi_{1}}\ket{\psi_{2}}_{\text{cell}}$, where
$N$ is the number of unit cells.

We use the relation,
$\hat{\bm{P}}(e^{i\bm{k}\cdot(\bm{r}-\bm{r}_{0})}\ket{u_{\bm{k}_{2}}})=e^{i\bm{k}\cdot(\bm{r}-\bm{r}_{0})}\tilde{\bm{P}}(\bm{k}_{2})\ket{u_{\bm{k}_{2}}}$.
Then,

$\displaystyle\mathcal{S}=$ $\displaystyle{}_{\text{all
space}}\bra{u_{\bm{k}_{1},n}}e^{-i\bm{k}\cdot(\hat{\bm{r}}-\bm{r}_{0})}(\bm{r}-\bm{r}_{0})_{\mu}\hat{\bm{P}}_{\nu}e^{i\bm{k}\cdot(\bm{r}-\bm{r}_{0})}\ket{u_{\bm{k}_{2},n}}_{\text{all
space}}$ $\displaystyle={}_{\text{all
space}}\bra{u_{\bm{k}_{1},n}}e^{i(\bm{k}_{2}-\bm{k}_{1})\cdot(\hat{\bm{r}}-\bm{r}_{0})}(\bm{r}-\bm{r}_{0})_{\mu}\tilde{\bm{P}}_{\nu}(\bm{k}_{2})\ket{u_{\bm{k}_{2},n}}_{\text{all
space}}$

Since $\ket{u_{\bm{k}_{2},n}}$ for different values of $n$ are eigenstates of
$H(\bm{k}_{2})$ for fixed value of $\bm{k}_{2}$, they form a complete set.
Then, $\mathcal{I}=\sum_{n}\ket{u_{\bm{k}_{2},n}}{}_{\text{cell
}}{}_{\text{cell}}\bra{u_{\bm{k}_{2},n}}$, and due to periodicity,
$\frac{1}{N}\sum_{n}\ket{u_{\bm{k}_{2},n}}{}_{\text{all space }}{}_{\text{all
space}}\bra{u_{\bm{k}_{2},n}}=\mathcal{I}$. We insert this in the expression
above.

$\displaystyle\mathcal{S}$ $\displaystyle={}_{\text{all
space}}\bra{u_{\bm{k}_{1},n}}e^{i(\bm{k}_{2}-\bm{k}_{1})\cdot(\hat{\bm{r}}-\bm{r}_{0})}(\bm{r}-\bm{r}_{0})_{\mu}\tilde{\bm{P}}_{\nu}(\bm{k}_{2})\ket{u_{\bm{k}_{2},n}}_{\text{all
space}}$ $\displaystyle={}_{\text{all
space}}\bra{u_{\bm{k}_{1},n}}i\left(\frac{\partial}{\partial{k_{1}}_{\mu}}e^{i(\bm{k}_{2}-\bm{k}_{1})\cdot(\bm{r}-\bm{r}_{0})}\right)\tilde{\bm{P}}_{\nu}(\bm{k}_{2})\ket{u_{\bm{k}_{2},n}}_{\text{all
space}}$ $\displaystyle=\frac{1}{N}\sum_{n^{\prime}}{}_{\text{all
space}}\bra{u_{\bm{k}_{1},n}}i\left(\frac{\partial}{\partial{k_{1}}_{\mu}}e^{i(\bm{k}_{2}-\bm{k}_{1})\cdot(\bm{r}-\bm{r}_{0})}\right)\ket{u_{\bm{k}_{2},n^{\prime}}}{}_{\text{all
space }}{}_{\text{all
space}}\bra{u_{\bm{k}_{2},n^{\prime}}}\tilde{\bm{P}}_{\nu}(\bm{k}_{2})\ket{u_{\bm{k}_{2},n}}_{\text{all
space}}$

where we used the periodicity of the Bloch wavefunctions to absorb the
normalization factor. Now we would manipulate the derivative with respect to
${\bm{k}_{1}}_{\mu}$.

$\displaystyle\mathcal{S}$
$\displaystyle=\frac{i}{N}\sum_{n^{\prime}}\left(\frac{\partial}{\partial{k_{1}}_{\mu}}{}_{\text{all
space}}\bra{u_{\bm{k}_{1},n}}e^{i(\bm{k}_{2}-\bm{k}_{1})\cdot(\bm{r}-\bm{r}_{0})}\right)\ket{u_{\bm{k}_{2},n^{\prime}}}_{\text{all
space}}{}_{\text{all
space}}\bra{u_{\bm{k}_{2},n^{\prime}}}\tilde{\bm{P}}_{\nu}(\bm{k}_{2})\ket{u_{\bm{k}_{2},n}}{}_{\text{all
space}}$ $\displaystyle-\frac{i}{N}\sum_{n^{\prime}}{}_{\text{all
space}}\bra{\frac{\partial}{\partial{k_{1}}_{\mu}}u_{\bm{k}_{1},n}}e^{i(\bm{k}_{2}-\bm{k}_{1})\cdot(\bm{r}-\bm{r}_{0})}\ket{u_{\bm{k}_{2},n^{\prime}}}{}_{\text{all
space }}{}_{\text{all
space}}\bra{u_{\bm{k}_{2},n^{\prime}}}\tilde{\bm{P}}_{\nu}(\bm{k}_{2})\ket{u_{\bm{k}_{2},n}}{}_{\text{all
space}}$

In the first term, we use the relation $\frac{1}{N}{}_{\text{all
space}}\bra{u_{\bm{k}_{2},n^{\prime}}}\tilde{\bm{P}}_{\nu}(\bm{k}_{2})\ket{u_{\bm{k}_{2},n}}{}_{\text{all
space}}={}_{\text{cell}}\bra{u_{\bm{k}_{2},n^{\prime}}}\tilde{\bm{P}}_{\nu}(\bm{k}_{2})\ket{u_{\bm{k}_{2},n}}{}_{\text{cell}}$
and in the second term, we again use
$\frac{1}{N}\sum_{n}\ket{u_{\bm{k}_{2},n}}{}_{\text{all space }}{}_{\text{all
space}}\bra{u_{\bm{k}_{2},n}}=\mathcal{I}$.

$\displaystyle\mathcal{S}$
$\displaystyle=i\sum_{n^{\prime}}\frac{\partial}{\partial{k_{1}}_{\mu}}\left({}_{\text{all
space}}\bra{u_{\bm{k}_{1},n}}e^{i(\bm{k}_{2}-\bm{k}_{1})\cdot(\bm{r}-\bm{r}_{0})}\ket{u_{\bm{k}_{2},n^{\prime}}}_{\text{all
space}}\right){}_{\text{cell}}\bra{u_{\bm{k}_{2},n^{\prime}}}\tilde{\bm{P}}_{\nu}(\bm{k}_{2})\ket{u_{\bm{k}_{2},n}}{}_{\text{cell
}}$ $\displaystyle-i{}_{\text{ all
space}}\bra{\frac{\partial}{\partial{k_{1}}_{\mu}}u_{\bm{k}_{1},n}}e^{i(\bm{k}_{2}-\bm{k}_{1})\cdot(\bm{r}-\bm{r}_{0})}\tilde{\bm{P}}_{\nu}(\bm{k}_{2})\ket{u_{\bm{k}_{2},n}}{}_{\text{all
space}}$

When the number of sites in a lattice is very large, the points in reciprocal
space are dense, and we can approximate the sum over the real space lattice
sites,
$\sum_{\bm{R}}e^{i(\bm{k}_{1}-\bm{k}_{2})\cdot\bm{R}}\rightarrow\frac{1}{V_{c}}\int
d\bm{r}e^{i(\bm{k}_{1}-\bm{k}_{2})\cdot\bm{r}}=\frac{(2\pi)^{d}}{V_{c}}\delta(\bm{k}_{1}-\bm{k}_{2})$,
where $V_{c}$ is the volume of the $d$ dimensional unit cell. We can further
write $\frac{(2\pi)^{d}}{V_{c}}=V_{\text{BZ}}$, where $V_{\text{BZ}}$ is the
volume of the ${1^{\textit{st}}}$ Brillouin zone.

For any periodic function over the unit cells,

$\int_{\text{all
space}}e^{i(\bm{k}_{1}-\bm{k}_{2})\cdot\bm{r}}f(\text{periodic})=\left(\sum_{\bm{R}\in\text{
lattice
sites}}e^{i(\bm{k}_{1}-\bm{k}_{2})\cdot\bm{R}}\right)\left(\int_{\text{unit
cell}}f(\text{periodic})\right)=\left[V_{\text{BZ}}\delta(\bm{k}_{1}-\bm{k}_{2})\right]\left(\int_{\text{unit
cell}}f(\text{periodic})\right).$

Now we choose units such that $V_{\text{BZ}}=1$ for the simplicity of the
subsequent calculations (anyway it can always be absorbed in the
normalization). We get,

$\displaystyle\mathcal{S}$
$\displaystyle=i\sum_{n^{\prime}}\frac{\partial}{\partial{k_{1}}_{\mu}}\left(\delta(\bm{k}_{1}-\bm{k}_{2}){}_{\text{cell}}\bra{u_{\bm{k}_{1},n}}\ket{u_{\bm{k}_{2},n^{\prime}}}_{\text{cell}}\right){}_{\text{cell}}\bra{u_{\bm{k}_{2},n^{\prime}}}\tilde{\bm{P}}_{\nu}(\bm{k}_{2})\ket{u_{\bm{k}_{2},n}}{}_{\text{cell
}}$ $\displaystyle-i\delta(\bm{k}_{1}-\bm{k}_{2}){}_{\text{
cell}}\bra{\frac{\partial}{\partial{k_{1}}_{\mu}}u_{\bm{k}_{1},n}}\tilde{\bm{P}}_{\nu}(\bm{k}_{2})\ket{u_{\bm{k}_{2},n}}{}_{\text{cell}}$
$\displaystyle=i\sum_{n^{\prime}}\frac{\partial}{\partial{k_{1}}_{\mu}}\left(\delta(\bm{k}_{1}-\bm{k}_{2})\right){}_{\text{cell}}\bra{u_{\bm{k}_{1},n}}\ket{u_{\bm{k}_{2},n^{\prime}}}_{\text{cell}}{}_{\text{cell}}\bra{u_{\bm{k}_{2},n^{\prime}}}\tilde{\bm{P}}_{\nu}(\bm{k}_{2})\ket{u_{\bm{k}_{2},n}}{}_{\text{cell
}}$
$\displaystyle+i\sum_{n^{\prime}}\delta(\bm{k}_{1}-\bm{k}_{2}){}_{\text{cell}}\bra{\frac{\partial}{\partial{k_{1}}_{\mu}}u_{\bm{k}_{1},n}}\ket{u_{\bm{k}_{2},n^{\prime}}}_{\text{cell}}{}_{\text{cell}}\bra{u_{\bm{k}_{2},n^{\prime}}}\tilde{\bm{P}}_{\nu}(\bm{k}_{2})\ket{u_{\bm{k}_{2},n}}{}_{\text{cell
}}$ $\displaystyle-i\delta(\bm{k}_{1}-\bm{k}_{2}){}_{\text{
cell}}\bra{\frac{\partial}{\partial{k_{1}}_{\mu}}u_{\bm{k}_{1},n}}\tilde{\bm{P}}_{\nu}(\bm{k}_{2})\ket{u_{\bm{k}_{2},n}}{}_{\text{cell}}$

After summing over $n^{\prime}$ in the second term of the above expression,
the second and the third terms cancel each other.

$\displaystyle\mathcal{S}$
$\displaystyle=i\sum_{n^{\prime}}\left[\frac{\partial}{\partial{k_{1}}_{\mu}}\delta(\bm{k}_{1}-\bm{k}_{2})\right]{}_{\text{cell}}\bra{u_{\bm{k}_{1},n}}\ket{u_{\bm{k}_{2},n^{\prime}}}_{\text{cell}}{}_{\text{cell}}\bra{u_{\bm{k}_{2},n^{\prime}}}\tilde{\bm{P}}_{\nu}(\bm{k}_{2})\ket{u_{\bm{k}_{2},n}}{}_{\text{cell
}}$
$\displaystyle=\boxed{i\left[\frac{\partial}{\partial{k_{1}}_{\mu}}\delta(\bm{k}_{1}-\bm{k}_{2})\right]{}_{\text{cell}}\bra{u_{\bm{k}_{1},n}}\tilde{\bm{P}}_{\nu}(\bm{k}_{2})\ket{u_{\bm{k}_{2},n}}{}_{\text{cell
}}}$

####  C.1.4 To show that
$\bra{\psi_{\bm{k},\bm{r}_{0}}}\\{(\hat{\bm{r}}-\bm{r}_{0})_{\mu},\hat{P}_{\nu}\\}\ket{\psi_{\bm{k},\bm{r}_{0}}}$
is antisymmetric in $\mu$, $\nu$ for a Bloch wavepacket
$\psi_{\bm{k},\bm{r}_{0}}$, which implies $\mathcal{M}_{ij}$ is totally
antisymmetric

It is mentioned in Ref.[32] that
$\mathcal{M}_{\mu\nu}=\bra{\psi_{\bm{k},\bm{r}_{0}}}\frac{1}{2}\\{(\bm{r}-\bm{r}_{0})_{\mu},\hat{P}_{\nu}\\}\ket{\psi_{\bm{k},\bm{r}_{0}}}$
is a completely antisymmetric tensor. Here we prove it. A similar calculation
to find the orbital angular momentum of a wavepacket can be found at Appendix
B of Ref.[4]. First we calculate the quantity without the anticommutator.

$\displaystyle\bra{\psi_{\bm{k},\bm{r}_{0}}}(\bm{r}-\bm{r}_{0})_{\mu}\hat{P}_{\nu}\ket{\psi_{\bm{k},\bm{r}_{0}}}$
$\displaystyle=\int\int
d\bm{k}_{1}d\bm{k}_{2}w(\bm{k}_{1}-\bm{k})w(\bm{k}_{2}-\bm{k})e^{-i\bm{A}(\bm{k})\cdot(\bm{k}_{1}-\bm{k})}e^{i\bm{A}(\bm{k})\cdot(\bm{k}_{2}-\bm{k})}$
$\displaystyle\times{}_{\text{all
space}}\bra{u_{\bm{k}_{1}}}e^{-i\bm{k}\cdot(\hat{\bm{r}}-\bm{r}_{0})}(\bm{r}-\bm{r}_{0})_{\mu}\hat{\bm{P}}_{\nu}e^{i\bm{k}\cdot(\bm{r}-\bm{r}_{0})}\ket{u_{\bm{k}_{2}}}_{\text{all
space}}$

Now we use the boxed result from the previous section,

$\displaystyle=\int\int
d\bm{k}_{1}d\bm{k}_{2}w(\bm{k}_{1}-\bm{k})w(\bm{k}_{2}-\bm{k})e^{-i\bm{A}(\bm{k})\cdot(\bm{k}_{1}-\bm{k})}e^{i\bm{A}(\bm{k})\cdot(\bm{k}_{2}-\bm{k})}$
$\displaystyle\times
i\left[\frac{\partial}{\partial{k_{1}}_{\mu}}\delta(\bm{k}_{1}-\bm{k}_{2})\right]{}_{\text{cell}}\bra{u_{\bm{k}_{1},n}}\tilde{\bm{P}}_{\nu}(\bm{k}_{2})\ket{u_{\bm{k}_{2},n}}{}_{\text{cell
}}$

After integrating by parts to shift the derivative from the delta function to
the other quantities,

$\displaystyle=\int\int
d\bm{k}_{1}d\bm{k}_{2}w(\bm{k}_{1}-\bm{k})w(\bm{k}_{2}-\bm{k}_{0})e^{-i\bm{A}(\bm{k})\cdot(\bm{k}_{1}-\bm{k})}e^{i\bm{A}(\bm{k})\cdot(\bm{k}_{2}-\bm{k})}$
$\displaystyle\times
i\left[\frac{\partial}{\partial{k_{1}}_{\mu}}\delta(\bm{k}_{1}-\bm{k}_{2})\right]{}_{\text{cell}}\bra{u_{\bm{k}_{1},n}}\tilde{\bm{P}}_{\nu}(\bm{k}_{2})\ket{u_{\bm{k}_{2},n}}{}_{\text{cell
}}$ $\displaystyle=-i\int\int
d\bm{k}_{1}d\bm{k}_{2}w(\bm{k}_{1}-\bm{k})w(\bm{k}_{2}-\bm{k})e^{-i\bm{A}(\bm{k}_{0})\cdot(\bm{k}_{1}-\bm{k}_{0})}e^{i\bm{A}(\bm{k})\cdot(\bm{k}_{2}-\bm{k})}$
$\displaystyle\times\delta(\bm{k}_{1}-\bm{k}_{2})\left[{}_{\text{cell}}\bra{\frac{\partial}{\partial{k_{1}}_{\mu}}u_{\bm{k}_{1},n}}\tilde{\bm{P}}_{\nu}(\bm{k}_{2})\ket{u_{\bm{k}_{2},n}}{}_{\text{cell
}}-i\bm{A}_{\mu}(\bm{k}){}_{\text{
cell}}\bra{u_{\bm{k}_{1},n}}\tilde{\bm{P}}_{\nu}(\bm{k}_{2})\ket{u_{\bm{k}_{2},n}}{}_{\text{cell
}}\right]$

Integrating over the delta function, and using the result
${}_{\text{cell}}\bra{\frac{\partial}{\partial{k_{1}}_{\mu}}u_{\bm{k}_{1},n}}\tilde{\bm{P}}_{\nu}(\bm{k}_{1})\ket{u_{\bm{k}_{1},n}}{}_{\text{cell
}}=\frac{m}{\hbar}\left[\bra{\frac{\partial}{\partial{k_{1}}_{\mu}}u}\varepsilon_{{}_{\bm{k_{1}}}}-H(\bm{k_{1}})\ket{\frac{\partial}{\partial{k_{1}}_{\nu}}u}\right]+i{A}_{\mu}(\bm{k}_{1})\expectationvalue{{P}_{\nu}(\bm{k}_{1})}$
from the previous section, we get,

$\displaystyle-i\int
d\bm{k}_{1}[w(\bm{k}_{1}-\bm{k})]^{2}\frac{m}{\hbar}\left[\bra{\frac{\partial}{\partial{k_{1}}_{\mu}}u}\varepsilon_{{}_{\bm{k}_{1}}}-H(\bm{k}_{1})\ket{\frac{\partial}{\partial{k_{1}}_{\nu}}u}\right]$
$\displaystyle\boxed{=-i\frac{m}{\hbar}\left[\bra{\frac{\partial}{\partial
k_{\mu}}u}\varepsilon_{{}_{\bm{k}}}-H(\bm{k})\ket{\frac{\partial}{\partial
k_{\nu}}u}\right]}$

The last line follows because $w(\bm{k}_{1}-\bm{k})$ being sharply peaked at
$\bm{k}_{1}=\bm{k}$, picks up the value of the integrand at that point, and
due to normalization of $\psi(\bm{r})_{n,\bm{k}}$ in Eq. (28), $\int
d\bm{k}_{1}[w(\bm{k}_{1}-\bm{k})]^{2}=1$ (In other words,
$[w(\bm{k}_{1}-\bm{k})]^{2}$ effectively behaves like
$\delta{(\bm{k}_{1}-\bm{k})}$).

Case: $\mu\neq\nu$.

In this case, $(\bm{r}-\bm{r}_{0})_{\mu}\hat{P}_{\nu}$ Hermitian (because
$(\bm{r}-\bm{r}_{0})_{\mu}$ and $\hat{P}_{\nu}$ commute, and they are
individually Hermitian), and its expectation value must be real. Then,
$\bra{\frac{\partial}{\partial
k_{\mu}}u}\varepsilon_{{}_{\bm{k}}}-H(\bm{k})\ket{\frac{\partial}{\partial
k_{\nu}}u}=\frac{i\hbar}{m}\bra{\psi_{\bm{k},\bm{r}_{0}}}\\{(\bm{r}-\bm{r}_{0})_{\mu}\hat{P}_{\nu}\\}\ket{\psi_{\bm{k},\bm{r}_{0}}}$
must be purely imaginary.

Now, $\varepsilon_{{}_{\bm{k}}}-H(\bm{k})$ is Hermitian. Then,
$\bra{\frac{\partial}{\partial
k_{\mu}}u}\varepsilon_{{}_{\bm{k}}}-H(\bm{k})\ket{\frac{\partial}{\partial
k_{\nu}}u}^{*}=\bra{\frac{\partial}{\partial
k_{\nu}}u}\varepsilon_{{}_{\bm{k}}}-H(\bm{k})\ket{\frac{\partial}{\partial
k_{\mu}}u}$.

But since this is purely imaginary, $\bra{\frac{\partial}{\partial
k_{\mu}}u}\varepsilon_{{}_{\bm{k}}}-H(\bm{k})\ket{\frac{\partial}{\partial
k_{\nu}}u}^{*}=-\bra{\frac{\partial}{\partial
k_{\mu}}u}\varepsilon_{{}_{\bm{k}}}-H(\bm{k})\ket{\frac{\partial}{\partial
k_{\nu}}u}$.

Combining the two, we get, $\bra{\frac{\partial}{\partial
k_{\mu}}u}\varepsilon_{{}_{\bm{k}}}-H(\bm{k})\ket{\frac{\partial}{\partial
k_{\nu}}u}=-\bra{\frac{\partial}{\partial
k_{\nu}}u}\varepsilon_{{}_{\bm{k}}}-H(\bm{k})\ket{\frac{\partial}{\partial
k_{\mu}}u}$.

Also, since $(\bm{r}-\bm{r}_{0})_{\mu}$ and $\hat{P}_{\nu}$ commute,
$(\bm{r}-\bm{r}_{0})_{\mu}\hat{P}_{\nu}=\frac{1}{2}\\{(\bm{r}-\bm{r}_{0})_{\mu},\hat{P}_{\nu}\\}$.

Consequently, for $\mu\neq\nu$,
$\bra{\psi_{\bm{k},\bm{r}_{0}}}\frac{1}{2}\\{(\bm{r}-\bm{r}_{0})_{\mu},\hat{P}_{\nu}\\}\ket{\psi_{\bm{k},\bm{r}_{0}}}=-\bra{\psi_{\bm{k},\bm{r}_{0}}}\frac{1}{2}\\{(\bm{r}-\bm{r}_{0})_{\nu},\hat{P}_{\mu}\\}\ket{\psi_{\bm{k},\bm{r}_{0}}}$,
i.e., $\boxed{\mathcal{M}_{\mu\nu}=-\mathcal{M}_{\nu\mu}}$.

Case: $\mu=\nu$. In this subsection, no sum is implied for repeated indices

In this case, $(\bm{r}-\bm{r}_{0})_{\mu}\hat{P}_{\mu}$ is not anymore
Hermitian. However, $\left[\bra{\frac{\partial}{\partial
k_{\mu}}u}\varepsilon_{{}_{\bm{k}}}-H(\bm{k})\ket{\frac{\partial}{\partial
k_{\mu}}u}\right]$ must be real. Consequently,
$-i\frac{m}{\hbar}\left[\bra{\frac{\partial}{\partial
k_{\mu}}u}\varepsilon_{{}_{\bm{k}}}-H(\bm{k})\ket{\frac{\partial}{\partial
k_{\mu}}u}\right]=\expectationvalue{(\bm{r}-\bm{r}_{0})_{\mu}\hat{P}_{\mu}}$
is purely imaginary, say, $iy$. Then,
$\expectationvalue{\hat{P}_{\mu}(\bm{r}-\bm{r}_{0})_{\mu}}=-iy$, and

$\expectationvalue{(\bm{r}-\bm{r}_{0})_{\mu}\hat{P}_{\mu}+\hat{P}_{\mu}(\bm{r}-\bm{r}_{0})_{\mu}}=0.$

Therefore, $\boxed{\mathcal{M}_{\mu\mu}=0}$ for all $\mu$.

Note: Since we know that
$\expectationvalue{(\bm{r}-\bm{r}_{0})_{\mu}\hat{P}_{\mu}}$ is purely
imaginary, and also,
$[(\bm{r}-\bm{r}_{0})_{\mu}\hat{P}_{\mu}-\hat{P}_{\mu}(\bm{r}-\bm{r}_{0})_{\mu}]=i\hbar$,
it follows that
$\expectationvalue{(\bm{r}-\bm{r}_{0})_{\mu}\hat{P}_{\mu}}=\frac{i\hbar}{2}$.

Therefore, $\mathcal{M}$ is a completely antisymmetric tensor.

###  C.2 $\tilde{m}_{i}=\frac{1}{2}\mathcal{M}_{jk}\epsilon_{ijk}$ is
identical to the magnetic moment defined in Eq. (10)

Since only the $j\neq k$ terms contribute in $\mathcal{M}_{jk}\epsilon_{ijk}$,
and $\hat{p}_{k}$ commutes with $\hat{r}_{j}$ for $j\neq k$, we can get rid of
the anticommutator,

$\displaystyle\tilde{m}_{i}=\frac{1}{2}\bra{\psi_{\bm{k},\bm{r}_{0}}}\frac{1}{2}\\{-e\frac{\hat{p}_{k}}{m},(\bm{\hat{r}}-\bm{r}_{0})_{j}\\}\ket{\psi_{\bm{k},\bm{r}_{0}}}\epsilon_{ijk}$
$\displaystyle=-\frac{e}{2m}\epsilon_{ijk}\bra{\psi_{\bm{k},\bm{r}_{0}}}(\bm{\hat{r}}-\bm{r}_{0})_{j}\hat{p}_{k}\ket{\psi_{\bm{k},\bm{r}_{0}}}$
$\displaystyle=-\frac{e}{2m}\bra{\psi_{\bm{k},\bm{r}_{0}}}\left((\bm{\hat{r}}-\bm{r}_{0})\times\hat{\bm{p}}\right)_{i}\ket{\psi_{\bm{k},\bm{r}_{0}}}$

which is the $i$-th component of magnetic moment $\bm{m}_{\bm{k}}$, defined in
Eq. (10), that is, $\tilde{m}_{i}=m_{i}$.

##  Appendix D Derivation of the transport electric current density and
derivation of the charge magnetization

The total electric current is,

$\displaystyle\bm{j}^{e}_{\text{total }}$
$\displaystyle=-e\int\frac{2d\bm{k}}{(2\pi)^{d}}\left[{g}_{\bm{k}}\frac{1}{\hbar}\frac{\partial\varepsilon_{0}({{\bm{k}}})}{\partial\bm{k}}+{f}_{\bm{k}}\frac{e}{\hbar}(\bm{E}\times\bm{\Omega}(\bm{k}))\right]$
(32)
$\displaystyle+\bm{\nabla}\times\int\frac{2d\bm{k}}{(2\pi)^{d}}f_{\bm{k}}\bm{m}_{\bm{k}}$
$\displaystyle=-e\int\frac{2d\bm{k}}{(2\pi)^{d}}\left[{g}_{\bm{k}}\frac{1}{\hbar}\frac{\partial\varepsilon_{0}({{\bm{k}}})}{\partial\bm{k}}+{f}_{\bm{k}}\frac{e}{\hbar}(\bm{E}\times\bm{\Omega}(\bm{k}))\right]$
$\displaystyle+\gradient{\mu}\times\int\frac{2d\bm{k}}{(2\pi)^{d}}\frac{\partial
f_{\bm{k}}}{\partial\mu}\bm{m}_{\bm{k}}+\gradient{T}\times\int\frac{2d\bm{k}}{(2\pi)^{d}}\frac{\partial
f_{\bm{k}}}{\partial T}\bm{m}_{\bm{k}}$

The electric magnetization current is,

$\displaystyle\bm{j}^{e}_{M}=\curl{\bm{M}}^{e}=\gradient{\mu}\times\frac{\partial\bm{M}^{e}}{\partial\mu}+\gradient{T}\times\frac{\partial\bm{M}^{e}}{\partial
T}$ (33)

The electric transport current is,

$\displaystyle\bm{j}^{e}_{\text{transport }}$
$\displaystyle=\bm{j}^{e}_{\text{total }}-\bm{j}^{e}_{M}$ (34)
$\displaystyle=-e\int\frac{2d\bm{k}}{(2\pi)^{d}}\left[{g}_{\bm{k}}\frac{1}{\hbar}\frac{\partial\varepsilon_{0}({{\bm{k}}})}{\partial\bm{k}}+{f}_{\bm{k}}\frac{e}{\hbar}(\bm{E}\times\bm{\Omega}(\bm{k}))\right]$
$\displaystyle+\gradient{\mu}\times\int\frac{2d\bm{k}}{(2\pi)^{d}}\frac{\partial
f_{\bm{k}}}{\partial\mu}\bm{m}_{\bm{k}}+\gradient{T}\times\int\frac{2d\bm{k}}{(2\pi)^{d}}\frac{\partial
f_{\bm{k}}}{\partial T}\bm{m}_{\bm{k}}$
$\displaystyle-\gradient{\mu}\times\frac{\partial\bm{M}^{e}}{\partial\mu}-\gradient{T}\times\frac{\partial\bm{M}^{e}}{\partial
T}$

It is shown in section VI that $g_{\bm{k}}$ is already a function of
$(e\bm{E}+\gradient{\mu})$. Then, for the Einstein relation to hold for the
electric transport current,

$\frac{\partial\bm{M}^{e}}{\partial\mu}=\int\frac{2d\bm{k}}{(2\pi)^{d}}\frac{\partial
f_{\bm{k}}}{\partial\mu}\bm{m}_{\bm{k}}+\int\frac{2d\bm{k}}{(2\pi)^{d}}{f}_{\bm{k}}\frac{e}{\hbar}\bm{\Omega}(\bm{k}).$

Integrating,

$\bm{M}^{e}=\int\frac{2d\bm{k}}{(2\pi)^{d}}\left[\bm{m}(\bm{k})f_{\bm{k}}+k_{B}T\frac{e\bm{\Omega}}{\hbar}\log(1+e^{-\beta(\varepsilon_{0}(\bm{k})-\mu)})+\bm{c}(\bm{k},T)\right],$

where $\bm{c}(\bm{k},T)$ is a $\mu$-independent function. Now, the
magnetization should be zero for an empty band, i.e., at any temperature,
$\bm{M}^{e}\rightarrow 0$ as $\mu\rightarrow-\infty$. Applying this condition,
we get $\bm{c}(\bm{k},T)=0$.

##  Appendix E Derivation of transport heat current density

Since the energy magnetization $\bm{M}_{0}^{E}$ is a function of the
temperature $T$ and the chemical potential $\mu$, it follows from Eq. (15)
that the magnetization energy current is,

$\displaystyle\bm{j}^{E}_{M}$
$\displaystyle=\gradient{\mu}\times\frac{\partial\bm{M}_{0}^{E}}{\partial\mu}+\gradient{T}\times\frac{\partial\bm{M}_{0}^{E}}{\partial
T}-\bm{E}\times\bm{M}^{e}$
$\displaystyle=\gradient{\mu}\times\frac{\partial\bm{M}_{0}^{E}}{\partial\mu}+\gradient{T}\times\frac{\partial\bm{M}_{0}^{E}}{\partial
T}$
$\displaystyle-\bm{E}\times\int\frac{2d\bm{k}}{(2\pi)^{d}}\left[\bm{m}(\bm{k})f_{\bm{k}}+k_{B}T\frac{e\bm{\Omega}}{\hbar}\log(1+e^{-\beta(\varepsilon_{0}(\bm{k})-\mu)})\right]$

Then, the transport energy current density is,

$\displaystyle\bm{j}_{\text{transport }}^{E}$
$\displaystyle=\bm{j}_{\text{total }}^{E}-\bm{j}^{E}_{M}$ (35)
$\displaystyle=\int\frac{2d\bm{k}}{(2\pi)^{d}}\varepsilon_{0}(\bm{k})\left[{g}_{\bm{k}}\frac{1}{\hbar}\frac{\partial\varepsilon_{0}(\bm{k})}{\partial\bm{k}}+\frac{e}{\hbar}{f}_{\bm{k}}(\bm{E}\times\bm{\Omega}(\bm{k}))\right]$
$\displaystyle+\bm{E}\times\int\frac{2d\bm{k}}{(2\pi)^{d}}k_{B}T\frac{e\bm{\Omega}}{\hbar}\log(1+e^{-\beta(\varepsilon_{0}(\bm{k})-\mu)})$
$\displaystyle-\gradient{\mu}\times\frac{\partial\bm{M}_{0}^{E}}{\partial\mu}-\gradient{T}\times\frac{\partial\bm{M}_{0}^{E}}{\partial
T}$

The charge magnetization has been shown to be [48, 55],
$\bm{M}^{e}=\frac{1}{\beta}\int\frac{2d\bm{k}}{(2\pi)^{d}}\frac{e}{\hbar}\bm{\Omega}(\bm{k})\log(1+e^{-\beta(\varepsilon_{0}(\bm{k})-\mu)})+\int\frac{2d\bm{k}}{(2\pi)^{d}}f_{\bm{k}}\bm{m}_{\bm{k}}$,
and the corresponding circulating magnetization electric current is,

$\displaystyle\bm{j}^{e}_{M}$ $\displaystyle=\curl{\bm{M}^{e}}$
$\displaystyle=\bm{\nabla}\times\int\frac{2d\bm{k}}{(2\pi)^{d}}f\bm{m}_{\bm{k}}+\gradient{\mu}\times\int\frac{2d\bm{k}}{(2\pi)^{d}}\frac{e}{\hbar}\bm{\Omega}(\bm{k})f+\frac{\gradient{T}}{T}\times\int\frac{2d\bm{k}}{(2\pi)^{d}}\frac{e}{\hbar}\bm{\Omega}(\bm{k})f_{\bm{k}}(\varepsilon_{0}(\bm{k})-\mu)$
$\displaystyle+\frac{\gradient{T}}{T}\times\int\frac{2d\bm{k}}{(2\pi)^{d}}\frac{e}{\hbar}\bm{\Omega}(\bm{k})k_{B}T\log(1+e^{-\beta(\varepsilon_{0}(\bm{k})-\mu)})$

The electric transport current is [48],

$\displaystyle{{\bm{j}}^{e}}_{\text{transport}}$
$\displaystyle=\int\frac{2d\bm{k}}{(2\pi)^{d}}(-e){g}_{\bm{k}}\frac{1}{\hbar}\frac{\partial\varepsilon_{0}({{\bm{k}}})}{\partial\bm{k}}-{f}_{\bm{k}}\frac{e^{2}}{\hbar}\left(\left[\bm{E}+\frac{\gradient{\mu}}{e}\right]\times\bm{\Omega}(\bm{k})\right)$
$\displaystyle-\frac{\gradient{T}}{T}\times\int\frac{2d\bm{k}}{(2\pi)^{d}}\frac{e}{\hbar}\bm{\Omega}(\bm{k})\left[f_{\bm{k}}(\varepsilon_{0}(\bm{k})-\mu)+k_{B}T\log(1+e^{-\beta(\varepsilon_{0}(\bm{k})-\mu)})\right]$

The number density transport current is
$\bm{j}^{N}_{\text{transport}}=\frac{\bm{j}^{e}_{\text{transport}}}{-e}$, and
the full expression for the transport heat current density (Eq. (23)) is,

$\displaystyle\bm{j}_{\text{transport }}^{Q}$
$\displaystyle=\bm{j}_{\text{transport }}^{E}-\mu\bm{j}_{\text{transport
}}^{N}$
$\displaystyle=\int\frac{2d\bm{k}}{(2\pi)^{d}}(\varepsilon_{0}(\bm{k})-\mu)\frac{1}{\hbar}\frac{\partial\varepsilon_{0}(\bm{k})}{\partial\bm{k}}g_{k}$
$\displaystyle+\bm{E}\times\int\frac{2d\bm{k}}{(2\pi)^{d}}(\varepsilon_{0}(\bm{k})-\mu)\frac{e}{\hbar}\bm{\Omega}f_{\bm{k}}$
$\displaystyle+\bm{E}\times\int\frac{2d\bm{k}}{(2\pi)^{d}}k_{B}T\frac{e\bm{\Omega}}{\hbar}\log(1+e^{-\beta(\varepsilon_{0}(\bm{k})-\mu)})$
$\displaystyle+\int\frac{2d\bm{k}}{(2\pi)^{d}}(-\mu)\frac{1}{\hbar}\gradient{\mu}\times\bm{\Omega}f_{\bm{k}}$
$\displaystyle-\frac{\mu}{\hbar}\int\frac{2d\bm{k}}{(2\pi)^{d}}\frac{\gradient{T}}{T}\times\bm{\Omega}f_{k}\left(\varepsilon_{0}(\bm{k})_{k}-\mu\right)$
$\displaystyle-\frac{\mu}{\hbar}\int\frac{2d\bm{k}}{(2\pi)^{d}}\frac{\gradient{T}}{T}\times\bm{\Omega}k_{B}T\log(1+e^{-\beta(\varepsilon_{0}(\bm{k})-\mu)})$
$\displaystyle-\gradient{\mu}\times\frac{\partial\bm{M}_{0}^{E}}{\partial\mu}-\gradient{T}\times\frac{\partial\bm{M}_{0}^{E}}{\partial
T}$

##  Appendix F Simplification of the condition on energy magnetization at low
temperature

At very low temperature (in the limit $\beta\mu\gg 1$), for $|\bm{k}|<k_{f}$,
we have $f_{\bm{k}}\approx 1$, and
$\log(1+e^{-\beta(\varepsilon_{0}(\bm{k})-\mu)})\approx\log(e^{-\beta(\varepsilon_{0}(\bm{k})-\mu)})=\beta(\mu-\varepsilon_{0}(\bm{k}))$.
Then,
$\varepsilon_{0}(\bm{k})f_{\bm{k}}+k_{B}T\log(1+e^{-\beta(\varepsilon_{0}(\bm{k})-\mu)})\approx\mu$.
And for $|\bm{k}|>k_{f}$, $f_{\bm{k}}\approx
e^{-\beta(\varepsilon_{0}(\bm{k})-\mu)}\ll 1$, and
$\log(1+e^{-\beta(\varepsilon_{0}(\bm{k})-\mu)})\approx
e^{-\beta(\varepsilon_{0}(\bm{k})-\mu)}\ll 1$. Then, the quantity
$\left[\varepsilon_{0}(\bm{k})f_{\bm{k}}+k_{B}T\log(1+e^{-\beta(\varepsilon_{0}(\bm{k})-\mu)})\right]$
is negligibly small compared to $\mu$, and it can be neglected.
Combining both the relations, we can write

$\varepsilon_{0}(\bm{k})f_{\bm{k}}+k_{B}T\log(1+e^{-\beta(\varepsilon_{0}(\bm{k})-\mu)})\approx\mu\Theta(\mu-\varepsilon_{0}(\bm{k}))\approx\mu
f_{\bm{k}},$

where $\Theta$ is the Heaviside step function.

##  Appendix G Discussion about the absolute value of the chemical potential
in the expression of the circulating energy current at the edge

Note that the condition on the energy magnetization (Eq. (19)) depends on the
absolute value of the chemical potential (not something like its relative
value with respect to the minima of the band). However, this is not
unphysical. When we change the chemical potential by a small amount (i.e.
$|\frac{\Delta\mu}{\mu-\varepsilon_{\text{min}}}|\ll 1$), new electrons are
added to the system, and they have energy in the order of $\mu$. Even if these
newly added electrons belong to a bulk state (i.e., not a chiral edge state),
the circulating energy currents due to these bulk states (which cancel in the
bulk) would contribute to the change in the energy current circulating in the
edges, and this change is in the order of (the absolute value of) $\mu$, and
the net contribution is, as if, only chiral edge states are contributing to
the circulating current (A similar result holds for the circulating electric
current [56]). For example, if the chemical potential is at zero energy (then
the band minima would be at some negative energy), and new electrons are added
to the system, they would have zero energy, and the circulating energy current
at the edges would not change, a result which agrees with Eq. (20).

##  Appendix H Calculations for the Boltzmann Transport Equation

In this appendix, the subscripts $\bm{k}$ of $\varepsilon$, $g$ and $\tau$ are
not explicitly written. It is to be understood that they are functions of the
crystal momentum.

###  H.1 Simplification of RHS of Eq. (25)

Since $f=\frac{1}{e^{\beta(\varepsilon-\mu)}+1}$, it follows that
$\frac{\partial f}{\partial\bm{r}}=\frac{\partial
f}{\partial\varepsilon}\left[-\frac{\gradient{T}}{T}(\varepsilon-\mu)-\gradient{\mu}\right]$,
and $\frac{\partial f}{\partial\bm{k}}=\frac{\partial
f}{\partial\varepsilon}\frac{\partial\varepsilon}{\partial\bm{k}}$. Using
${\dot{\bm{r}}=\frac{\frac{1}{\hbar}\frac{\partial\varepsilon}{\partial\bm{k}}+\frac{e}{\hbar}(\bm{E}\times\bm{\Omega})}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}}$,
and
$\dot{\bm{k}}=-\frac{\frac{e}{\hbar}\bm{E}+\frac{e}{\hbar^{2}}\frac{\partial\varepsilon}{\partial\bm{k}}\times\bm{B}}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}$,
we now calculate
$\dot{\bm{r}}\cdot\frac{\partial}{\partial\bm{r}}f+\dot{\bm{k}}\cdot\frac{\partial}{\partial\bm{k}}f$.
Here, $\frac{\partial f}{\partial\bm{r}}$ is linear in external fields, and
$\dot{\bm{r}}$ has a term $\bm{E}\times\bm{\Omega}$, and we neglect the
product of these terms, which is a second order quantity.

Then we get,
$-\dot{\bm{r}}\cdot\frac{\partial}{\partial\bm{r}}f-\dot{\bm{k}}\cdot\frac{\partial}{\partial\bm{k}}f=\frac{1}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\frac{1}{\hbar}\frac{\partial\varepsilon}{\partial\bm{k}}\cdot\left[e\bm{E}+\gradient{\mu}+\frac{\gradient{T}}{T}(\varepsilon-\mu)\right]$.

###  H.2 Solution of the BTE for $\bm{E}\neq 0$,
$\gradient{\mu}=\gradient{T}=0$

In this case, the equation becomes,

$~{}\frac{g}{\tau}-\frac{\frac{e}{\hbar^{2}}}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\frac{\partial\varepsilon}{\partial\bm{k}}\times\bm{B}\cdot\frac{\partial}{\partial\bm{k}}g=\frac{\partial
f}{\partial\varepsilon}\frac{1}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\frac{1}{\hbar}\frac{\partial\varepsilon}{\partial\bm{k}}\cdot
e\bm{E}$ (36)

We solve this equation, treating the second term in LHS as a perturbation (See
Appendix I for justification). We write $g=g_{0}+g_{1}$, such that
$\frac{g_{0}}{\tau}=\frac{\partial
f}{\partial\varepsilon}\frac{1}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\frac{1}{\hbar}\frac{\partial\varepsilon}{\partial\bm{k}}\cdot
e\bm{E}$, and
$\frac{g_{1}}{\tau}-\frac{\frac{e}{\hbar^{2}}}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\frac{\partial\varepsilon}{\partial\bm{k}}\times\bm{B}\cdot\frac{\partial}{\partial\bm{k}}g_{0}=0$.
Then, $g_{1}$ is a linear function of $\bm{B}$, and it is justified to discard
the term
$\frac{\frac{e}{\hbar^{2}}}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\frac{\partial\varepsilon}{\partial\bm{k}}\times\bm{B}\cdot\frac{\partial}{\partial\bm{k}}g_{1}$,
as it would be quadratic in $\bm{B}$.

Then, ${g_{0}}=\frac{\partial
f}{\partial\varepsilon}\frac{1}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\frac{e\tau}{\hbar}\frac{\partial\varepsilon}{\partial\bm{k}}\cdot\bm{E}$.
To find $g_{1}$, we need to calculate $\frac{\partial g_{0}}{\partial\bm{k}}$,
i.e., the quantity $\frac{\partial}{\partial\bm{k}}\left[\frac{\partial
f}{\partial\varepsilon}\frac{1}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\frac{e\tau}{\hbar}\frac{\partial\varepsilon}{\partial\bm{k}}\cdot\bm{E}\right]$.

To calculate it, let us first calculate
$\gradient{\left[\phi\bm{A}\cdot\bm{C}\right]}$, where $\phi$ is a scalar
function, $\bm{A}$ is a vector function, and $\bm{C}$ is a constant vector.

$\gradient{\left[\phi\bm{A}\cdot\bm{C}\right]}=(\gradient{\phi})(\bm{A}\cdot\bm{C})+\phi(\bm{C}\cdot\bm{\nabla}){A}+\phi(\bm{C}\times(\curl{\bm{A}}))$.
In our calculation, $\phi\sim\frac{\partial
f}{\partial\varepsilon}\frac{1}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\frac{e\tau}{\hbar}$,
$\bm{A}\sim\frac{\partial\varepsilon}{\partial\bm{k}}$, and
$\bm{C}\sim\bm{E}$.

Then,

$\frac{\partial}{\partial\bm{k}}\left[\frac{\partial
f}{\partial\varepsilon}\frac{1}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\frac{e\tau}{\hbar}\frac{\partial\varepsilon}{\partial\bm{k}}\cdot\bm{E}\right]=\frac{\partial}{\partial\bm{k}}\left[\frac{\frac{\partial
f}{\partial\varepsilon}\tau}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\right]\frac{e\bm{E}}{\hbar}\cdot\frac{\partial\varepsilon}{\partial\bm{k}}+\frac{\frac{\partial
f}{\partial\varepsilon}\tau}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\left[\left(\frac{e}{\hbar}\bm{E}\cdot\frac{\partial}{\partial\bm{k}}\right)\frac{\partial\varepsilon}{\partial\bm{k}}+\frac{e}{\hbar}\bm{E}\times\left(\frac{\partial}{\partial\bm{k}}\times\frac{\partial\varepsilon}{\partial\bm{k}}\right)\right].$

The last term is zero because it is the curl of a gradient,
$\frac{\partial}{\partial\bm{k}}\times\frac{\partial\varepsilon}{\partial\bm{k}}=\nabla_{{}_{\bm{k}}}\times(\nabla_{{}_{\bm{k}}}\varepsilon)=0$.

Finally, the solution is, obtained by adding $g_{0}$ and $g_{1}$,

$g=\frac{\partial
f}{\partial\varepsilon}\frac{1}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\frac{\tau}{\hbar}\frac{\partial\varepsilon}{\partial\bm{k}}\cdot
e\bm{E}+\frac{\frac{e\tau}{\hbar^{2}}}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\frac{\partial\varepsilon}{\partial\bm{k}}\cdot\bm{B}\times\left[\frac{\partial}{\partial\bm{k}}\left[\frac{\frac{\partial
f}{\partial\varepsilon}\tau}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\right]\left(\frac{e\bm{E}}{\hbar}\cdot\frac{\partial\varepsilon}{\partial\bm{k}}\right)+\frac{\frac{\partial
f}{\partial\varepsilon}\tau}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\left[\left(\frac{e}{\hbar}\bm{E}\cdot\frac{\partial}{\partial\bm{k}}\right)\frac{\partial\varepsilon}{\partial\bm{k}}\right]\right]$
(37)

###  H.3 Solution of the BTE for $\bm{E}=0$, $\gradient{\mu}\neq 0$,
$\gradient{T}\neq 0$

Under these conditions, the equation becomes,

$\frac{g}{\tau}-\frac{\frac{e}{\hbar^{2}}\frac{\partial\varepsilon}{\partial\bm{k}}\times\bm{B}}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\cdot\frac{\partial}{\partial\bm{k}}g=\frac{\partial
f}{\partial\varepsilon}\frac{1}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\frac{1}{\hbar}\frac{\partial\varepsilon}{\partial\bm{k}}\cdot\left[\gradient{\mu}+\gradient{T}\frac{\varepsilon-\mu}{T}\right]$

Let $g=g_{0}+g_{1}$, with $g_{0}=\frac{\partial
f}{\partial\varepsilon}\frac{1}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\frac{\tau}{\hbar}\frac{\partial\varepsilon}{\partial\bm{k}}\cdot\left[\gradient{\mu}+\gradient{T}\frac{\varepsilon-\mu}{T}\right]$
and
$\frac{g_{1}}{\tau}-\frac{\frac{e}{\hbar^{2}}\frac{\partial\varepsilon}{\partial\bm{k}}\times\bm{B}}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\cdot\frac{\partial}{\partial\bm{k}}g_{0}=0$.

The form is otherwise similar to the form of the equations in Appendix H.2,
but here we have a factor of $(\varepsilon-\mu)$ with $\gradient{T}$. As a
result, we would get additional factors like
$\frac{\partial\varepsilon}{\partial\bm{k}}\cdot\frac{\gradient{T}}{T}$ when
we calculate $\frac{\partial g_{0}}{\partial\bm{k}}$ in order to calculate
$g_{1}$. One might think that such additional factors may violate the Onsager
relation. But such factors would cancel due to a vector triple product being
zero, and the Onsager relation continues to hold.

We have,

$\frac{\partial
g_{0}}{\partial\bm{k}}=\frac{\partial}{\partial\bm{k}}\left[\frac{\frac{\partial
f}{\partial\varepsilon}\tau}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\right]\frac{\gradient{\mu}+\frac{\varepsilon-\mu}{T}}{\hbar}\cdot\frac{\partial\varepsilon}{\partial\bm{k}}+\frac{\frac{\partial
f}{\partial\varepsilon}\tau}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\left[\left(\frac{\gradient{\mu}+\frac{\varepsilon-\mu}{T}}{\hbar}\cdot\frac{\partial}{\partial\bm{k}}\right)\frac{\partial\varepsilon}{\partial\bm{k}}+\frac{1}{\hbar}\left(\gradient{\mu}+\frac{\varepsilon-\mu}{T}\right)\times\underbrace{\left(\frac{\partial}{\partial\bm{k}}\times\frac{\partial\varepsilon}{\partial\bm{k}}\right)}_{0}\right]$
$+\frac{\frac{\partial
f}{\partial\varepsilon}\tau}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\left[\frac{1}{\hbar}\left(\frac{\partial\varepsilon}{\partial\bm{k}}\cdot\frac{\partial}{\partial\bm{k}}\right)\left(\frac{\varepsilon\gradient{T}}{T}\right)+\frac{1}{\hbar}\frac{\partial\varepsilon}{\partial\bm{k}}\times\left(\frac{\partial}{\partial\bm{k}}\times\left(\frac{\varepsilon\gradient{T}}{T}\right)\right)\right]$

The term inside the last third bracket can be further simplified.

$\left(\frac{\partial\varepsilon}{\partial\bm{k}}\cdot\frac{\partial}{\partial\bm{k}}\right)\left(\frac{\varepsilon\gradient{T}}{T}\right)+\frac{\partial\varepsilon}{\partial\bm{k}}\times\left(\frac{\partial}{\partial\bm{k}}\times\left(\frac{\varepsilon\gradient{T}}{T}\right)\right)=\left(\frac{\partial\varepsilon}{\partial\bm{k}}\cdot\frac{\partial\varepsilon}{\partial\bm{k}}\right)\left(\frac{\gradient{T}}{T}\right)+\frac{\partial\varepsilon}{\partial\bm{k}}\times\left(\frac{\partial\varepsilon}{\partial\bm{k}}\times\left(\frac{\gradient{T}}{T}\right)\right)$
$=\left(\frac{\partial\varepsilon}{\partial\bm{k}}\cdot\frac{\partial\varepsilon}{\partial\bm{k}}\right)\left(\frac{\gradient{T}}{T}\right)+\left(\frac{\partial\varepsilon}{\partial\bm{k}}\cdot\frac{\gradient{T}}{T}\right)\frac{\partial\varepsilon}{\partial\bm{k}}-(\frac{\partial\varepsilon}{\partial\bm{k}}\cdot\frac{\partial\varepsilon}{\partial\bm{k}})(\frac{\gradient{T}}{T})$
$=\left(\frac{\partial\varepsilon}{\partial\bm{k}}\cdot\frac{\gradient{T}}{T}\right)\frac{\partial\varepsilon}{\partial\bm{k}}.$

When we calculate
$g_{1}=\tau\frac{\frac{e}{\hbar^{2}}\frac{\partial\varepsilon}{\partial\bm{k}}\times\bm{B}}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\cdot\frac{\partial}{\partial\bm{k}}g_{0}$,
this term cancels because
$(\frac{\partial\varepsilon}{\partial\bm{k}}\times\bm{B})\cdot\frac{\partial\varepsilon}{\partial\bm{k}}=0$.

Then, the full solution is, $g=g_{0}+g_{1}$

$\displaystyle=\frac{\partial
f}{\partial\varepsilon}\frac{1}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\frac{\tau}{\hbar}\frac{\partial\varepsilon}{\partial\bm{k}}\cdot\left[\gradient{\mu}+\gradient{T}\frac{\varepsilon-\mu}{T}\right]$
$\displaystyle+\frac{\frac{e\tau}{\hbar^{2}}}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\frac{\partial\varepsilon}{\partial\bm{k}}\cdot\bm{B}\times\left[\frac{\partial}{\partial\bm{k}}\left[\frac{\frac{\partial
f}{\partial\varepsilon}\tau}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\right]\frac{\gradient{\mu}+\gradient{T}\frac{\varepsilon-\mu}{T}}{\hbar}\cdot\frac{\partial\varepsilon}{\partial\bm{k}}+\frac{\frac{\partial
f}{\partial\varepsilon}\tau}{1+\frac{e}{\hbar}\bm{B}\cdot\bm{\Omega}}\left[\left(\frac{\gradient{\mu}+\gradient{T}\frac{\varepsilon-\mu}{T}}{\hbar}\cdot\frac{\partial}{\partial\bm{k}}\right)\frac{\partial\varepsilon}{\partial\bm{k}}\right]\right]$

Since Eq. (26) is linear, when the electric field, the temperature gradient
and the chemical potential gradient are each non-zero, we can add the two
solutions in appendix H.2 and appendix H.3, and we get the solution in Eq.
(27).

##  Appendix I Validity of the perturbation theory on Boltzmann Transport
Equation

Let us do an order of the magnitude calculation of the term containing the
magnetic field. Classically,
$\dot{\bm{k}}\cdot\frac{\partial}{\partial\bm{k}}g\sim\dot{\bm{p}}\cdot\frac{\partial}{\partial\bm{p}}g\sim
e(\bm{v}\times{B})\cdot\frac{\partial}{m\partial\bm{v}}g\sim\bm{\omega}\times\bm{v}\cdot\frac{\partial
g}{\partial\bm{v}}\sim\omega g$.

Then, if $\omega\ll\frac{1}{\tau}$ (in the limit of low magnetic field), it is
justified to treat $\omega g$ as a perturbation over $\frac{g}{\tau}$.

###  I.1 How low magnetic is ‘low’?

For $\omega=\frac{1}{\tau}$, $B_{crit}=\frac{m^{*}}{\tau e}$, and a magnetic
field much lower than this qualifies as a “low” magnetic field. If we take
$m^{*}=m_{e}$ (the free electron mass) and $\tau\sim 10^{-14}s$ (chapter 3 of
Ref.[62]), $B_{crit}\sim 570T$ (which is found in a typical metal), which is
order of magnitudes above the magnetic field accessible in labs. In typical
metals, $m^{*}>m_{e}$, the bare electron mass, and the critical magnetic field
is even higher.

In very clean samples, $\tau$ can be as high as $10^{-12}s$, and even there,
the critical magnetic field would be $5.7T$, which too is hard to achieve in
most labs. In less clean samples, the critical magnetic field is much more.

Conclusion: Unless the sample is very clean, the magnetic fields produced in
labs qualify as “low”, and in this regime, this perturbation theory remains
valid.

## References

  * Berry [1984] M. V. Berry, Quantal phase factors accompanying adiabatic changes, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 392, 45 (1984).
  * Karplus and Luttinger [1954] R. Karplus and J. M. Luttinger, Hall effect in ferromagnetics, Phys. Rev. 95, 1154 (1954).
  * Chang and Niu [1995] M.-C. Chang and Q. Niu, Berry phase, hyperorbits, and the hofstadter spectrum, Phys. Rev. Lett. 75, 1348 (1995).
  * Chang and Niu [1996] M.-C. Chang and Q. Niu, Berry phase, hyperorbits, and the hofstadter spectrum: Semiclassical dynamics in magnetic bloch bands, Phys. Rev. B 53, 7010 (1996).
  * Sundaram and Niu [1999] G. Sundaram and Q. Niu, Wave-packet dynamics in slowly perturbed crystals: Gradient corrections and berry-phase effects, Phys. Rev. B 59, 14915 (1999).
  * Chong [2010] Y. D. Chong, Berry’s phase and the anomalous velocity of bloch wavepackets, Phys. Rev. B 81, 052303 (2010).
  * Ralph [2020] D. C. Ralph, Berry curvature, semiclassical electron dynamics, and topological materials: Lecture notes for introduction to solid state physics (2020), arXiv:2001.04797 [cond-mat.mes-hall] .
  * 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).
  * Mandal  _et al._ [2020] D. Mandal, K. Das, and A. Agarwal, Magnus nernst and thermal hall effect, Phys. Rev. B 102, 205414 (2020).
  * Xiao  _et al._ [2021] R.-C. Xiao, Z. Wang, Z.-Q. Zhang, J. Liu, and H. Jiang, Magnus hall effect in two-dimensional materials, Chinese Physics Letters 38, 057301 (2021).
  * Xiao  _et al._ [2010] D. Xiao, M.-C. Chang, and Q. Niu, Berry phase effects on electronic properties, Rev. Mod. Phys. 82, 1959 (2010).
  * Behnia and Aubin [2016] K. Behnia and H. Aubin, Nernst effect in metals and superconductors: a review of concepts and experiments, Reports on Progress in Physics 79, 046502 (2016).
  * Yang  _et al._ [2020] H. Yang, W. You, J. Wang, J. Huang, C. Xi, X. Xu, C. Cao, M. Tian, Z.-A. Xu, J. Dai, and Y. Li, Giant anomalous nernst effect in the magnetic weyl semimetal ${\mathrm{co}}_{3}{\mathrm{sn}}_{2}{\mathrm{s}}_{2}$, Phys. Rev. Materials 4, 024202 (2020).
  * Saha and Tewari [2018] S. Saha and S. Tewari, Anomalous nernst effect in type-ii weyl semimetals, The European Physical Journal B 91, 4 (2018).
  * Ivanov  _et al._ [2018] Y. V. Ivanov, A. T. Burkov, and D. A. Pshenay-Severin, Thermoelectric properties of topological insulators, physica status solidi (b) 255, 1800020 (2018).
  * Caio  _et al._ [2016] M. D. Caio, N. R. Cooper, and M. J. Bhaseen, Hall response and edge current dynamics in chern insulators out of equilibrium, Phys. Rev. B 94, 155104 (2016).
  * Ma  _et al._ [2019] D. Ma, H. Jiang, H. Liu, and X. C. Xie, Planar hall effect in tilted weyl semimetals, Phys. Rev. B 99, 115121 (2019).
  * Dantas  _et al._ [2018] R. M. A. Dantas, F. Peña-Benitez, B. Roy, and P. Surówka, Magnetotransport in multi-weyl semimetals: a kinetic theory approach, Journal of High Energy Physics 2018, 69 (2018).
  * Burkov [2014] A. A. Burkov, Anomalous hall effect in weyl metals, Phys. Rev. Lett. 113, 187202 (2014).
  * Das and Agarwal [2019a] K. Das and A. Agarwal, Linear magnetochiral transport in tilted type-i and type-ii weyl semimetals, Phys. Rev. B 99, 085405 (2019a).
  * Das and Agarwal [2019b] K. Das and A. Agarwal, Berry curvature induced thermopower in type-i and type-ii weyl semimetals, Phys. Rev. B 100, 085406 (2019b).
  * Lu and Shen [2017] H.-Z. Lu and S.-Q. Shen, Quantum transport in topological semimetals under magnetic fields, Frontiers of Physics 12, 127201 (2017).
  * König  _et al._ [2016] E. J. König, P. M. Ostrovsky, M. Dzero, and A. Levchenko, Anomalous hall effect on the surface of topological kondo insulators, Phys. Rev. B 94, 041403 (2016).
  * Xiao  _et al._ [2016] C. Xiao, D. Li, and Z. Ma, Thermoelectric response of spin polarization in rashba spintronic systems, Frontiers of Physics 11, 117201 (2016).
  * Kapri  _et al._ [2021] P. Kapri, B. Dey, and T. K. Ghosh, Role of berry curvature in the generation of spin currents in rashba systems, Phys. Rev. B 103, 165401 (2021).
  * Price and Cooper [2012] H. M. Price and N. R. Cooper, Mapping the berry curvature from semiclassical dynamics in optical lattices, Phys. Rev. A 85, 033620 (2012).
  * Spurrier and Cooper [2018] S. Spurrier and N. R. Cooper, Semiclassical dynamics, berry curvature, and spiral holonomy in optical quasicrystals, Phys. Rev. A 97, 043603 (2018).
  * Wang  _et al._ [2021] Z. Wang, L. Dong, C. Xiao, and Q. Niu, Berry curvature effects on quasiparticle dynamics in superconductors, Phys. Rev. Lett. 126, 187001 (2021).
  * Höckendorf  _et al._ [2020] B. Höckendorf, A. Alvermann, and H. Fehske, Topological origin of quantized transport in non-hermitian floquet chains, Phys. Rev. Research 2, 023235 (2020).
  * Silberstein  _et al._ [2020] N. Silberstein, J. Behrends, M. Goldstein, and R. Ilan, Berry connection induced anomalous wave-packet dynamics in non-hermitian systems, Phys. Rev. B 102, 245147 (2020).
  * Xiao  _et al._ [2012] D. Xiao, G.-B. Liu, W. Feng, X. Xu, and W. Yao, Coupled spin and valley physics in monolayers of ${\mathrm{mos}}_{2}$ and other group-vi dichalcogenides, Phys. Rev. Lett. 108, 196802 (2012).
  * Dong  _et al._ [2020] L. Dong, C. Xiao, B. Xiong, and Q. Niu, Berry phase effects in dipole density and the mott relation, Phys. Rev. Lett. 124, 066601 (2020).
  * Xiao and Niu [2017] C. Xiao and Q. Niu, Semiclassical theory of spin-orbit torques in disordered multiband electron systems, Phys. Rev. B 96, 045428 (2017).
  * Friedlan and Dignam [2021] A. Friedlan and M. M. Dignam, Valley polarization in biased bilayer graphene using circularly polarized light, Phys. Rev. B 103, 075414 (2021).
  * Zhang  _et al._ [2009] C. Zhang, S. Tewari, and S. Das Sarma, Berry-phase-mediated topological thermoelectric transport in gapped single and bilayer graphene, Phys. Rev. B 79, 245424 (2009).
  * Xiao  _et al._ [2007] D. Xiao, W. Yao, and Q. Niu, Valley-contrasting physics in graphene: Magnetic moment and topological transport, Phys. Rev. Lett. 99, 236809 (2007).
  * Ostrovsky  _et al._ [2008] P. M. Ostrovsky, I. V. Gornyi, and A. D. Mirlin, Theory of anomalous quantum hall effects in graphene, Phys. Rev. B 77, 195430 (2008).
  * Sodemann and Fu [2015] I. Sodemann and L. Fu, Quantum nonlinear hall effect induced by berry curvature dipole in time-reversal invariant materials, Phys. Rev. Lett. 115, 216806 (2015).
  * Gao and Xiao [2018] Y. Gao and D. Xiao, Orbital magnetic quadrupole moment and nonlinear anomalous thermoelectric transport, Phys. Rev. B 98, 060402 (2018).
  * Zhang  _et al._ [2018] Y. Zhang, J. van den Brink, C. Felser, and B. Yan, Electrically tuneable nonlinear anomalous hall effect in two-dimensional transition-metal dichalcogenides WTe 2 and MoTe 2, 2D Mater. 5, 044001 (2018).
  * Xiao  _et al._ [2019a] C. Xiao, H. Zhou, and Q. Niu, Scaling parameters in anomalous and nonlinear hall effects depend on temperature, Phys. Rev. B 100, 161403 (2019a).
  * Xiao  _et al._ [2019b] C. Xiao, Z. Z. Du, and Q. Niu, Theory of nonlinear hall effects: Modified semiclassics from quantum kinetics, Phys. Rev. B 100, 165422 (2019b).
  * Yu  _et al._ [2019] X.-Q. Yu, Z.-G. Zhu, J.-S. You, T. Low, and G. Su, Topological nonlinear anomalous nernst effect in strained transition metal dichalcogenides, Phys. Rev. B 99, 201410 (2019).
  * Ashcroft and Mermin [1976] N. W. Ashcroft and N. D. Mermin,  _Solid State Physics_ (Holt-Saunders, 1976).
  * Onsager [1931a] L. Onsager, Reciprocal relations in irreversible processes. i., Phys. Rev. 37, 405 (1931a).
  * Onsager [1931b] L. Onsager, Reciprocal relations in irreversible processes. ii., Phys. Rev. 38, 2265 (1931b).
  * 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).
  * Xiao  _et al._ [2006] D. Xiao, Y. Yao, Z. Fang, and Q. Niu, Berry-phase effect in anomalous thermoelectric transport, Phys. Rev. Lett. 97, 026603 (2006).
  * Jackson [1999] J. D. Jackson,  _Classical electrodynamics_, 3rd ed. (Wiley, New York, NY, 1999).
  * Kim  _et al._ [2014] K.-S. Kim, H.-J. Kim, and M. Sasaki, Boltzmann equation approach to anomalous transport in a weyl metal, Phys. Rev. B 89, 195137 (2014).
  * Pal  _et al._ [2021] O. Pal, B. Dey, and T. K. Ghosh, Berry curvature induced magnetotransport in 3d noncentrosymmetric metals, J. Phys.: Condens. Matter 34, 025702 (2021).
  * Xiao  _et al._ [2005] D. Xiao, J. Shi, and Q. Niu, Berry phase correction to electron density of states in solids, Phys. Rev. Lett. 95, 137204 (2005).
  * Thonhauser  _et al._ [2005] T. Thonhauser, D. Ceresoli, D. Vanderbilt, and R. Resta, Orbital magnetization in periodic insulators, Phys. Rev. Lett. 95, 137205 (2005).
  * Ceresoli  _et al._ [2006] D. Ceresoli, T. Thonhauser, D. Vanderbilt, and R. Resta, Orbital magnetization in crystalline solids: Multi-band insulators, chern insulators, and metals, Phys. Rev. B 74, 024408 (2006).
  * Shi  _et al._ [2007] J. Shi, G. Vignale, D. Xiao, and Q. Niu, Quantum theory of orbital magnetization and its generalization to interacting systems, Phys. Rev. Lett. 99, 197202 (2007).
  * Resta [2010] R. Resta, Electrical polarization and orbital magnetization: the modern theories, Journal of Physics: Condensed Matter 22, 123201 (2010).
  * Hatsugai [1993] Y. Hatsugai, Chern number and edge states in the integer quantum hall effect, Phys. Rev. Lett. 71, 3697 (1993).
  * Mong and Shivamoggi [2011] R. S. K. Mong and V. Shivamoggi, Edge states and the bulk-boundary correspondence in dirac hamiltonians, Phys. Rev. B 83, 125109 (2011).
  * Mukerjee and Huse [2004] S. Mukerjee and D. A. Huse, Nernst effect in the vortex-liquid regime of a type-ii superconductor, Phys. Rev. B 70, 014506 (2004).
  * Ziman [1972] J. M. Ziman,  _Principles of the theory of solids_ , 2nd ed. (CUP, 1972).
  * Culcer  _et al._ [2004] D. Culcer, J. Sinova, N. A. Sinitsyn, T. Jungwirth, A. H. MacDonald, and Q. Niu, Semiclassical spin transport in spin-orbit-coupled bands, Phys. Rev. Lett. 93, 046602 (2004).
  * Simon [2013] S. H. Simon,  _The Oxford Solid State Basics_ (Oxford University Press, 2013).