# Quantum kinetic equation and thermal conductivity tensor for bosons

Léo Mangeolle  École Normale Supérieure de Lyon, CNRS, Laboratoire de
physique, 46, allée d’Italie, 69007 Lyon  Kavli Institute for Theoretical
Physics, University of California, Santa Barbara, CA 93106-4030  Lucile Savary
École Normale Supérieure de Lyon, CNRS, Laboratoire de physique, 46, allée
d’Italie, 69007 Lyon  Kavli Institute for Theoretical Physics, University of
California, Santa Barbara, CA 93106-4030  Leon Balents  Kavli Institute for
Theoretical Physics, University of California, Santa Barbara, CA 93106-4030
Canadian Institute for Advanced Research, Toronto, Ontario, Canada

###### Abstract

We systematically derive the quantum kinetic equation in full phase space for
any quadratic hamiltonian of bosonic fields, including in the absence of
translational invariance. This enables the treatment of boundaries,
inhomogeneous systems and states with non-trivial textures, such as skyrmions
in the context of magnetic bosons. We relate the evolution of the distribution
of bosons in phase space to single-electron, band-diagonal, physical
quantities such as Berry curvature and energy magnetization, by providing a
procedure to “diagonalize” the Hamiltonian in phase space, using the formalism
of the Moyal product. We obtain exact equations, which can be expanded order
by order, for example in the “smallness” of the spatial gradients, providing a
“semiclassical” approximation. In turn, at first order, we recover the usual
full Boltzmann equation and give a self-contained and exact derivation of the
intrinsic thermal Hall effect of bosons. The formulation clarifies the
contribution from “energy magnetization” in natural manner, and does not
require the inclusion of Luttinger’s pseudo-gravitational field to obtain
thermal transport quantities.

##  I Introduction

It is now widely recognized that the dynamics of electrons in energy bands is
fundamentally affected by topological and geometric aspects of the Bloch
eigenstates. This is described semi-classically by Berry phases and
curvatures, which enter the equations of motion for a wave packet Sundaram and
Niu (1999). The same ideas apply to any smooth collection of non-degenerate
single particle states indexed by continuous parameters, and in particular to
bosonic elementary excitations like phonons or magnons Qin  _et al._ (2012);
Matsumoto  _et al._ (2014). These concepts are particularly powerful in
describing out of, but near equilibrium phenomena such as charge and energy
transport: incorporating the semi-classical topological dynamics into a
Boltzmann description provides the most compact and intuitive understanding of
the anomalous Hall effect, as well as many other transport properties Karplus
and Luttinger (1954); Adams and Blount (1959); Haldane (2004); Nagaosa (2006);
Sinitsyn (2007); Mishchenko and Starykh (2014). Being semi-classical, i.e.
allowing simultaneous consideration of position and momentum, such a method
also makes it possible to treat inhomogeneous systems, including boundaries,
textures, etc., in a compact and natural manner. It sometimes offers
advantages even for calculation of bulk transport properties like thermal
conductivity, because it allows to transparently isolate spurious
contributions due to bound “magnetization” currents, which have in the past
obscured correct results Katsura  _et al._ (2010).

The semi-classical approach is known to reproduce the exact result for the
intrinsic anomalous Hall effect Haldane (2004), as is often the case for
transport coefficients. This is because the semi-classical approximation is
controlled by the smallness of spatial and temporal gradients – of the
Hamiltonian and of the distribution function – which are indeed small in a
perturbative response to uniform DC fields. Theoretically, the smallness of
gradients is required both to justify the semi-classical equations of motion,
which describe the motion of one particle in phase space, the Boltzmann
equation itself, which describes the evolution of the distribution of
particles, and which is an approximation to the fully quantum evolution of the
density matrix. The derivation of the Boltzmann equation semi-classically from
the full quantum kinetic equation (QKE) is an old problem studied extensively
prior to the widespread incorporation of Berry phase effects into band
dynamics Keldysh (1964); Al’Tshuler (1978); Rammer and Smith (1986). Despite
the common practice of combining these two semi-classical approximations, they
are almost always treated independently, with attention focusing primarily on
the single particle equations of motion and the Boltzmann equation adopted
without justification.

A unified derivation for the electronic case was provided by Wickles and
Belzig Gosselin  _et al._ (2007); Wickles and Belzig (2013). They showed that
the equation of motion for the “one-particle density matrix” of a multi-band
electron system can be systematically reduced to the Boltzmann one, with the
renormalized single particle equations of motion emerging in the same
treatment via a single semi-classical expansion. Here we extend this treatment
to bosons, which, contrary to electrons, do not have a conserved
number/charge, but only energy. The difference of statistics, but more
importantly the latter lack of conserved charge, lead to some significant
differences from the Wickles and Belzig treatment. Nevertheless, we obtain a
full derivation of the leading semi-classical kinetic equation and associated
observable quantities such as energy and current densities, and the
formulation allows a clear route to extend to higher orders in the semi-
classical expansion.

The derivation we present assumes non-interacting bosons, and hence neglects
any scattering (but this can of course be added). Aside from the neglect of
interactions, the derived quantum kinetic equation is asymptotically exact in
the limit in which the parameters of the theory are slowly varying in space.
In particular, we obtain exactly the thermal conductivity tensor, which
describes the heat current induced by a small temperature gradient. Our
approach therefore provides a self-contained and exact derivation of the
intrinsic thermal Hall effect of non-interacting bosons such as phonons and
magnons, and in the last section we show that the quantum kinetic equation
(QKE) indeed is in agreement with Kubo formula calculations Qin  _et al._
(2012); Matsumoto  _et al._ (2014). The advantage is that the QKE calculation
is considerably more intuitive, and naturally avoids subtleties associated
with magnetization currents Smrcka and Streda (1977); Cooper  _et al._ (1997);
Bashan and Auerbach (2022), which plagued early Kubo calculations Katsura  _et
al._ (2010), as mentioned above. For concreteness we further compute and plot
the local (i.e. position-resolved) energy currents for a system of chiral
phonons (excitations of linear elasticity theory with a phonon Hall viscosity
term) in a finite geometry in the presence of a temperature gradient.

Figure 1: Local momentum-integrated energy currents
${J}_{X_{y}}^{(1)}(x)\equiv\int_{p}\mathcal{J}_{X_{y}}^{(1)}$ (first term in
Eq. (17), blue) and
${J}_{X_{y}}^{(2)}(x)\equiv\int_{p}\mathcal{J}_{X_{y}}^{(2)}$ (magnetization
current, second term in Eq. (17), red) in arbitrary units for a system of
chiral phonons described by the Lagrangian Eq. (66), calculated for $g(x)={\rm
exp}\left\\{-\frac{1}{4}\Theta[|x|-L/2]\left(|x|-L/2\right)^{3}/\xi^{3}\right\\}$
(yellow, mimicking the existence of a boundary) in the presence of a constant
temperature gradient $\delta T/L$, i.e. for a temperature profile
$T(x)=T_{0}+(x/L)\delta T$, using the numerical values given in Table 1.

##  II Overview

The principal result of this paper, the derivation of the semi-classical QKE
and that of observables, is somewhat technical and subtle. Hence we provide in
this section an overview of the derivation, highlighting the key steps,
leaving full details to Sec. III.

###  II.1 Basic formulation

Our starting point is a general free boson system described by $2N$ free
hermitian fields $\Phi_{a}(r)$, $r\in\mathbb{R}^{d}$, where $d$ is the
dimension of our system, $a=1,..,2N$. For free fields the commutator

$[{\Phi}_{a}(r_{1}),{\Phi}_{b}(r_{2})]\equiv\hbar\,\hat{\sf\Gamma}_{ab}(r_{1},r_{2}),$
(1)

is a c-number $2N\times 2N$ matrix. For example, in the context of three-
dimensional elasticity, $\Phi_{1,2,3}$ might be the displacement fields
$u_{x,y,z}$, and $\Phi_{4,5,6}$ represent the conjugate momenta $\Pi_{x,y,z}$.
Since the commutator is antisymmetric and the fields are hermitian, we have
$\hat{\sf\Gamma}_{ab}(r_{1},r_{2})=-\hat{\sf\Gamma}_{ba}(r_{2},r_{1})$ and
$\hat{\sf\Gamma}_{ba}(r_{2},r_{1})=\left(\hat{\sf\Gamma}_{ab}(r_{1},r_{2})\right)^{*}$.

We take an arbitrary quadratic Hamiltonian

$\textrm{H}=\frac{1}{2}\int_{r_{1},r_{2}}\sum_{a,b}\hat{\sf
H}_{ab}(r_{1},r_{2})\Phi_{a}(r_{1})\Phi_{b}(r_{2}),$ (2)

where $\hat{\sf H}(r_{1},r_{2})$ is a c-number $2N\times 2N$ matrix. The
hermiticity of ${\rm H}$ requires $\hat{\sf H}(r_{2},r_{1})=\left(\hat{\sf
H}(r_{1},r_{2})\right)^{*}$, and to remove any redundancy we impose $\hat{\sf
H}_{ab}(r_{1},r_{2})=\hat{\sf H}_{ba}(r_{2},r_{1})$. Both $\hat{\sf H}$ and
$\hat{\sf\Gamma}$, regarded as matrices in the index ($a,b$) $\times$
coordinate space ($r_{1},r_{2}$), are hermitian matrices.

Apart from hermiticity and the symmetry/antisymmetry of
$\hat{\mathsf{H}}$/$\hat{\mathsf{\Gamma}}$, respectively, we need only assume
that $\hat{\mathsf{H}}$ has a spectrum which is bounded below, so that the
system is stable.

The expectation value of any quadratic observable can be obtained as a linear
combination of the expectation values

$\hat{\sf
F}_{ab}(r_{1},r_{2})\equiv\frac{1}{2}\left\langle\left\\{\Phi_{a}(r_{1}),\Phi_{b}(r_{2})\right\\}\right\rangle.$
(3)

Here the expectation value is the expectation value with respect to an initial
many-body density matrix $\rho_{0}$, $\langle A(t)\rangle={\rm
Tr}[A(t)\rho_{0}]$, with $A(t)$ an operator in the Heisenberg picture. The
corresponding c-number matrix $\hat{\sf F}$, which captures any two-point
expectation value of any two observables, is sometimes called ‘density
matrix,’ and is symmetric and hermitian in the index $\times$ coordinate
space.

###  II.2 Transition to phase space

From Eq. (1) and Eq. (2) the unitary time evolution of any operator is
determined. To frame this in phase space, we introduce the Wigner transform

$\displaystyle f(X,p)$ $\displaystyle=[\hat{f}]^{W}(X,p)$
$\displaystyle\equiv\int_{x}e^{-i\frac{px}{\hbar}}\hat{f}\left(X+\frac{x}{2},X-\frac{x}{2}\right),$
(4)

where $p,X,x\in\mathbb{R}^{d}$ (we can identify $X=\frac{1}{2}(r_{1}+r_{2})$
and $x=r_{1}-r_{2}$) and $\hbar$ is Planck’s constant divided by $2\pi$. We
denote $\int_{x}=\int d^{d}x$ for spatial integrals, and later we will use
$\int_{p}=\int d^{d}p/(2\pi\hbar)^{d}$ for momentum integrals (then
$\int_{x,p}$ is dimensionless). We will denote Wigner transformed matrices by
removing the hats, e.g. $[\hat{\sf F}]^{W}={\sf F}$. These objects,
$\mathsf{F},\mathsf{H},\mathsf{\Gamma}$ etc. are functions on phase space
($X,p$), and are all hermitian as matrices at fixed $X,p$. Note that we will
use the convention that $p$ denotes momentum while $k$ denotes wavenumber,
i.e. $p=\hbar k$. This is useful when we keep track of factors of $\hbar$ and
to obtain the classical limit.

In terms of Wigner transformed objects, one can then derive the equation of
motion. Namely, the equation of motion for the ‘density matrix’ in phase space
is

$\partial_{t}{\sf F}(X,p)=\frac{-i}{\hbar}\left({\sf K}\star{\sf F}-{\sf
F}\star{\sf K}^{\dagger}\right),$ (5)

where

${\sf K}(X,p)=\hbar\,{\sf\Gamma}\star{\sf H}$ (6)

is sometimes called the “dynamical matrix.” Note that $\mathsf{K}$ is not
hermitian.

Here the $\star$ represents the “star” or “Moyal” product, defined as, for any
two matrix functions ${\sf O}_{1,2}$,

${\sf O}_{1}\star{\sf O}_{2}\equiv{\sf
O}_{1}\exp\left(i\frac{\hslash}{2}\epsilon_{\alpha\beta}\reflectbox{$\vec{\reflectbox{$\partial$}}$}_{\alpha}\vec{\partial}_{\beta}\right){\sf
O}_{2},$ (7)

where repeated indices are summed, $\alpha,\beta=1,..,2d$,
$\partial_{\alpha}\equiv\partial_{q_{\alpha}}$, where $q_{\mu}=X_{\mu}$ and
$q_{d+\mu}=p_{\mu}$ for $\mu=1,..,d$, and where we have defined
$\epsilon_{\alpha\beta}$ to be such that

$\epsilon_{X_{\mu}p_{\nu}}=-\epsilon_{p_{\mu}X_{\nu}}=\delta_{\mu\nu},$ (8)

with $\mu=1,..,d$. Note that $\hslash=\hbar$ but we are formally
distinguishing the two here to treat $\hslash$ as an expansion parameter when
it stems from an expansion of Eq. (7) Wickles and Belzig (2013). We provide
more details in Appendix A. See also Refs. Rammer and Smith (1986);
Polkovnikov (2010); Wickles and Belzig (2013).

###  II.3 Semiclassical solution

A proper semi-classical limit is obtained by first formally diagonalizing the
dynamical matrix and distribution function, i.e. reducing them to scalars
within each band, and then performing a small $\hslash$ expansion, which
expands the star products Eq. (7). Diagonalization is accomplished via a
similarity transformation $\mathsf{S}$ such that
$\mathsf{S}^{-1}\star\mathsf{K}\star\mathsf{S}=\mathsf{K}_{d}$ is a diagonal
matrix of mode frequencies. Similarly,
$\mathsf{F}=\mathsf{S}\star\mathsf{F}_{d}\star\mathsf{S}^{\dagger}$, with
$\mathsf{F}_{d}$ a diagonal distribution function and the inverse ${\sf
S}^{-1}$ is defined with respect to the star product (${\sf S}^{-1}\star{\sf
S}={\sf S}\star{\sf S}^{-1}=1$, see Sec. III.2 for details). Note that
$\mathsf{F}_{d}$ is real by definition, and we show later that
$\mathsf{K}_{d}$ is also real.

There is some phase ambiguity in the diagonalization, which can be regarded as
a gauge freedom since the diagonalization exists at every point in phase
space. Physical quantities must be gauge invariant. We find that the diagonal
gauge invariant quantities $\underline{\sf O}_{d}$ take the form

$\displaystyle\underline{\sf O}_{d}(q)\equiv{\sf
O}_{d}[\underline{\smash{q}}],$ (9)

with ${\sf O}_{d}=\mathsf{K}_{d},\mathsf{F}_{d}$, and

$\displaystyle\underline{\smash{q}}_{\alpha}$
$\displaystyle=q_{\alpha}+\hslash\epsilon_{\alpha\beta}{\sf
A}_{\beta}+O(\hslash^{2}),$ (10)

where $\mathsf{A}_{\alpha}$ is a Berry gauge field:

$\displaystyle\mathsf{A}_{\alpha}$
$\displaystyle=\textrm{Im}\,\left(\mathsf{\Lambda}_{\alpha}\right)^{({\rm
d})},$ $\displaystyle\mathsf{\Lambda}_{\alpha}$
$\displaystyle=\mathsf{S}^{-1}\star\partial_{\alpha}\mathsf{S}.$ (11)

Here $({\rm d})$ denotes the diagonal part and Im the imaginary part. Note
that ${\sf A}$ is diagonal so these equations, e.g. Eq. (9), are unambiguous:
see Eq. (44) for an explicit expression for a diagonal matrix function of a
diagonal matrix.

Applying this prescription we obtain the gauge invariant kinetic equation to
first order in $\hslash$,

$\displaystyle\partial_{t}\underline{\sf F}_{d}+\upsilon\partial_{t}{\sf
q}_{\alpha}\partial_{\alpha}\underline{\sf F}_{d}=0+O(\upsilon\hslash^{2}),$
(12)

where we work in the collisionless limit of free particles without short-
distance scatterers although scattering effects could be straightforwardly
included on the right-hand-side of Eq. (12) in the form of a “collision
integral” $\mathcal{I}_{\rm coll}$ Mangeolle  _et al._ (2022a, b). Here
$\upsilon=\hslash/\hbar$ is a symbolic parameter that counts the order of the
semi-classical expansion, shorn of any dimension: taken literally $\upsilon=1$
since $\hslash=\hbar$. The total number of powers of $\hslash$ and $\upsilon$
in an expression gives the semi-classical order of the corresponding term in
the kinetic equation, and is equal to the number of spatial gradients.

In Eq. (12) we have defined the object

$\partial_{t}{\sf
q}_{\alpha}\equiv\epsilon_{\alpha\beta}\left(\partial_{\beta}\underline{\sf
K}_{d}+\hslash\epsilon_{\gamma\sigma}{\sf\Omega}_{\beta\gamma}\partial_{\sigma}\underline{\sf
K}_{d}\right)+O(\hslash^{2}).$ (13)

Mathematically $\partial_{t}{\sf q}_{\alpha}$ is just a function of phase
space given by the right hand side of Eq. (13), i.e. there is no need to solve
a dynamical equation for ${\sf q}_{\alpha}(t)$ to solve the Boltzmann
equation. However, we can identify physically Eq. (13) as the single particle
equation of motion in phase space of a particle with phase space coordinates
${\sf q}$. We also introduced

$\displaystyle{\sf\Omega}_{\alpha\beta}$
$\displaystyle=\partial_{\alpha}\mathsf{A}_{\beta}-\partial_{\beta}\mathsf{A}_{\alpha}$
(14)

which is the Berry curvature in phase space. Note that this quantity is band-
diagonal by definition and ${\sf\Omega}_{\alpha\beta}={\rm
Im}\left[\mathsf{\Lambda}_{\beta}\,\overset{\star}{,}\,\mathsf{\Lambda}_{\alpha}\right]^{\rm(d)}$,
where we defined the star/Moyal bracket as $[{\sf A}\,\overset{\star}{,}\,{\sf
B}]\equiv{\sf A}\star{\sf B}-{\sf B}\star{\sf A}$. This result, which is
derived systematically using only the assumptions in Sec. II.1, reproduces the
expected form from the heuristic wave packet theory Sundaram and Niu (1999).

By re-expressing physical quantities in the gauge invariant diagonal
variables, we analogously obtain explicitly gauge-invariant expressions for
them. The simplest is the energy density in phase space,

$\displaystyle\mathcal{H}(q)$ $\displaystyle=\frac{1}{2}{\rm
Tr}\Big{[}\,\mathfrak{J}\,\underline{\sf F}_{d}\,\Big{]}+O(\hbar^{2}),$ (15)

where

$\displaystyle\mathfrak{J}=1+\hslash\,\sOmega_{X_{\mu}p_{\mu}}+O(\hslash^{2})$
(16)

is a Jacobian in phase space, see Sec. D. It has the physical implication of
inducing on the density of states a dependence on Berry curvature Xiao  _et
al._ (2005).

Similarly, we obtain the current density in phase space

$\mathcal{J}_{\alpha}(q)=\frac{\upsilon}{2}{\rm
Tr}\left[\,\mathfrak{J}\,\partial_{t}{\sf q}_{\alpha}\,\underline{\sf
F}_{d}+\hslash\epsilon_{\alpha\beta}\epsilon_{\gamma\lambda}\partial_{\gamma}\left(\mathsf{M}_{\beta\lambda}\underline{\sf
F}_{d}\right)\right]+O(\upsilon\hslash^{2}).$ (17)

One should keep in mind that $\alpha$ here is a general phase space index, so
that this is more general than the usual current density. Specifically
$\mathcal{J}_{X_{\mu}}(X,p)$ is the contribution of states at momentum $p$ to
the energy current density at position $X$, while $\mathcal{J}_{p_{\mu}}(X,p)$
gives the analogous contribution to the “force density”.

In Eq. (17), we introduced

$\mathsf{M}_{\alpha\beta}=\frac{1}{2}{\rm
Im}\left\\{\mathsf{\Lambda}_{\beta},\left[\mathsf{\Lambda}_{\alpha},{\sf
K}_{d}\right]\right\\}^{({\rm d})},$ (18)

which has the interpretation of an energy magnetization Cooper  _et al._
(1997); Qin  _et al._ (2011).

###  II.4 Practical remarks

Note that in Eq. (17) we have separated out two contributions to the current
density, in which the second is manifestly a total derivative (phase space
curl). The curl term gives zero when integrated over phase space, and most
importantly gives zero net (momentum integrated) flux through a surface in
real space for which the magnetization $\mathsf{M}_{\beta\lambda}$ vanishes on
the boundary (i.e. a surface which cuts the entire sample volume) — see also
Appendix E. Thus the energy magnetization does not contribute to the total
‘transport current’, and for transport purposes only the first term
(containing $\partial_{t}{\sf q}_{\alpha}$) need be considered.

To use the above results in a specific problem, one needs to specify the
diagonalized band energy ${\sf K}_{d}$ and the connection
$\mathsf{\Lambda}_{\alpha}$. These are obtained by solving the diagonalization
problem perturbatively in $\hslash$. For the energy, there is a simple general
result to first order,

$\displaystyle\underline{\sf K}_{d}(q)$ $\displaystyle={\sf
K}_{d}[\underline{\smash{q}}]$ (19) $\displaystyle={\sf
K}_{0,d}-\hslash\frac{\epsilon_{\alpha\beta}}{2}\mathsf{M}_{\alpha\beta}+\hslash\,{\rm
Re}\left[\left({\sf S}_{0}^{-1}{\sf K}_{1}{\sf S}_{0}\right)^{({\rm
d})}\right]+O(\hslash^{2}),$

where ${\sf K}_{0,d}$ is a diagonal matrix containing the conventional
eigenvalues of the matrix $\mathsf{K}_{0}=\hbar\mathsf{\Gamma}\mathsf{H}$, and
the connection $\mathsf{\Lambda}_{\alpha}$ may be approximated at the same
accuracy by the connection obtained from the this eigenvalue problem at zeroth
order, i.e.
$\mathsf{\Lambda}_{\alpha}\approx\mathsf{S}_{0}^{-1}\partial_{\alpha}\mathsf{S}_{0}$,
with $\mathsf{S}_{0}^{-1}\mathsf{K}_{0}\mathsf{S}_{0}={\sf K}_{0,d}$. The
final term in the square brackets in Eq. (19) is a correction arising from the
semi-classical expansion of $\mathsf{K}$ itself, with ${\sf K}={\sf
K}_{0}+\hslash{\sf K}_{1}+O(\hslash^{2})$ (see App. B.2). We will show that
the spectrum of $\mathsf{K}$ is symmetric with positive and negative energy
eigenstates related by complex conjugation and momentum reversal (see Appendix
B.3).

It can be useful both for intuition and as a check to consider the dimensions
of the various objects encountered in this treatment. This is somewhat problem
dependent, owing to different possible choices of the fundamental $\Phi_{a}$
fields, which might be position and momentum densities for phonons, or
transverse magnetization densities in a ferromagnet, for example. Thus
different components might even have different dimensions. However, by
appropriate rescaling, it is generally possible to bring them to a form in
which all $\Phi_{a}$ fields have dimension of the square root of density,
$[\Phi_{a}]=L^{-d/2}$. Then one obtains

$\displaystyle[X_{\mu}]=L,\quad[k_{\mu}]=1/L,\quad[p_{\mu}]=Et/L,\quad[\mathfrak{J}]=1,$
$\displaystyle[\mathcal{H}]=[{\sf K}]=[{\sf K}_{d}]=[{\sf F}_{d}]=E,$
$\displaystyle[\mathcal{J}_{X}]=EL/t,\quad[\mathcal{J}_{p}]=E^{2}/L,$
$\displaystyle[{\sf\Lambda}_{\alpha}]=[{\sf A}_{\alpha}]=1/[\alpha],$
$\displaystyle[{\sf\Omega}_{\alpha\beta}]=1/([\alpha][\beta]),\quad[{\sf
M}_{\alpha\beta}]=E/([\alpha][\beta]),$ $\displaystyle[{\sf H}][{\sf
F}]=E,\quad[{\sf H}][{\sf\Gamma}]=1/t,\quad[{\sf S}]^{2}=1/E,$ (20)

where $E$, $L$ and $t$ stand for energy, length and time, respectively. Note
that in Eq. (II.4), all quantities from $\mathfrak{J}$ onward are defined in
phase space. This brings a factor of $L^{d}$ to many quantities owing to the
Wigner transform, which contains a spatial integral, for example
$[\hat{\mathsf{H}}]=E/L^{d}$ is an energy density while $[\mathsf{H}]=E$ is an
energy (if $[\Phi_{a}]=L^{-d/2}$). Note that all the quantities in Eq. (II.4)
are of order $O(\hslash^{0})$. This is irrespective of our choice of
definitions for $p_{\mu}$ and ${\sf K}$ for example, which we chose to be
momentum $p_{\mu}=\hbar k_{\mu}=O(\hslash^{0})$ and energy ${\sf
K}=\hbar\,{\sf\Gamma}\star{\sf H}=O(\hslash^{0})$ as opposed to wavevector and
frequency, respectively.

##  III Detailed derivation

Here we provide the derivation of the results derived in Section II.

###  III.1 Dynamics and continuity equation

Together, Eqs. (1-2) determine the dynamics of the fields according to the
Heisenberg equation of motion,

$\partial_{t}\Phi_{a}=-\frac{i}{\hbar}[\Phi_{a},{\rm
H}]=\frac{1}{\hbar}\left(\hat{\mathsf{K}}\otimes\Phi\right)_{a},$ (21)

where convolution $\otimes$ is defined as

$(\hat{\sf O}_{1}\otimes\hat{\sf
O}_{2})_{ab}(r_{1},r_{2})\equiv\int_{r}\sum_{c}\hat{\sf
O}_{1,ac}(r_{1},r)\hat{\sf O}_{2,cb}(r,r_{2}),$ (22)

and the dynamical matrix is

$\hat{\sf K}\equiv\hbar\,\hat{\sf\Gamma}\otimes\hat{\sf H}.$ (23)

In turn, this implies that

$\partial_{t}\hat{\sf F}=\frac{-i}{\hbar}(\hat{\sf K}\otimes\hat{\sf
F}-\hat{\sf F}\otimes\hat{\sf K}^{\dagger}),$ (24)

where $(\hat{\sf K}^{\dagger})_{ab}(r_{1},r_{2})\equiv\left(\hat{\sf
K}_{ba}(r_{2},r_{1})\right)^{*}$. In the same notation the total energy is

$H\equiv\langle{\rm H}\rangle=\frac{1}{2}\int_{r}{\rm Tr}[(\hat{\sf
H}\otimes\hat{\sf F})(r,r)],$ (25)

where ${\rm Tr}$ is the usual matrix trace, i.e. ${\rm Tr}[\hat{\sf
O}]=\sum_{a}\hat{\sf O}_{aa}$.

Now we pass to the Wigner transform, and use the fact that the Wigner
transform of a convolution is the star (or Moyal) product of the Wigner
transforms, i.e. $[\hat{f}_{1}\otimes\hat{f}_{2}]^{W}(X,p)=f_{1}(X,p)\star
f_{2}(X,p)$ (see Appendix A). This leads directly to Eq. (5) and Eq. (6).

In this way, the total energy can also be written as

$\langle{\rm H}\rangle\equiv H=\frac{1}{2}\int_{X,p}{\rm Tr}[{\sf H}\star{\sf
F}],$ (26)

where ${\rm Tr}[{\sf O}]=\sum_{a}{\sf O}_{aa}$ for any matrix ${\sf O}$ is the
regular matrix trace.

Conservation of energy is ensured by Heisenberg evolution, Eq. (21), which is
unitary, and when the Hamiltonian is local, this gives rise to a continuity
equation. To obtain it we first define from Eq. (26) a phase space energy
density

$\mathcal{H}(X,p)=\frac{1}{4}{\rm Tr}\left({\sf H}\star{\sf F}+{\sf
F}\star{\sf H}\right)\equiv\frac{1}{4}{\rm Tr}\left\\{{\sf
H}\,\overset{\star}{,}\,{\sf F}\right\\}.$ (27)

Here we symmetrized by hand the argument of the trace, so that the energy
density is real, $\mathcal{H}=\mathcal{H}^{*}$ (note that the choice of local
energy density is not unique Kapustin and Spodyneiko (2020); Hardy (1963)).

Taking the time derivative of the energy density, we find

$\partial_{t}\mathcal{H}=\frac{1}{2\hbar}{\rm Im}{\rm Tr}\left([{\sf
K}\,\overset{\star}{,}\,{\sf F}\star{\sf H}]\right).$ (28)

Very generally, the phase space integral of the trace of a Moyal bracket
vanishes (see Appendix A). Hence, the right-hand-side of Eq. (28) can can be
written as a phase-space divergence, i.e. there exists
$\boldsymbol{\mathcal{J}}(X,p)$ such that, in the collisionless limit,

$\partial_{t}\mathcal{H}+\partial_{\alpha}\mathcal{J}_{\alpha}=0.$ (29)

Up to this point, all the above formulas are exact.

Expanding the “center” of the Moyal bracket (i.e. “$\overset{\star}{,}$”) in
$\hslash$ in Eq. (28), we find

$\mathcal{J}_{\alpha}(X,p)=\frac{\upsilon}{2}\epsilon_{\alpha\beta}{\rm
Re}\,{\rm Tr}\left(\partial_{\beta}{\sf K}({\sf F}\star{\sf
H})\right)+O(\upsilon\hslash^{2}),$ (30)

where ${\rm Re}$ is the real part. We note that the local phase-space current
density $\boldsymbol{\mathcal{J}}$ depends on the precise choice of local
energy density and enjoys an additional freedom in the form of an additive
phase-space curl,
$\mathcal{J}_{\alpha}\rightarrow\mathcal{J}_{\alpha}+\partial_{\beta}\mathcal{M}_{\alpha\beta}$
for any differentiable phase-space function $\mathcal{M}_{\alpha\beta}$ odd
under $\alpha\leftrightarrow\beta$, provided that the curl is gauge invariant.

To summarize, the energy density is defined by Eq. (27), and the current
density by Eq. (30). Now by expressing these quantities in the diagonal band
basis, we will eventually arrive at Eq. (15) and Eq. (17).

###  III.2 Diagonalization

The effects of curvature and geometry arise from the projection into the
manifold of a band Vanderbilt (2018). Formally the projection is carried out
in the phase space representation by diagonalizing the dynamical matrix using
a similarity transformation $\mathsf{S}$, introduced in Sec. II.3

$\displaystyle\mathsf{S}^{-1}\star\mathsf{K}\star\mathsf{S}$
$\displaystyle=\mathsf{K}_{d},$ (31)

where $\mathsf{K}_{d}$ is a diagonal matrix (whose entries are real, as shown
in Appendix B.1).

Because $\mathsf{K}$ is not hermitian, $\mathsf{S}$ cannot be chosen star-
unitary, but instead can be chosen to satisfy (see Appendix B.1)

$\displaystyle\mathsf{S}^{\dagger}\star\mathsf{H}\star\mathsf{S}$
$\displaystyle=1,$ (32)

where $1$ is the identity matrix. Eq. (32) is the analog of the star-unitarity
condition for a non-unitary matrix. Note that we define the inverse here in
the sense of the star product, i.e.
$\mathsf{S}^{-1}\star\mathsf{S}=\mathsf{S}\star\mathsf{S}^{-1}=1$, as
mentioned in Sec. II.3.

Eqs. (31,32) are somewhat formal, due to the non-trivial nature of the star
product, but can be solved order by order in $\hslash$, starting from a zeroth
order solution which is a standard matrix diagonalization problem, and which
is guaranteed to exist by the properties of $\mathsf{K}$.

The similarity transformation $\mathsf{S}$ is therefore determined (up to a
gauge freedom we will return to below) by $\mathsf{K}$, and once it is found,
we can, as mentioned in Sec. II.3, define the transformed density matrix
$\mathsf{F}_{d}$ as well via

$\displaystyle\mathsf{F}=\mathsf{S}\star\mathsf{F}_{d}\star\mathsf{S}^{\dagger}.$
(33)

Note that because $\mathsf{S}$ is not unitary but instead satisfies Eq. (32),
$\mathsf{F}_{d}$ differs dimensionally from $\mathsf{F}$ by a factor of
energy. In this “diagonal frame,” and in the collisionless limit, one has

$\partial_{t}{\sf F}_{d}+\frac{i}{\hbar}[{\sf K}_{d}\,\overset{\star}{,}\,{\sf
F}_{d}]=0.$ (34)

This equation, which is exact, allows a diagonal and time-independent solution
for $\mathsf{F}_{d}$. This is the standard form of the distribution function
at and near equilibrium in the following sense. There is nothing prohibiting
off-diagonal terms in $\mathsf{F}_{d}$, but such solutions necessarily
oscillate on the scale of the differences between elements in
$\mathsf{K}_{d}/\hbar$, i.e. the mode frequencies. This means that slowly
time-varying solutions are predominantly diagonal. Moreover, when weak
scattering is included, such that the associated relaxation time $\tau$ is
longer than these oscillation periods, the effects of such oscillations
largely average out. In any case, in the present formulation without
scattering, there is no mixing of the diagonal and off-diagonal parts of
$\mathsf{F}_{d}$, so we can focus on the former consistently.

###  III.3 Gauge transformations

####  III.3.1 General

As mentioned above, Eqs. (31,32) leave some freedom in the choice of
$\mathsf{S}$. In particular, if $\mathsf{S}$ satisfies these equations in the
sense that it produces a diagonal $\mathsf{K}_{d}$, then so too does

${\sf S}^{\prime}={\sf S}\star{\sf\Theta},$ (35)

produce a diagonal ${\sf K}_{d}^{\prime}$, where ${\sf\Theta}$ is a diagonal
(star-) unitary matrix,

$\displaystyle\mathsf{\Theta}^{\dagger}\star\mathsf{\Theta}=\mathsf{\Theta}\star\mathsf{\Theta}^{\dagger}=1.$
(36)

Note that the diagonalized matrix $\mathsf{K}^{\prime}_{d}$ is not generally
equal to $\mathsf{K}_{d}$. Rather, if we view the change from $\mathsf{S}$ to
$\mathsf{S}^{\prime}$ as a transformation,

$\displaystyle\mathsf{S}\rightarrow\mathsf{S}^{\prime}={\sf
S}\star{\sf\Theta},$ (37)

then

$\displaystyle\mathsf{K}_{d}\rightarrow\mathsf{K}^{\prime}_{d}=\mathsf{\Theta}^{\dagger}\star\mathsf{K}_{d}\star\mathsf{\Theta},$
(38)

which is manifestly also diagonal and hermitian. The same transformation law
holds for the density matrix,

$\displaystyle\mathsf{F}_{d}\rightarrow\mathsf{\Theta}^{\dagger}\star\mathsf{F}_{d}\star\mathsf{\Theta}.$
(39)

Physical quantities defined in terms of
$\mathsf{H},\mathsf{\Gamma},\mathsf{F}$ must, by construction, be invariant
under these gauge transformations.

We will now show, perturbatively in $\hslash$, that such gauge transformations
can be interpreted as coordinate changes in phase space.

####  III.3.2 First order in $\hbar$

We can express explicitly

$\displaystyle\mathsf{\Theta}$
$\displaystyle=\textrm{diag}\,\left(e^{i\theta_{1}},\cdots,e^{i\theta_{2N}}\right),$
(40)

where the condition of unitarity, Eq. (36) perturbatively reduces to

$\displaystyle{\rm Im}\,\theta_{a}=O(\hslash^{2}).$ (41)

Hence to $O(\hslash)$, we can use Eq. (40) with real-valued
$\theta_{a}(q)=\theta_{a}(X,p)$, which in general are arbitary functions on
phase space.

Using the explicit form of $\mathsf{\Theta}$ in Eq. (40), Eqs. (38,39) reduce
to

$\displaystyle{\sf O}_{d}(q)$ $\displaystyle\rightarrow{\sf
O}^{\prime}_{d}(q)={\sf
O}_{d}(q)-\hslash\epsilon_{\alpha\beta}\partial_{\beta}\stheta\partial_{\alpha}{\sf
O}_{d}(q)+O(\hslash^{2}),$ $\displaystyle={\sf
O}_{d}[q^{\prime}]+O(\hslash^{2}),$ (42)

for $\mathsf{O}_{d}=\mathsf{K}_{d},\mathsf{F}_{d}$, where
$\stheta=\textrm{diag}\,(\theta_{1},\cdots,\theta_{2N})$, and

$\displaystyle q^{\prime}_{\alpha}$
$\displaystyle=q_{\alpha}-\hslash\epsilon_{\alpha\beta}\partial_{\beta}\stheta.$
(43)

In Eq. (III.3.2) and Eq. (43), we have noticed that the change of ${\sf
O}_{d}$ at order $\hslash$ under a gauge transformation is proportional to a
derivative of ${\sf O}_{d}$, and hence can be absorbed into an equivalent
shift of its argument.

Note that $q^{\prime}_{\alpha}$ in Eq. (43) contains $\stheta$ which is a
diagonal matrix, i.e. a band-dependent quantity, and so Eq. (III.3.2) defines
a matrix function of a matrix. This is not ambiguous because both the matrix
and the matrix argument are diagonal, and so can simply be evaluated for each
band. More explicitly, the diagonal matrix function $\mathsf{O}$ of a diagonal
matrix argument $\mathsf{C}$ is defined as

${\sf O}[{\sf C}]\equiv{\rm diag}\Big{(}\big{\\{}{\sf O}_{aa}({\sf
C}_{aa})\big{\\}}_{a}\Big{)}.$ (44)

###  III.4 Gauge connection and gauge invariant quantities

From Eq. (III.3.2) it is evident that
$\mathsf{O}_{d}=\mathsf{F}_{d},\mathsf{K}_{d}$ are not gauge invariant. This
ultimately reflects the presence of Berry phases. In particular, we see that
the effect of a gauge transformation is to shift the “canonical” variable
$q_{\alpha}$ within the argument of ${\sf O}_{d}$ by a quantity proportional
to $\partial_{\beta}\stheta$. This motivates a procedure to restore gauge
invariance by modifying the argument of $\mathsf{O}_{d}$ to compensate for
this shift.

This compensation is obtained from the Berry gauge field introduced in Eq.
(II.3). Under the gauge transformation Eq. (37), to zeroth order

${\sf A}_{\alpha}\rightarrow{\sf
A}_{\alpha}+\partial_{\alpha}\stheta+O(\hslash).$ (45)

We will see that $\mathsf{A}_{\alpha}$ as well as $\mathsf{\Lambda}_{\alpha}$
emerge naturally in calculations of physical quantities. Comparing Eq. (45)
and Eq. (43) shows how to compensate the transformation of $\mathsf{O}_{d}$ by
changing its argument.

By combining Eqs. (9,10, III.3.2,45), we can see that $\underline{\sf
O}_{d}(q)$ is gauge invariant to first order,

$\displaystyle\underline{\sf O}^{\prime}_{d}(q)$ $\displaystyle=\underline{\sf
O}_{d}(q)+O(\hslash^{2}).$ (46)

It is expected that physical quantities can be expressed as manifestly gauge
invariant expressions when written in terms of such functions, i.e.
$\underline{\mathsf{K}}_{d},\underline{\mathsf{F}}_{d}$.

###  III.5 Gauge invariant kinetic equation

To derive the gauge-invariant band-diagonal kinetic equation, we
semiclassically expand the Moyal bracket in Eq. (34):

$\displaystyle\partial_{t}\mathsf{F}_{d}$
$\displaystyle=\upsilon\epsilon_{\alpha\beta}\partial_{\alpha}{\sf
K}_{d}\partial_{\beta}\mathsf{F}_{d}+O(\upsilon\hslash^{2}).$ (47)

Next we take the time derivative of Eq. (9) with
$\mathsf{O}_{d}=\mathsf{F}_{d}$, which gives

$\displaystyle\partial_{t}\underline{\mathsf{F}}_{d}(q)$
$\displaystyle=\partial_{t}\mathsf{F}_{d}[\underline{\smash{q}}],$ (48)

because $\mathsf{A}_{\beta}$ is time independent by assumption since we did
not include any time-dependence in ${\rm H}$. Hence from Eq. (47)

$\displaystyle\partial_{t}\underline{\mathsf{F}}_{d}(q)$
$\displaystyle=\upsilon\epsilon_{\alpha\beta}\left.\partial_{\alpha}{\sf
K}_{d}(q)\partial_{\beta}\mathsf{F}_{d}(q)\right|_{q\rightarrow\underline{\smash{q}}}+O(\upsilon\hslash^{2}).$
(49)

Now we need to convert the functions on the right hand side above to
underlined quantities. To do so, we take the phase space derivative of Eq.
(9),

$\displaystyle\partial_{\alpha}\underline{\mathsf{O}}_{d}(q)$
$\displaystyle=\partial_{\alpha}\mathsf{O}_{d}[\underline{\smash{q}}]+\hslash\epsilon_{\gamma\lambda}\partial_{\gamma}\mathsf{O}_{d}\,\partial_{\alpha}\mathsf{A}_{\lambda}+O(\hslash^{2}),$
$\displaystyle\Rightarrow\partial_{\alpha}\mathsf{O}_{d}[\underline{\smash{q}}]$
$\displaystyle=\partial_{\alpha}\underline{\mathsf{O}}_{d}(q)-\hslash\epsilon_{\gamma\lambda}\partial_{\gamma}\underline{\mathsf{O}}_{d}\,\partial_{\alpha}\mathsf{A}_{\lambda}+O(\hslash^{2}).$
(50)

Using Eq. (III.5) for both ${\sf K}_{d}$ and $\mathsf{F}_{d}$ in the right
hand side of Eq. (49), we arrange terms to obtain

$\displaystyle\partial_{t}\underline{\mathsf{F}}_{d}$
$\displaystyle=\upsilon\bigg{[}\epsilon_{\alpha\beta}\partial_{\alpha}\underline{\sf
K}_{d}\partial_{\beta}\mathsf{F}_{d}-\hslash\epsilon_{\alpha\beta}\epsilon_{\gamma\lambda}\partial_{\gamma}\underline{\sf
K}_{d}\partial_{\beta}\underline{\mathsf{F}}_{d}\partial_{\alpha}\mathsf{A}_{\lambda}$
$\displaystyle\quad-\hslash\epsilon_{\alpha\beta}\epsilon_{\gamma\lambda}\partial_{\alpha}\underline{\sf
K}_{d}\partial_{\gamma}\underline{\mathsf{F}}_{d}\partial_{\beta}\mathsf{A}_{\lambda}+O(\hslash^{2})\bigg{]}.$
(51)

Now relabeling dummy indices in the last term above according to
$\alpha\rightarrow\gamma$, $\beta\rightarrow\lambda$,
$\gamma\rightarrow\beta$, $\lambda\rightarrow\alpha$, one obtains Eq. (12)
with $\partial_{t}{\sf q}_{\alpha}$ as in Eq. (13).

###  III.6 Gauge invariant energy and current density

Now we proceed to derive diagonal and gauge invariant expressions for the
energy density, Eq. (27) and energy current density, Eq. (30).

####  III.6.1 Energy density

Let us begin with the energy density. Using properties of the trace and others
recapped in Appendix A, we write Eq. (27) as

$\displaystyle\mathcal{H}$
$\displaystyle=\frac{1}{2}\textrm{Re}\,\textrm{Tr}\,\left(\mathsf{F}\star\mathsf{H}\right)$
$\displaystyle=\frac{1}{2}\textrm{Re}\,\textrm{Tr}\,\left(\mathsf{S}\star\mathsf{F}_{d}\star\mathsf{S}^{-1}\right),$
(52)

using Eqs. (31,32,33).

To simplify this, we use the identity, valid for any operator ${\sf O}$,

$\displaystyle\textrm{Tr}\,\left(\mathsf{S}\star\mathsf{O}\star\mathsf{S}^{-1}\right)$
$\displaystyle=\textrm{Tr}\,\mathsf{O}+i\hslash\epsilon_{\alpha\beta}\partial_{\beta}\textrm{Tr}\,\left(\mathsf{\Lambda}_{\alpha}\mathsf{O}\right)+O(\hslash^{2}),$
(53)

which is derived straightforwardly by the expansion of the star product.
Taking $\mathsf{O}=\mathsf{F}_{d}$ and applying this to Eq. (III.6.1), we then
obtain

$\displaystyle\mathcal{H}$ $\displaystyle=\frac{1}{2}\left({\rm Tr}{\sf
F}_{d}-\hslash\epsilon_{\alpha\beta}\partial_{\beta}{\rm Tr}\left[{\sf
A}_{\alpha}{\sf F}_{d}\right]\right)+O(\hslash^{2})$ (54)
$\displaystyle=\frac{1}{2}{\rm
Tr}\Big{[}\left(1+\frac{\hslash}{2}\epsilon_{\gamma\lambda}{\sf\Omega}_{\gamma\lambda}\right)\left({\sf
F}_{d}+\hslash\epsilon_{\alpha\beta}{\sf A}_{\beta}\partial_{\alpha}{\sf
F}_{d}\right)\Big{]}+O(\hslash^{2}),$

where we defined (the band-diagonal quantity) ${\sf\Omega}_{\alpha\beta}$ in
Eq. (14). This expression is diagonalized and gauge-invariant but the latter
gauge invariance is not explicit. To make it so, we use the definition in Eq.
(9) and Taylor expand in $\hslash$ at fixed $q_{\alpha}$. This gives

$\displaystyle\underline{\mathsf{O}}_{d}(q)$
$\displaystyle=\mathsf{O}_{d}(q)+\hslash\epsilon_{\alpha\beta}\mathsf{A}_{\beta}\partial_{\alpha}\mathsf{O}_{d}(q)+O(\hslash^{2}).$
(55)

Applying Eq. (55) with $\mathsf{O}_{d}=\mathsf{F}_{d}$ immediately gives Eqs.
(15,16) in Sec. II. The fact that $\mathfrak{J}$, defined in Eq. (16),
identifies with the jacobian of the $q\mapsto\underline{\smash{q}}$
transformation is shown in App. D.

####  III.6.2 Energy current

Now we turn to the energy current, written in Eq. (30). It contains two
factors inside the trace, $\partial_{\beta}{\sf K}$, and
$\mathsf{F}\star\mathsf{H}$, which we want to write in diagonal form. To
express the former, we note the useful identity

$\displaystyle\partial_{\alpha}\left(\mathsf{S}\star\mathsf{O}\star\mathsf{S}^{-1}\right)$
$\displaystyle=\mathsf{S}\star\left(\partial_{\alpha}\mathsf{O}+\left[{\sf\Lambda}_{\alpha}\overset{\star}{,}\,\mathsf{O}\right]\right)\star\mathsf{S}^{-1}.$
(56)

Applying Eq. (56) with $\mathsf{O}={\sf K}_{d}=\mathsf{S}^{-1}\star{\sf
K}\star\mathsf{S}$ gives

$\partial_{\alpha}{\sf K}={\sf S}\star{\sf Q}_{\alpha}\star{\sf S}^{-1},$ (57)

where we defined for convenience

${\sf Q}_{\alpha}\equiv\partial_{\alpha}{\sf K}_{d}+\left[{\sf
S}^{-1}\star\partial_{\alpha}{\sf S}\,\overset{\star}{,}\,{\sf K}_{d}\right].$
(58)

The second multiplicative factor in the trace in Eq. (30) is, using Eqs.
(32,33),

$\displaystyle\mathsf{F}\star\mathsf{H}$
$\displaystyle=\mathsf{S}\star\mathsf{F}_{d}\star\mathsf{S}^{-1}.$ (59)

Using Eq. (57) and Eq. (59), one can convert Eq. (30) to diagonal form using
manipulations similar to those for the energy density, though with more
involved algebra. The details are left to Appendix C. One obtains in this way

$\mathcal{J}_{\alpha}=\mathcal{J}_{\alpha}^{(1)}+\mathcal{J}_{\alpha}^{(2)},$
(60)

in which each term is separately gauge invariant:

$\displaystyle\mathcal{J}_{\alpha}^{(1)}$
$\displaystyle=\frac{\upsilon}{2}\epsilon_{\alpha\beta}{\rm
Tr}\Big{[}\big{(}1+\frac{\hslash}{2}\epsilon_{\gamma\lambda}{\sf\Omega}_{\gamma\lambda}\big{)}$
$\displaystyle\quad\times\big{(}\partial_{\beta}({\sf
K}_{d}+\hslash\epsilon_{\rho\sigma}{\sf A}_{\sigma}\partial_{\rho}{\sf
K}_{d})+\hslash\epsilon_{\rho\sigma}{\sf\Omega}_{\beta\rho}\partial_{\sigma}{\sf
K}_{d}\big{)}$ $\displaystyle\quad\times\big{(}{\sf
F}_{d}+\hslash\epsilon_{\xi\eta}{\sf A}_{\eta}\partial_{\xi}{\sf
F}_{d}\big{)}\Big{]}+O(\upsilon\hslash^{2}),$ (61a)
$\displaystyle\mathcal{J}_{\alpha}^{(2)}$
$\displaystyle=\frac{\upsilon\hslash}{2}\epsilon_{\alpha\beta}\epsilon_{\gamma\lambda}\partial_{\gamma}{\rm
Tr}\big{[}\mathsf{M}_{\beta\lambda}{\sf F}_{d}\big{]}+O(\upsilon\hslash^{2}),$
(61b)

where $\mathsf{M}_{\alpha\beta}$ was defined in Eq. (18). Now we can again use
Eq. (55) to identify gauge-invariant functions $\underline{\sf K}_{d}$ and
$\underline{\mathsf{F}}_{d}$, and we directly thereby obtain Eq. (17) in Sec.
II.

Let us discuss the physical meaning and relevance of
$\mathcal{J}_{\alpha}^{(2)}$. This term satisfies the continuity equation _by
definition_ , as it is a phase space “curl”. Such terms appear to be ambiguous
from the introduction of the heat current, which occurred in passing from Eq.
(28) to Eq. (30). That ambiguity arises because the quantity which enters the
continuity equation is only the divergence of the phase space current, so that
a shift
$\mathcal{J}_{\alpha}\rightarrow\mathcal{J}_{\alpha}+\partial_{\beta}\mathcal{M}_{\alpha\beta}$
with antisymmetric $\mathcal{M}$, leaves the continuity equation unchanged.

The above derivation certainly provides an unambiguous result for
$\mathcal{J}_{\alpha}^{(2)}$ given the definition in Eq. (30). However, we
would like to ask whether the result in Eq. (61b) can be considered to resolve
the intrinsic ambiguity. To resolve it, one must clearly impose some other
conditions. Related conditions are discussed for the energy magnetization
itself in Ref. Qin  _et al._ (2011). For the energy magnetization current, we
are not sure of the full set of conditions that should be imposed, but we do
note that the above definition satisfies many reasonable ones:

  1. 1.

The quantity ${\sf M}_{\alpha\beta}$ itself is a gauge-invariant band-diagonal
quantity which can be identified as an energy magnetization. Then,
$\mathcal{J}^{(2)}_{\alpha}$ is naturally the gauge-invariant divergenceless
current which can be built from it, and it indeed directly appears when
extracting the divergence from the time-derivative of the local energy density
$\partial_{t}\mathcal{H}$.

  2. 2.

The magnetization contribution, $\mathcal{J}_{\alpha}^{(2)}$ vanishes in
equilibrium when the system is locally uniform. We may define equilibrium to
be such that the energy density ${\sf F}_{d}$ depends only on the band energy
${\sf K}_{d}$, i.e. ${\sf F}_{d}(q)=f[{\sf K}_{d}(q)]$, then
$\partial_{\gamma}{\sf F}_{d}=\partial_{\gamma}{\sf K}_{d}f^{\prime}[{\sf
K}_{d}(q)]$. Then $\mathcal{J}^{(2)}_{\alpha}$ vanishes if $\partial_{X}{\sf
K}_{d}=0$ and $\partial_{X}{\sf S}=0$. It does not, however, a priori vanish
otherwise, even in equilibrium. A local nonvanishing current in equilibrium is
in fact how the magnetization current has been defined in the literature
Cooper  _et al._ (1997).

  3. 3.

One might expect that in equilibrium, even the local thermal current
$J_{X_{\mu}}=\int_{p}(\mathcal{J}_{X_{\mu}}^{(1)}+\mathcal{J}_{X_{\mu}}^{(2)})$,
vanishes at zero temperature. This would seem natural since the entropy
vanishes at zero temperature. We are able to show that this is the case in the
special but important case of a $(p,X)$-separable problem such as we describe
below and in Appendix F.

  4. 4.

As noted in Refs. Qin  _et al._ (2011, 2012), the momentum integral defining
$J_{X_{\mu}}$ is rid of a UV divergence when evaluated in a continuum theory,
in the $(p,X)$-separable case. This is not the case for the integrands
$\mathcal{J}_{X_{\mu}}^{(1)}$ and $\mathcal{J}_{X_{\mu}}^{(2)}$ considered
separately. Thus Eq. (61b) is “physical” in the sense that it provides a
compensating term to an otherwise unphysical divergence in Eq. (61a).

##  IV Application: thermal Hall current in a finite system of chiral bosons

We now show how our QKE formalism can be applied to any general system with a
boundary but homogeneous away from the boundary, and further specialize to a
system inhomogeneous linear elasticity theory with a phonon Hall viscosity
term.

###  IV.1 Uniform finite system with a boundary: formal result for thermal
conductivity

One natural application of our formalism is to a finite system, in which the
presence of a boundary is included via variations of the Hamiltonian in space.
Specifically, we assume that the (2d) system is constrained within the
interval $I_{x}=\left\\{x\in\left[-\frac{L}{2},+\frac{L}{2}\right]\right\\}$,
and we consider a Hamiltonian of the form ${\sf H}(X,p)=g(X){\sf H}^{\rm
h}(p)$, with $g$ a c-number function equal to 1 in the $I_{x}$ interval and
decaying over distances $O(L^{0})$ to 0 outside, mimicking the
presence/absence of physical degrees of freedom in these regions. This is just
a specific instance of the more general problem of a $(X,p)$-separable theory,
which we consider in detail in Appendix F.

There, we show that solving the separable problem amounts to solving the
homogeneous one (i.e. $g(X)=1$ everywhere), which we label with ${\rm h}$
superscripts:

$\displaystyle(\mathsf{S}^{\rm h})^{-1}\mathsf{K}^{\rm h}\mathsf{S}^{\rm h}$
$\displaystyle=\mathsf{K}_{d}^{\rm h},$ (62a) $\displaystyle(\mathsf{S}^{\rm
h})^{\dagger}\mathsf{H}^{\rm h}\mathsf{S}^{\rm h}$ $\displaystyle=1,$ (62b)

where the $\star$-product was replaced by the standard matrix product in this
$X$-independent problem, and $\mathsf{K}^{\rm h}=\hbar\,\mathsf{\Gamma}^{\rm
h}\mathsf{H}^{\rm h}$. This is now a standard linear algebra problem, yielding
in particular the homogeneous energies $\mathsf{K}_{d}^{\rm h}$ and Berry
curvatures ${\sf\Omega}^{\rm h}_{p_{y}p_{x}}$.

We assume translational invariance along $y$ and consider a temperature
gradient $\partial_{x}T$ along the $x$ direction. The (“transport”) thermal
conductivity relates the latter to the net energy flux flowing in response—via
Fourier’s law $\boldsymbol{J}_{\rm tot}=-\boldsymbol{\kappa}^{\rm
tr}\cdot\boldsymbol{\nabla}T$— which is obtained by integrating the energy
current along the $x$ direction and taking the thermodynamic limit:

$\displaystyle\kappa_{xy}^{\rm tr}$
$\displaystyle=\frac{1}{\partial_{x}T}\lim_{L\rightarrow\infty}\frac{1}{L}\int_{-\infty}^{+\infty}\text{d}x\int_{p}\mathcal{J}_{X_{y}}(X,p).$
(63)

Note that we denote here $(X_{x},X_{y})\equiv(x,y)$, and we dropped the
$O(\hslash^{2})$ correction, since $\kappa_{xy}^{\rm tr}$ is by definition a
linear response quantity. The explicit form of both currents from the
decomposition Eq. (60) is derived in App. F.3.

Because $\mathcal{J}_{X_{y}}^{(2)}$ is a total real-space derivative, it does
not contribute to the transport current. Similarly, the zeroth-order term
$g\partial_{p_{y}}{\sf K}_{d}^{\rm h}$ in $\mathcal{J}_{X_{y}}^{(1)}$ drops
out as it is a total momentum-space derivative. The only contribution which
does not vanish after performing the momentum integral is

$\displaystyle\mathcal{J}_{X_{y}}^{(\rm eff)}$
$\displaystyle\equiv-\frac{\hslash}{2}\psi_{x}{\rm Tr}\Big{[}{\sf\Omega}^{\rm
h}_{p_{y}p_{x}}g{\sf K}_{d}^{\rm h}f(g{\sf K}_{d}^{\rm
h},T)\Big{]}+O(\hslash^{2}),$ (64)

where $\psi_{x}\equiv\partial_{x}g/g$ and
$f(\varepsilon,T)=\varepsilon\left(n_{\rm
B}(\varepsilon,T)+\tfrac{1}{2}\right)$, with $n_{\rm
B}(\varepsilon,T)=1/(e^{\varepsilon/(k_{\rm B}T)}-1)$ the Bose function
($k_{\rm B}$ is Boltzmann’s constant), is the equilibrium energy distribution
evaluated at $\varepsilon$ of a system of bosons at temperature $T$. Further
decomposing explicitly $\rm Tr=Tr_{+}+Tr_{-}$ where $\rm Tr_{+}$ (resp. $\rm
Tr_{-}$) sums over positive (resp. negative) energy bands (see Appendix B.3),
and the fact that ${\sf\Omega}^{\rm h}_{p_{y}p_{x}},{\sf K}_{d}^{\rm h}$ are
odd upon spectrum reflection while $f$ is even (note that
$f(-\varepsilon,T)=f(\varepsilon,T)$ and see Appendix B.3), allows one to
replace $\rm Tr\mapsto 2Tr_{+}$ in Eq. (64).

Importantly, the gradient of temperature is not introduced “artificially” in
the form of a gravitational field like in Luttinger’s trick Luttinger (1964)
(in particular it is unrelated to $\psi_{x}$), but simply arises from the fact
that the function $f$ depends implicitly on $x$ since the Bose function must
be computed at the _local equilibrium_ temperature $T(x)$ Cooper  _et al._
(1997).

Eqs. (63,64) are our solution to the problem, which is readily amenable to
numerical evaluation. It is also possible to reconcile it with existing
literature with just one extra step. Performing the spatial integration with a
few simple tricks to turn it into an integral over energies $\varepsilon$ (cf.
Appendix G), we arrive at the result

$\displaystyle\kappa_{xy}^{\rm tr}$ $\displaystyle=-\frac{1}{T}{\rm
Tr}_{+}\int_{p}\mathsf{\Omega}_{p_{x}p_{y}}^{\rm
h}\int\text{d}\varepsilon\,\Theta(\varepsilon-{\sf K}_{d}^{\rm
h})\,\varepsilon^{2}\partial_{\varepsilon}n_{B}\left(\varepsilon,T\right),$
(65)

where $\Theta$ is the step function and $T$ the mean temperature within
$I_{x}$. This formula Eq. (65) agrees with prior results Qin  _et al._ (2012);
Matsumoto  _et al._ (2014) which were obtained using the Kubo formalism.

###  IV.2 Local energy current in a chiral phonon system

We now consider a concrete example. The model is that of elasticity with a
time-reversal breaking term, described by the (euclidian) Lagrangian

$\mathcal{L}_{\rm
ph}=\frac{\rho}{2}\left(\partial_{\tau}\boldsymbol{u}\right)^{2}+\frac{1}{2}c_{ij\mu\nu}\partial_{\mu}u^{i}\partial_{\nu}u^{j}+i\eta_{ij\mu\nu}\partial_{\tau}u^{i}\partial^{2}_{\mu\nu}u^{j},$
(66)

where $\boldsymbol{u}=(u^{i})_{i}$ is the displacement field ($i=x,y,z$),
$\rho$ is the mass density, $c_{ij\mu\nu}$ is the bulk modulus
($j,\mu,\nu=x,y,z$), and $\eta_{ij\mu\nu}$ is the lattice (or phonon) Hall
viscosity. To remove ambiguity we take
$c_{ij\nu\mu}=c_{ij\mu\nu}=c_{ji\mu\nu}$ and
$\eta_{ij\nu\mu}=\eta_{ij\mu\nu}=-\eta_{ji\mu\nu}$.

The Hall viscosity term breaks time reversal symmetry explicitly, and accounts
for the chiral nature of bosonic excitations in the system (chiral phonons).
The physical origin of this term has been discussed elsewhere in various
contexts such as quantum Hall systems Barkeshli  _et al._ (2012), magnetic
insulators Ye  _et al._ , ionic crystals Saito  _et al._ (2019), etc. He we
simply use Eq. (66) as our starting point, as an effective theory for chiral
phonons.

The full inhomogeneous problem, where all parameters
$\rho,c_{ij\mu\nu},\eta_{ij\mu\nu}$ are arbitrary functions of position, can
be recast into the form of Eqs. (1,2) and can be solved explicitly following
the procedure of Secs. II, III. We leave the details of calculations and
explicit expressions to Appendix H.

Here in the main text, we restrict ourselves to the spatial dependence
described in the previous subsection. The details are left to Appendix F, but
it suffices to then study the homogeneous problem, with position-independent
parameters and matrices $\mathsf{H}^{\rm h},\mathsf{\Gamma}^{\rm h}$
($\Phi=(u,\Pi/\rho)$, with $\Pi_{i}$ the canonically conjugate momentum to
$u_{i}$, see Eq. (136)) which in the case of the Lagrangian Eq. (66) take the
simple forms

$\displaystyle{\sf\Gamma}^{\rm h}=\frac{i}{\hbar\rho}\begin{bmatrix}0&1\\\
-1&-\frac{2\eta_{\mu\nu}}{\hbar^{2}\rho}p_{\mu}p_{\nu}\end{bmatrix},\quad{\sf
H}^{\rm h}=\begin{bmatrix}\frac{1}{\hbar^{2}}c_{\mu\nu}p_{\mu}p_{\nu}&0\\\
0&\rho\end{bmatrix},$ (67)

where $(c_{\mu\nu})_{ij}\equiv c_{ij\mu\nu}$,
$(\eta_{\mu\nu})_{ij}\equiv\eta_{ij\mu\nu}$, and so each block in Eq. (67) is
a $3\times 3$ matrix, and $\mathsf{K}^{\rm h}=\hbar\,\mathsf{\Gamma}^{\rm
h}\mathsf{H}^{\rm h}$. For concreteness, we choose the following expressions
for the Hall viscosity and bulk modulus:

$\displaystyle\eta_{ij\mu\nu}$
$\displaystyle=\epsilon^{ijz}\delta_{\mu\nu}\left(\eta_{1}\overline{\delta}_{\nu
z}+\eta_{2}\delta_{\nu z}\right),$ $\displaystyle c_{ij\mu\nu}$
$\displaystyle=c_{1}\delta_{ij}\delta_{\mu\nu}+\delta_{i\mu}\delta_{j\nu}\left(c_{2}\overline{\delta}_{iz}\overline{\delta}_{jz}+c_{3}\delta_{iz}\delta_{jz}\right)$
$\displaystyle+\delta_{\mu\nu}\left(c_{4}\overline{\delta}_{\nu
z}\delta_{iz}\delta_{jz}+c_{5}\delta_{\nu
z}\overline{\delta}_{iz}\overline{\delta}_{jz}\right),$ (68)

where $\epsilon^{ijl}$ is the 3d Levi-Civita tensor,
$\overline{\delta}_{ab}\equiv 1-\delta_{ab}$, and
$\eta_{1},\eta_{2},c_{a},a=1..5$, are real parameters. This makes
$\mathcal{L}_{\rm ph}$ the most general possible Lagrangian preserving $O(2)$
symmetry around the $\mathbf{\hat{z}}$ axis, except that, for simplicity, we
do not include terms which couple the $u^{z}$ to the $(u^{x},u^{y})$
components, even if they preserve $O(2)$ symmetry. This allows us to reduce
the problem to be effectively a $4\times 4$ problem in matrix space.

We leave details of the solution of the diagonalization problem to the
appendix, and show in Figure 1 the position-resolved momentum-integrated
energy currents for a given profile of $g$ and the parameter values provided
in Table 1.

$T_{0}$ | $\delta T$ | $\rho$ | $c_{1}$ | $c_{2}$ | $c_{3}$ | $c_{4}$ | $c_{5}$ | $\eta_{1}$ | $\eta_{2}$ | $\xi$ | $L$
---|---|---|---|---|---|---|---|---|---|---|---
${\mathtt{1.0}}$ | $\mathtt{0.8}$ | $\mathtt{1.0}$ | $\mathtt{1.5}$ | $\mathtt{1.2}$ | $\mathtt{0.8}$ | $\mathtt{0.4}$ | $\mathtt{0.0}$ | $\mathtt{-0.25}$ | $\mathtt{-0.15}$ | $\mathtt{1.0}$ | $\mathtt{10}$
Table 1: Values used for the parameters defined in Eqs. (66,68) and in the
caption of Figure 1.

As is evident in Figure 1, the magnetization current
$\mathcal{J}_{X_{y}}^{(2)}(x)$ is not zero in the bulk (the exaggerated
temperature gradient we used is made to render this visible), but its integral
over a slice of material vanishes exactly. Meanwhile, the chiral current
$\mathcal{J}_{X_{y}}^{(1)}(x)$ is exactly zero in the bulk, but the
counterpropagating edge densities yield a finite integral, i.e. a nonzero Hall
conductivity.

##  V Conclusion

In this paper, we presented a systematic derivation of the quantum kinetic
equation in the full phase space for any quadratic Hamiltonian of bosonic
fields. Translation invariance is not assumed, and boundaries may be included.
The treatment provides a derivation of how the quasi-classical distribution
function descends from the full quantum density matrix, and how the evolution
of these quantities is related to single-particle characteristics such as the
spectrum, Berry curvature, and energy magnetization. To obtain these results,
we follow the exact statistical formulation of quantum dynamics using the
Moyal product in phase space, and thereby obtain exact equations that can be
expanded systematically order by order in a semiclassical parameter such as
the smallness of spatial gradients.

From this procedure, in the leading approximation we recover the full
Boltzmann equation and provide a self-contained and exact derivation of the
intrinsic thermal Hall effect of bosons. Moreover, the procedure allows a
spatially resolved examination of the associated thermal currents, treating
boundaries and the spatial profile they generate. For the intrinsic thermal
Hall effect, this method clarifies the separation of transport and
“magnetization” currents transparently, and without the need to introduce any
artificial gravitational field which is required in a Kubo formulation
following Luttinger Luttinger (1964).

Aside from pedagogical value, we believe this formalism will be useful in
simplifying and clarifying future calculations. It allows a direct application
to intrinsically inhomogeneous systems such as skyrmionic textures in chiral
magnets, which will be addressed in a forthcoming publication. It also offers
a natural extension to higher orders in the semi-classical expansion, in which
deviations from Boltzmann transport become evident. This is important for
example in non-linear transport effects (i.e. beyond linear order in applied
electric fields and thermal gradients), since such effects are indeed higher
order in gradients. Existing theory of such effects reveals a role for
geometric properties of bands such as Berry curvature and quantum metric, but
because it is based on Boltzmann transport which is only formally exact up to
_first_ order in phase space gradients, the theory is incomplete and should be
revisited.

A natural important future direction is the inclusion of scattering and
interactions into the present formalism. Scattering can be straightforwardly
added in an _ad hoc_ manner, through a collision integral, such as that
provided in Refs. Mangeolle  _et al._ (2022a, b). It would also be interesting
to derive it rigorously from the Keldysh formalism, which is a natural
extension with some precedents Keldysh (1964); Prange and Kadanoff (1964);
Rammer (2007); Kamenev (2011). Another relevant extension is to time-dependent
Hamiltonians, such as would occur in a sliding texture like a driven spin
density wave or skyrmion lattice, which leads to Berry phase effects in
spacetime. Moreover, the same methods can be applied to fermions, as in Ref.
Wickles and Belzig (2013) for normal electrons, and extended to Bogoliubov-de
Gennes quasiparticles for the superconducting case. This potentially provides
a new way to address phenomena in topological superconductors.

###### Acknowledgements.

This project was funded by the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation program (Grant agreement
No. 853116, acronym TRANSPORT). L.B. was supported by the DOE, Office of
Science, Basic Energy Sciences under Award No. DE-FG02-08ER46524. We
acknowledge the hospitality of the KITP, where we carried out part of this
project, as one does, funded under NSF Grant NSF PHY-1748958.

## References

  * Sundaram and Niu (1999) Ganesh Sundaram and Qian Niu, “Wave-packet dynamics in slowly perturbed crystals: Gradient corrections and Berry-phase effects,” Phys. Rev. B 59, 14915–14925 (1999).
  * Qin  _et al._ (2012) Tao Qin, Jianhui Zhou, and Junren Shi, “Berry curvature and the phonon Hall effect,” Phys. Rev. B 86, 104305 (2012).
  * Matsumoto  _et al._ (2014) Ryo Matsumoto, Ryuichi Shindou, and Shuichi Murakami, “Thermal Hall effect of magnons in magnets with dipolar interaction,” Phys. Rev. B 89, 054420 (2014).
  * Karplus and Luttinger (1954) Robert Karplus and J. M. Luttinger, “Hall effect in ferromagnetics,” Phys. Rev. 95, 1154–1160 (1954).
  * Adams and Blount (1959) EN Adams and EI Blount, “Energy bands in the presence of an external force field—II: Anomalous velocities,” Journal of Physics and Chemistry of Solids 10, 286–303 (1959).
  * Haldane (2004) FDM Haldane, “Berry curvature on the Fermi surface: Anomalous Hall effect as a topological Fermi-liquid property,” Physical review letters 93, 206602 (2004).
  * Nagaosa (2006) Naoto Nagaosa, “Anomalous Hall effect – a new perspective,” Journal of the Physical Society of Japan 75, 042001–042001 (2006).
  * Sinitsyn (2007) NA Sinitsyn, “Semiclassical theories of the anomalous Hall effect,” Journal of Physics: Condensed Matter 20, 023201 (2007).
  * Mishchenko and Starykh (2014) EG Mishchenko and OA Starykh, “Equilibrium currents in chiral systems with nonzero Chern number,” Physical Review B 90, 035114 (2014).
  * Katsura  _et al._ (2010) Hosho Katsura, Naoto Nagaosa, and Patrick A. Lee, “Theory of the thermal Hall effect in quantum magnets,” Phys. Rev. Lett. 104, 066403 (2010).
  * Keldysh (1964) Leonid V Keldysh, “Diagram technique for nonequilibrium processes,” Zh. Eksp. Teor. Fiz. 47, 1515 (1964).
  * Al’Tshuler (1978) BL Al’Tshuler, “Temperature dependence of impurity conductivity of metals at low temperatures,” Soviet Journal of Experimental and Theoretical Physics 48, 670 (1978).
  * Rammer and Smith (1986) J. Rammer and H. Smith, “Quantum field-theoretical methods in transport theory of metals,” Rev. Mod. Phys. 58, 323–359 (1986).
  * Gosselin  _et al._ (2007) P. Gosselin, A. Bérard, and H. Mohrbach, “Semiclassical diagonalization of quantum Hamiltonian and equations of motion with Berry phase corrections,” The European Physical Journal B 58, 137–148 (2007).
  * Wickles and Belzig (2013) Christian Wickles and Wolfgang Belzig, “Effective quantum theories for Bloch dynamics in inhomogeneous systems with nontrivial band structure,” Phys. Rev. B 88, 045308 (2013).
  * Smrcka and Streda (1977) L Smrcka and P Streda, “Transport coefficients in strong magnetic fields,” Journal of Physics C: Solid State Physics 10, 2153 (1977).
  * 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–2359 (1997).
  * Bashan and Auerbach (2022) Noga Bashan and Assa Auerbach, “Degeneracy-projected polarization formulas for Hall-type conductivities,” Phys. Rev. Lett. 128, 036601 (2022).
  * Polkovnikov (2010) Anatoli Polkovnikov, “Phase space representation of quantum dynamics,” Annals of Physics 325, 1790–1852 (2010).
  * Mangeolle  _et al._ (2022a) Léo Mangeolle, Leon Balents, and Lucile Savary, “Phonon thermal Hall conductivity from scattering with collective fluctuations,” Phys. Rev. X 12, 041031 (2022a).
  * Mangeolle  _et al._ (2022b) Léo Mangeolle, Leon Balents, and Lucile Savary, “Thermal conductivity and theory of inelastic scattering of phonons by collective fluctuations,” Phys. Rev. B 106, 245139 (2022b).
  * Xiao  _et al._ (2005) Di Xiao, Junren Shi, and Qian Niu, “Berry phase correction to electron density of states in solids,” Phys. Rev. Lett. 95, 137204 (2005).
  * Qin  _et al._ (2011) Tao Qin, Qian Niu, and Junren Shi, “Energy magnetization and the thermal Hall effect,” Phys. Rev. Lett. 107, 236601 (2011).
  * Kapustin and Spodyneiko (2020) Anton Kapustin and Lev Spodyneiko, “Thermal Hall conductance and a relative topological invariant of gapped two-dimensional systems,” Phys. Rev. B 101, 045137 (2020).
  * Hardy (1963) Robert J. Hardy, “Energy-flux operator for a lattice,” Phys. Rev. 132, 168–177 (1963).
  * Vanderbilt (2018) David Vanderbilt,  _Berry phases in electronic structure theory: electric polarization, orbital magnetization and topological insulators_ (Cambridge University Press, 2018).
  * Luttinger (1964) J. M. Luttinger, “Theory of thermal transport coefficients,” Phys. Rev. 135, A1505–A1514 (1964).
  * Barkeshli  _et al._ (2012) Maissam Barkeshli, Suk Bum Chung, and Xiao-Liang Qi, “Dissipationless phonon Hall viscosity,” Phys. Rev. B 85, 245107 (2012).
  * (29) Menxing Ye, Lucile Savary, and Leon Balents, “Phonon Hall viscosity in magnetic insulators,” arXiv 2103.04223.
  * Saito  _et al._ (2019) Takuma Saito, Kou Misaki, Hiroaki Ishizuka, and Naoto Nagaosa, “Berry phase of phonons and thermal Hall effect in nonmagnetic insulators,” Phys. Rev. Lett. 123, 255901 (2019).
  * Prange and Kadanoff (1964) Richard E Prange and Leo P Kadanoff, “Transport theory for electron-phonon interactions in metals,” Physical Review 134, A566 (1964).
  * Rammer (2007) Jørgen Rammer,  _Quantum Field Theory of Non-equilibrium States_ (Cambridge University Press, 2007).
  * Kamenev (2011) Alex Kamenev,  _Field theory of non-equilibrium systems_ (Cambridge University Press, 2011).

##  Appendix A Properties of the Moyal Product

The Moyal, or star, product is widely used in the phase space formulation of
quantum mechanics, and is particularly convenient for semi-classical analysis.
It is defined in the main text, Eq. (7), and can be applied to functions, or,
component-wise, to tensor multiplication. We review some key useful properties
of the Moyal product here.

Even for scalar functions, it is non-commutative, reflecting the non-
commutativity of operators in quantum mechanics. It is, however, associative,

${\sf O}_{1}\star{\sf O}_{2}\star{\sf O}_{3}=({\sf O}_{1}\star{\sf
O}_{2})\star{\sf O}_{3}={\sf O}_{1}\star({\sf O}_{2}\star{\sf O}_{3}).$ (69)

Under hermitian conjugation, we have

$\left({\sf O}_{1}\star{\sf O}_{2}\right)^{\dagger}={\sf
O}_{2}^{\dagger}\star{\sf O}_{1}^{\dagger},$ (70)

where the hermitian conjugate $\dagger$ is defined in the usual matrix sense,
and the phase space argument is unchanged.

The star product is not cyclic in the trace, but the combination of the trace
and phase space ($q\in\mathbb{R}^{2d}$) integral is cyclic if ${\rm
Tr}\left[{\sf O}_{1}\star{\sf O}_{2}\right]$ is integrable:

$\int_{q}{\rm Tr}\left[{\sf O}_{1}\star{\sf O}_{2}\right]=\int_{q}{\rm
Tr}\left[{\sf O}_{2}\star{\sf O}_{1}\right].$ (71)

This implies that the phase space integral of a Moyal bracket vanishes

$\displaystyle\int_{q}{\rm Tr}\left[{\sf O}_{1}\overset{\star}{,}\,{\sf
O}_{2}\right]$ $\displaystyle=0.$ (72)

If the functions ${\sf O}_{1,2}$ are differentiable, the chain rule applies,
i.e., for $\alpha=1,..,2d$,

$\partial_{q_{\alpha}}({\sf O}_{1}\star{\sf O}_{2})=\partial_{q_{\alpha}}{\sf
O}_{1}\star{\sf O}_{2}+{\sf O}_{1}\star\partial_{q_{\alpha}}{\sf O}_{2}.$ (73)

A matrix function ${\sf O}$ is star-invertible if there exists ${\sf O}^{-1}$
such that

${\sf O}\star{\sf O}^{-1}={\sf O}^{-1}\star{\sf O}=1.$ (74)

Note that, as for ordinary matrix multiplication, the left and right star-
inverses are identical. This follows from viewing the convolution as a matrix
product in the combinated coordinate-index space, and then using the fact that
the Wigner transform of a convolution is the star product of Wigner
transforms.

If ${\sf O}_{1}$ and ${\sf O}_{2}$ are star-invertible,

$({\sf O}_{1}\star{\sf O}_{2})^{-1}={\sf O}_{2}^{-1}\star{\sf O}_{1}^{-1}$
(75)

Up to first order in $\hslash$, we have:

${\sf O}_{1}\star{\sf O}_{2}={\sf O}_{1}{\sf
O}_{2}+i\hslash\frac{\epsilon_{\alpha\beta}}{2}\partial_{\alpha}{\sf
O}_{1}\partial_{\beta}{\sf O}_{2}+O(\hslash^{2}).$ (76)

The zeroth order term is just standard matrix multiplication while the first
order one is related to the Poisson bracket:
$\epsilon_{\alpha\beta}\partial_{\alpha}{\sf O}_{1}\partial_{\beta}{\sf
O}_{2}=\partial_{X_{\mu}}{\sf O}_{1}\partial_{p_{\mu}}{\sf
O}_{2}-\partial_{p_{\mu}}{\sf O}_{1}\partial_{X_{\mu}}{\sf O}_{2}\equiv\\{{\sf
O}_{1},{\sf O}_{2}\\}_{\rm p.b.}$. Eq. (76) can be rewritten, up to first
order:

${\sf O}_{1}{\sf O}_{2}={\sf O}_{1}\star{\sf
O}_{2}-i\frac{\hslash}{2}\epsilon_{\alpha\beta}\partial_{\alpha}{\sf
O}_{1}\partial_{\beta}{\sf O}_{2}+O(\hslash^{2}).$ (77)

From the above, we can find the ordinary product of a matrix with its star-
inverse

$\displaystyle{\sf O}^{-1}{\sf O}$
$\displaystyle=1-i\frac{\hslash}{2}\epsilon_{\alpha\beta}\partial_{\alpha}{\sf
O}^{-1}\partial_{\beta}{\sf O}+O(\hslash^{2}),$ (78)

which notably is not the identity.

##  Appendix B Diagonalization of ${\sf K}$

###  B.1 Formal relation to a hermitian problem

Here we show that the star-diagonalization of $\mathsf{K}$ is related to a
hermitian star-eigenvalue problem, provided the original physical problem is
well-defined. In particular, this requires that ${\sf H}$ is positive
definite. The latter implies that there exists a positive definite matrix
${\sf H}^{1/2}$ such that

${\sf H}^{1/2}\star{\sf H}^{1/2}={\sf H}.$ (79)

We define the matrix $\breve{\sf K}$

$\breve{\sf K}\equiv\hbar\,{\sf H}^{1/2}\star{\sf\Gamma}\star{\sf
H}^{1/2}=\breve{\sf K}^{\dagger}.$ (80)

Since $\breve{\sf K}$ is hermitian, one can find a star-unitary matrix ${\sf
U}$, i.e.

${\sf U}^{\dagger}\star{\sf U}={\sf U}\star{\sf U}^{\dagger}=1,$ (81)

which star-diagonalizes $\breve{\sf K}$ according to

${\sf U}^{\dagger}\star\breve{\sf K}\star{\sf U}=\breve{\sf K}_{d},$ (82)

where $\breve{\sf K}_{d}$ is a diagonal matrix with real entries (the
‘stargenvalues’).

Plugging in Eq. (80) into Eq. (82), we obtain

$\displaystyle\breve{\sf K}_{d}$ $\displaystyle=\hbar\,{\sf
U}^{\dagger}\star{\sf H}^{1/2}\star{\sf\Gamma}\star{\sf H}^{1/2}\star{\sf U}$
$\displaystyle=\hbar\,({\sf U}^{\dagger}\star{\sf
H}^{1/2})\star{\sf\Gamma}\star{\sf H}\star({\sf H}^{-1/2}\star{\sf U}).$ (83)

Using ${\sf K}=\hbar\,{\sf\Gamma}\star{\sf H}$ we see that this has the form
of Eq. (31) if we identify

$\breve{\sf K}_{d}={\sf K}_{d},$ (84)

which, in particular, means that ${\sf K}_{d}$ is real, and

${\sf S}={\sf H}^{-1/2}\star{\sf U},\quad{\sf S}^{-1}={\sf
U}^{\dagger}\star{\sf H}^{1/2},$ (85)

up to a gauge choice, as discussed in the main text.

From Eq. (85) we can immediately show Eq. (32) by direct calculation. Eq. (32)
and Eq. (85) can be considered a choice of normalization of $\mathsf{S}$.
Independent of that normalization, we note the identity

${\sf K}^{\dagger}={\sf H}\star{\sf K}\star{\sf H}^{-1}.$ (86)

###  B.2 Diagonalization to $O(\hbar)$

Now we show how to convert the formidable star-diagonalization problem defined
by Eq. (31) and Eq. (32) into conventional linear algebra in the semi-
classical expansion. We aim here to obtain $\mathsf{K}_{d}$ and $\mathsf{S}$
up to first order in $\hslash$ (and $\mathsf{S}$ up to gauge freedom, of
course).

To do so, we specify the expansions of $\mathsf{K}$ and $\mathsf{S}$. The
former is unambiguously determined by expanding the Moyal product, regarding
$\mathsf{H}$ and $\mathsf{\Gamma}$ as $\hslash$-independent. We have

${\sf K}=\hbar\,\left({\sf\Gamma}{\sf
H}+i\hslash\frac{\epsilon_{\alpha\beta}}{2}\partial_{\alpha}{\sf\Gamma}\partial_{\beta}{\sf
H}\right)+O(\hslash^{2}),$ (87)

and so

${\sf K}={\sf K}_{0}+\hslash{\sf K}_{1}+O(\hslash^{2}),$ (88)

with

$\displaystyle{\sf K}_{0}$ $\displaystyle=\hbar\,{\sf\Gamma}{\sf H},$
$\displaystyle{\sf K}_{1}$
$\displaystyle=i\hbar\frac{\epsilon_{\alpha\beta}}{2}\partial_{\alpha}{\sf\Gamma}\partial_{\beta}{\sf
H}.$ (89)

Next, we write the expansion of $\mathsf{S}$ according to the convenient form

$\displaystyle\mathsf{S}$
$\displaystyle=\mathsf{S}_{0}\left(1+\hslash\tilde{\mathsf{S}}_{1}\right)+O(\hslash^{2}).$
(90)

Here the first two terms in the expansion are specified by $\mathsf{S}_{0}$
and $\tilde{\mathsf{S}}_{1}$, both $O(\hslash^{0})$ a priori, which are to be
determined by enforcing the star-diagonalization conditions up to first order
in $\hslash$. The hermitian conjugate gives

$\displaystyle\mathsf{S}^{\dagger}$
$\displaystyle=\left(1+\hslash\tilde{\mathsf{S}}^{\dagger}_{1}\right)\mathsf{S}_{0}^{\dagger}+O(\hslash^{2}).$
(91)

Inserting Eq. (90) and Eq. (91) into Eq. (32), and equating terms at zeroth
and first order in $\hslash$ gives the conditions

$\displaystyle\mathsf{S}_{0}^{\dagger}\mathsf{H}\mathsf{S}_{0}$
$\displaystyle=1,$ (92a)
$\displaystyle\tilde{\mathsf{S}}_{1}^{\dagger}+\tilde{\mathsf{S}}_{1}$
$\displaystyle=-\frac{i\epsilon_{\alpha\beta}}{2}\left(\partial_{\alpha}\mathsf{S}_{0}^{\dagger}\mathsf{H}\partial_{\beta}\mathsf{S}_{0}+\mathsf{S}_{0}^{\dagger}\partial_{\beta}\mathsf{H}\mathsf{S}_{0}+\mathsf{S}_{0}^{\dagger}\partial_{\alpha}\mathsf{H}\partial_{\beta}\mathsf{S}_{0}\right).$
(92b)

We see that the hermitian part of $\tilde{\mathsf{S}}_{1}$ is fixed by the
second condition, once the zeroth order term is known.

To proceed, we also need $\mathsf{S}^{-1}$, which we obtain using Eq. (32),
Eq. (91) and expanding consistently to $O(\hslash)$:

$\displaystyle\mathsf{S}^{-1}$
$\displaystyle=\mathsf{S}^{\dagger}\star\mathsf{H}$
$\displaystyle=\left(1+\hslash\left(\tilde{\mathsf{S}}_{1}^{\dagger}+\frac{i\epsilon_{\alpha\beta}}{2}\partial_{\alpha}\mathsf{S}_{0}^{\dagger}\partial_{\beta}\mathsf{H}\mathsf{S}_{0}\right)\right)\mathsf{S}_{0}^{-1}+O(\hslash^{2})$
(93)

(take care that $\mathsf{S}^{-1}$ is the star-inverse not the ordinary matrix
inverse). We also used the zeroth order condition

$\displaystyle\mathsf{S}_{0}^{-1}$
$\displaystyle=\mathsf{S}_{0}^{\dagger}\mathsf{H},$ (94)

which follows from Eq. (92a) (here $\mathsf{S}_{0}^{-1}$ is the
$O(\hslash^{0})$ term in the expansion of $\mathsf{S}^{-1}$ as well as the
inverse of $\mathsf{S}_{0}$ for the _usual_ matrix product). Then one can use
Eq. (92b) to eliminate $\tilde{\mathsf{S}}_{1}^{\dagger}$ in favor of
$\tilde{\mathsf{S}}_{1}$ within Eq. (B.2), which yields

$\displaystyle\mathsf{S}^{-1}$
$\displaystyle=\left(1-\hslash\left(\tilde{\mathsf{S}}_{1}+\frac{i\epsilon_{\alpha\beta}}{2}\partial_{\alpha}\mathsf{S}_{0}^{-1}\partial_{\beta}\mathsf{S}_{0}\right)\right)\mathsf{S}_{0}^{-1}+O(\hslash^{2})$
$\displaystyle=\left(1-\hslash\left(\tilde{\mathsf{S}}_{1}-\frac{i\epsilon_{\alpha\beta}}{2}\mathsf{\Lambda}_{\alpha}\mathsf{\Lambda}_{\beta}\right)\right)\mathsf{S}_{0}^{-1}+O(\hslash^{2}).$
(95)

Eq. (B.2) follows from the constraint on ${\sf S}$ from the normalization
condition Eq. (32), but so far we have not used the actual diagonalization
condition in Eq. (31). Using Eq. (90) and Eq. (B.2) in Eq. (31) and collecting
terms to $O(\hslash)$ we find

$\displaystyle{\sf K}_{d}=$ $\displaystyle\,{\sf K}_{0,d}+\hslash\bigg{(}{\sf
S}_{0}^{-1}{\sf K}_{1}{\sf S}_{0}+\left[{\sf K}_{0,d},\tilde{\sf
S}_{1}\right]$
$\displaystyle\,-i\frac{\epsilon_{\alpha\beta}}{2}\left\\{{\sf\Lambda}_{\alpha},\partial_{\beta}{\sf
K}_{0,d}-{\sf K}_{0,d}{\sf\Lambda}_{\beta}\right\\}\bigg{)}+O(\hslash^{2}),$
(96)

where

$\displaystyle\mathsf{K}_{0,d}$
$\displaystyle=\mathsf{S}_{0}^{-1}\mathsf{K}_{0}\mathsf{S}_{0}.$ (97)

The $O(\hslash^{0})$ conditions, Eqs. (97,92a), taken together, define a
standard (but non-hermitian) eigenvalue problem, where we demand that
$\mathsf{K}_{0,d}$ is diagonal. By a classical limit, i.e. a regular matrix
product version, of the argument in Sec. B.1, it is straightforward to show
that a solution always exists and that $\mathsf{K}_{0,d}$ is not only diagonal
but hermitian, i.e. real diagonal.

To achieve the star-diagonalization to $O(\hslash)$, we must require that the
term in parenthesis in Eq. (B.2) be diagonal, i.e. that

$\displaystyle\bigg{(}{\sf S}_{0}^{-1}{\sf K}_{1}{\sf S}_{0}+\left[{\sf
K}_{0,d},\tilde{\sf S}_{1}\right]$ (98)
$\displaystyle\;-i\frac{\epsilon_{\alpha\beta}}{2}\left\\{{\sf\Lambda}_{\alpha},\partial_{\beta}{\sf
K}_{0,d}-{\sf K}_{0,d}{\sf\Lambda}_{\beta}\right\\}\bigg{)}_{ab}=0,$
$\displaystyle a\neq b.$

This is a condition on $\tilde{\mathsf{S}}_{1}$, and while an explicit form
for the latter can be found, we will actually not need it for our
calculations, so we suffice to say that Eq. (98) is well-defined and
consistent with Eq. (92b) and the hermiticity of $\mathsf{K}_{d}$.

Since $\mathsf{K}_{d}$ is diagonal and real, we can take the real part of the
diagonal projection of Eq. (B.2), which gives

$\displaystyle{\sf K}_{d}$ $\displaystyle={\sf K}_{0,d}+\hslash{\rm
Re}\left[\left({\sf S}_{0}^{-1}{\sf K}_{1}{\sf S}_{0}\right)^{({\rm
d})}\right]$ (99) $\displaystyle+\hslash\epsilon_{\alpha\beta}\left({\sf
A}_{\alpha}\partial_{\beta}{\sf K}_{0,d}-\frac{1}{4}{\rm
Im}\left\\{\mathsf{\Lambda}_{\beta},\left[\mathsf{\Lambda}_{\alpha},\mathsf{K}_{0,d}\right]\right\\}^{\rm(d)}\right)+O(\hslash^{2}).$

Transforming to the gauge-invariant function $\underline{\mathsf{K}}_{d}$ then
immediately gives Eq. (19) in the main text, noting that ${\rm
Im}\left\\{\mathsf{\Lambda}_{\beta},\left[\mathsf{\Lambda}_{\alpha},\mathsf{K}_{0,d}\right]\right\\}^{\rm(d)}$
is antisymmetric in $\alpha\leftrightarrow\beta$, and so $\sf
M_{\alpha\beta}=-\sf M_{\beta\alpha}$.

###  B.3 Symmetry of the spectrum

From their definitions, we know that
$\mathsf{\Gamma}(X,p)=-\mathsf{\Gamma}(X,-p)^{*}$ and
$\mathsf{H}(X,p)=\mathsf{H}(X,-p)^{*}$. This implies that the dynamical matrix
obeys $\mathsf{K}(X,p)=-\mathsf{K}(X,-p)^{*}$. The latter induces relations
amongst the star-eigenvalues of $\mathsf{K}_{d}$ at reversed momenta. Taking
the complex conjugate of Eq. (31) evaluated at momentum $-p$, and using the
fact that $\mathsf{K}_{d}$ is real, we deduce that

$\displaystyle\mathsf{K}_{d}(X,-p)=-\left[\mathsf{S}^{-1}(X,-p)\right]^{*}\star\mathsf{K}(X,p)\star\left[\mathsf{S}(X,-p)\right]^{*}.$
(100)

Taking the complex conjugate of Eq. (32) at momentum $-p$ gives

$\displaystyle\mathsf{S}(X,-p)^{\top}\star\mathsf{H}(X,p)\star\mathsf{S}(X,-p)^{*}=1.$
(101)

These two equations are solved by using again Eqs. (31,32), by

$\displaystyle\mathsf{S}(X,-p)=\mathsf{S}(X,p)^{*},$
$\displaystyle\mathsf{K}_{d}(X,-p)=-\mathsf{K}_{d}(X,p),$ (102)

_up to a gauge transformation_ , and up to permutations of the diagonal
entries of $\mathsf{K}_{d}$. We see that the negative energies at momentum $p$
appear with the opposite sign as positive energies at momentum $-p$, and vice
versa. Thus the modes of the spectrum of $\mathsf{K}$ come in positive (the
“physical” ones) and negative energy pairs,
$(\varepsilon_{a}(X,p),\varepsilon_{\bar{a}}(X,p)\equiv-\varepsilon_{a}(X,-p))$.
This defines the notation $\bar{a}$. Moreover, we can deduce that

$\displaystyle[\mathsf{\Omega}_{\alpha\beta}(X,-p)]_{\bar{a}\bar{a}}$
$\displaystyle=-[\mathsf{\Omega}_{\alpha\beta}(X,p)]_{aa},$
$\displaystyle[\mathsf{M}_{\alpha\beta}(X,-p)]_{\bar{a}\bar{a}}$
$\displaystyle=+[\mathsf{M}_{\alpha\beta}(X,p)]_{aa}.$ (103)

We make use of these last two properties in Sec. IV.1 and App. F.4.

##  Appendix C Diagonalization of $\mathcal{J}_{\alpha}$

Here we provide details of the calculations leading from Eq. (30) to Eqs.
(60,61a,61b) in Sec. III.6.2.

Beginning with Eq. (30), we use Eq. (57) and Eq. (59) to obtain

$\displaystyle\mathcal{J}_{\alpha}$
$\displaystyle=\frac{\upsilon}{2}\epsilon_{\alpha\beta}{\rm Re}{\rm
Tr}\Big{[}\left({\sf S}\star{\sf Q}_{\beta}\star{\sf S}^{-1}\right)\left({\sf
S}\star{\sf F}_{d}\star{\sf S}^{-1}\right)\Big{]}+O(\upsilon\hslash^{2}).$
(104)

We now use Eq. (77) with ${\sf O}_{1}={\sf S}\star{\sf Q}_{\beta}\star{\sf
S}^{-1}$ and ${\sf O}_{2}={\sf S}\star{\sf F}_{d}\star{\sf S}^{-1}$ to
“restore” a star product and as a second step keep only the zeroth order term
of the expansion of the star products in terms which are already of order
$\hslash$ (which amounts to turning the star products into regular matrix
products):

$\displaystyle\mathcal{J}_{\alpha}$
$\displaystyle=\frac{\upsilon}{2}\epsilon_{\alpha\beta}{\rm Re}{\rm
Tr}\Big{[}{\sf S}\star{\sf Q}_{\beta}\star{\sf F}_{d}\star{\sf
S}^{-1}\Big{]}+\frac{\upsilon\hslash}{4}\epsilon_{\alpha\beta}\epsilon_{\gamma\lambda}{\rm
Im}{\rm Tr}\left[\partial_{\gamma}\left({\sf S}{\sf Q}_{\beta}{\sf
S}^{-1}\right)\partial_{\lambda}\left({\sf S}{\sf F}_{d}{\sf
S}^{-1}\right)\right]+O(\upsilon\hslash^{2}).$ (105)

Now we use Eq. (53) with ${\sf O}={\sf Q}_{\beta}\star{\sf F}_{d}$, extract a
$\partial_{\gamma}$ derivative using
$\epsilon_{\gamma\lambda}\partial^{2}_{\gamma\lambda}=0$, expand the
$\partial_{\lambda}$ derivative using the chain rule, and use the cyclicity of
the trace:

$\displaystyle\mathcal{J}_{\alpha}$
$\displaystyle=\frac{\upsilon}{2}\epsilon_{\alpha\beta}{\rm Re}{\rm
Tr}\Big{[}{\sf Q}_{\beta}\star{\sf
F}_{d}\Big{]}-\frac{\upsilon}{2}\hslash\epsilon_{\alpha\beta}\epsilon_{\gamma\lambda}\partial_{\lambda}{\rm
Im}{\rm Tr}[{\sf S}^{-1}\partial_{\gamma}{\sf S}\left({\sf Q}_{\beta}{\sf
F}_{d}\right)]$
$\displaystyle\qquad+\frac{\upsilon\hslash}{4}\epsilon_{\alpha\beta}\epsilon_{\gamma\lambda}\partial_{\gamma}{\rm
Im}{\rm Tr}\left[{\sf S}{\sf Q}_{\beta}{\sf
S}^{-1}\left(\partial_{\lambda}{\sf S}{\sf F}_{d}{\sf S}^{-1}+{\sf
S}\partial_{\lambda}{\sf F}_{d}{\sf S}^{-1}+{\sf S}{\sf
F}_{d}\partial_{\lambda}{\sf S}^{-1}\right)\right]+O(\upsilon\hslash^{2})$
$\displaystyle=\frac{\upsilon}{2}\epsilon_{\alpha\beta}{\rm Re}{\rm
Tr}\Big{[}{\sf Q}_{\beta}\star{\sf
F}_{d}\Big{]}+\frac{\upsilon\hslash}{4}\epsilon_{\alpha\beta}\epsilon_{\gamma\lambda}\partial_{\gamma}{\rm
Im}{\rm Tr}\left[{\sf Q}_{\beta}\left(\left\\{{\sf
S}^{-1}\partial_{\lambda}{\sf S},{\sf F}_{d}\right\\}+\partial_{\lambda}{\sf
F}_{d}\right)\right]+O(\upsilon\hslash^{2}).$ (106)

Then, plugging in the expression for ${\sf Q}_{\beta}$ in Eq. (58), we get

$\displaystyle\mathcal{J}_{\alpha}$
$\displaystyle=\frac{\upsilon}{2}\epsilon_{\alpha\beta}\left({\rm Re}{\rm
Tr}\Big{[}\partial_{\beta}{\sf K}_{d}\star{\sf F}_{d}\Big{]}+{\rm Re}{\rm
Tr}\Big{[}[{\sf S}^{-1}\star\partial_{\beta}{\sf S}\,\overset{\star}{,}\,{\sf
K}_{d}]\star{\sf F}_{d}\Big{]}\right)$
$\displaystyle\quad+\frac{\upsilon\hslash}{4}\epsilon_{\alpha\beta}\epsilon_{\gamma\lambda}\partial_{\gamma}{\rm
Im}{\rm Tr}\left[(\partial_{\beta}{\sf K}_{d}+[{\sf
S}^{-1}\partial_{\beta}{\sf S},{\sf K}_{d}])\left(\partial_{\lambda}{\sf
F}_{d}+\left\\{{\sf S}^{-1}\partial_{\lambda}{\sf S},{\sf
F}_{d}\right\\}\right)\right]+O(\upsilon\hslash^{2}).$ (107)

In turn, expanding the star products yields

$\displaystyle\mathcal{J}_{\alpha}$
$\displaystyle=\frac{\upsilon}{2}\epsilon_{\alpha\beta}\left({\rm Re}{\rm
Tr}\Big{[}\partial_{\beta}{\sf K}_{d}{\sf F}_{d}\Big{]}+{\rm Re}{\rm
Tr}\Big{[}[{\sf S}^{-1}\partial_{\beta}{\sf S},{\sf K}_{d}]{\sf
F}_{d}\Big{]}\right)$
$\displaystyle\quad-\frac{\upsilon\hslash}{4}\epsilon_{\alpha\beta}\epsilon_{\gamma\lambda}{\rm
Im}{\rm Tr}\left[\partial^{2}_{\beta\gamma}{\sf K}_{d}\partial_{\lambda}{\sf
F}_{d}+\partial_{\gamma}[{\sf S}^{-1}\partial_{\beta}{\sf S},{\sf
K}_{d}]\partial_{\lambda}{\sf F}_{d}+\\{\partial_{\gamma}({\sf
S}^{-1}\partial_{\beta}{\sf S}),\partial_{\lambda}{\sf K}_{d}\\}{\sf
F}_{d}+[\partial_{\gamma}{\sf S}^{-1}\partial^{2}_{\beta\lambda}{\sf S},{\sf
K}_{d}]{\sf F}_{d}\right]$
$\displaystyle\quad+\frac{\upsilon\hslash}{4}\epsilon_{\alpha\beta}\epsilon_{\gamma\lambda}\partial_{\gamma}{\rm
Im}{\rm Tr}\left[(\partial_{\beta}{\sf K}_{d}+[{\sf
S}^{-1}\partial_{\beta}{\sf S},{\sf K}_{d}])\left(\partial_{\lambda}{\sf
F}_{d}+\left\\{{\sf S}^{-1}\partial_{\lambda}{\sf S},{\sf
F}_{d}\right\\}\right)\right]+O(\upsilon\hslash^{2}).$ (108)

The contributions from $\partial^{2}_{\beta\gamma}{\sf
K}_{d}\partial_{\lambda}{\sf F}_{d}$ and $\partial_{\beta}{\sf
K}_{d}\partial_{\lambda}{\sf F}_{d}$ vanish because the latter are real, and
those involving commutators, except for $[{\sf S}^{-1}\partial_{\beta}{\sf
S},{\sf K}_{d}]\left\\{{\sf S}^{-1}\partial_{\lambda}{\sf S},{\sf
F}_{d}\right\\}$, also vanish because the trace selects only the diagonal
elements. One is then left with

$\displaystyle\mathcal{J}_{\alpha}$
$\displaystyle=\frac{\upsilon}{2}\epsilon_{\alpha\beta}{\rm Re}{\rm
Tr}\Big{[}\partial_{\beta}{\sf K}_{d}{\sf
F}_{d}\Big{]}-\frac{\upsilon\hslash}{4}\epsilon_{\alpha\beta}\epsilon_{\gamma\lambda}{\rm
Im}{\rm
Tr}\Big{[}\\{\partial_{\gamma}\mathsf{\Lambda}_{\beta},\partial_{\lambda}{\sf
K}_{d}\\}{\sf F}_{d}\Big{]}$
$\displaystyle\quad+\frac{\upsilon\hslash}{4}\epsilon_{\alpha\beta}\epsilon_{\gamma\lambda}\partial_{\gamma}{\rm
Im}{\rm Tr}\Big{[}[\mathsf{\Lambda}_{\beta},{\sf
K}_{d}]\left\\{\mathsf{\Lambda}_{\lambda},{\sf
F}_{d}\right\\}+\partial_{\beta}{\sf
K}_{d}\left\\{\mathsf{\Lambda}_{\lambda},{\sf
F}_{d}\right\\}\Big{]}+O(\upsilon\hslash^{2}).$ (109)

Here we replaced
$\mathsf{S}^{-1}\partial_{\rho}\mathsf{S}\rightarrow\mathsf{\Lambda}_{\rho}$
which is correct to this order.

In almost all terms in Eq. (C), the imaginary diagonal part of
$\mathsf{\Lambda}_{\rho}$, which is simply $\mathsf{A}_{\rho}$, is selected by
the trace against other diagonal real functions. The one term not of this form
can be rewritten using the cyclicity of the trace as

$\displaystyle{\rm Tr}\left([\mathsf{\Lambda}_{\alpha},{\sf
K}_{d}]\\{\mathsf{\Lambda}_{\beta},{\sf F}_{d}\\}\right)$ $\displaystyle={\rm
Tr}\left(\left\\{\mathsf{\Lambda}_{\beta},\left[\mathsf{\Lambda}_{\alpha},{\sf
K}_{d}\right]\right\\}{\sf F}_{d}\right)=2{\rm
Tr}\left(\mathsf{M}_{\alpha\beta}{\sf F}_{d}\right),$ (110)

where $\mathsf{M}_{\alpha\beta}$ defined in Eq. (18). Consequently,

$\displaystyle\mathcal{J}_{\alpha}$
$\displaystyle=\frac{\upsilon}{2}\epsilon_{\alpha\beta}\left({\rm
Tr}\Big{[}\partial_{\beta}{\sf K}_{d}{\sf
F}_{d}\Big{]}+\frac{\hslash}{2}\epsilon_{\gamma\lambda}\left({\rm
Tr}\left[-\\{\partial_{\gamma}{\sf A}_{\beta},\partial_{\lambda}{\sf
K}_{d}\\}{\sf F}_{d}+\partial_{\gamma}(\partial_{\beta}{\sf K}_{d}\left\\{{\sf
A}_{\lambda},{\sf F}_{d}\right\\})\right]+2\partial_{\gamma}{\rm
Tr}\left[\mathsf{M}_{\beta\lambda}{\sf
F}_{d}\right]\right)\right)+O(\upsilon\hslash^{2})$
$\displaystyle=\frac{\upsilon}{2}\epsilon_{\alpha\beta}\left({\rm
Tr}\Big{[}\partial_{\beta}{\sf K}_{d}{\sf
F}_{d}\Big{]}+\hslash\epsilon_{\gamma\lambda}\left({\rm
Tr}\left[\partial_{\gamma}({\sf A}_{\lambda}\partial_{\beta}{\sf K}_{d}-{\sf
A}_{\beta}\partial_{\lambda}{\sf K}_{d}){\sf F}_{d}+{\sf
A}_{\lambda}\partial_{\beta}{\sf K}_{d}\partial_{\gamma}{\sf
F}_{d}\right]+\partial_{\gamma}{\rm Tr}\left[\mathsf{M}_{\beta\lambda}{\sf
F}_{d}\right]\right)\right)+O(\upsilon\hslash^{2}).$ (111)

From here, it is straightforward to separate the term containing
$\mathsf{M}_{\beta\lambda}$ from the others, and with some rearrangement to
first order in $\hslash$, obtain Eqs. (60,61a,61b) in Sec. III.6.2.

##  Appendix D Jacobian

When expressed in gauge invariant band diagonal form, physical quantities like
the energy and current density display explicit dependence on the “Jacobian”
factor defined in Eq. (16). Here we show that this can indeed be regarded as a
Jacobian, i.e. a determinant defining the measure for a change of variables in
phase space.

Consider the change of variables in Eq. (10) from $\underline{\smash{q}}$ to
$q$, which induces the change in measure

$\displaystyle d^{D}\underline{\smash{q}}$
$\displaystyle=d^{D}q\,\left|\textrm{det}_{D}\left[\left(\frac{\partial\underline{\smash{q}}_{\alpha}}{\partial
q_{\beta}}\right)_{\alpha\beta}\right]\right|.$ (112)

Here the measure is defined in phase space, with total dimension $D=2d$ twice
the spatial dimension, and the subscript $D$ on the determinant is meant to
indicate that the determinant is in this space. One should avoid confusion
with band space, which plays a trivial role here because the change of
variables is diagonal (since $\mathsf{A}$ in Eq. (10) is diagonal). The
determinant can be regarded as a scalar for a given band, or a diagonal matrix
in band space, evaluated separately for each band.

The Jacobian of the transformation is just this determinant. Using the
explicit form of Eq. (10),

$\displaystyle\mathfrak{J}$ $\displaystyle={\rm
det}_{D}\left[\left(\frac{\partial\underline{\smash{q}}_{\alpha}}{\partial
q_{\beta}}\right)_{\alpha\beta}\right]$ (113) $\displaystyle={\rm
det}_{D}\left[\left(\delta_{\alpha\beta}+\hslash\epsilon_{\alpha\gamma}\partial_{\beta}{\sf
A}_{\gamma}\right)_{\alpha\beta}\right]+O(\hslash^{2}).$

Note that we have dropped the absolute value of the determinant, because the
Jacobian is perturbatively close to $1$ in the semi-classical expansion and so
has definite sign. We have, with ${\rm Tr}_{D}$ the trace in phase space,

$\displaystyle\mathfrak{J}$ $\displaystyle=e^{{\rm
Tr}_{D}[\ln(\left(\delta_{\alpha\beta}+\hslash\epsilon_{\alpha\gamma}\partial_{\beta}{\sf
A}_{\gamma}\right)_{\alpha\beta})]}+O(\hslash^{2})$ $\displaystyle=e^{{\rm
Tr}_{D}[\left(\hslash\epsilon_{\alpha\gamma}\partial_{\beta}{\sf
A}_{\gamma}\right)_{\alpha\beta}]}+O(\hslash^{2})$ $\displaystyle=1+{\rm
Tr}_{D}[\left(\hslash\epsilon_{\alpha\gamma}\partial_{\beta}{\sf
A}_{\gamma}\right)_{\alpha\beta}]+O(\hslash^{2})$
$\displaystyle=1+\hslash\epsilon_{\alpha\beta}\partial_{\alpha}{\sf
A}_{\beta}+O(\hslash^{2}).$ (114)

This is precisely the definition in Eq. (16).

##  Appendix E Explicit separation into position and momentum coordinates

In the main text, we presented the quasiparticle equations of motion and
current density in compact forms combining positions and momenta into a single
phase space coordinate. To connect to more familiar forms, we expand these
explicitly here.

Writing $q_{X_{\mu}}=X_{\mu}$ and $q_{p_{\mu}}=p_{\mu}$, we obtain from Eq.
(13) the two analogs of Hamilton’s equations

$\displaystyle\partial_{t}{\sf X}_{\mu}$
$\displaystyle=\partial_{p_{\mu}}\underline{\sf
K}_{d}-\hslash\left({\sf\Omega}_{p_{\mu}p_{\nu}}\partial_{X_{\nu}}\underline{\sf
K}_{d}-{\sf\Omega}_{p_{\mu}X_{\nu}}\partial_{p_{\nu}}\underline{\sf
K}_{d}\right)+O(\hslash^{2}),$ (115a) $\displaystyle\partial_{t}{\sf p}_{\mu}$
$\displaystyle=-\partial_{X_{\mu}}\underline{\sf
K}_{d}+\hslash\left({\sf\Omega}_{X_{\mu}p_{\nu}}\partial_{X_{\nu}}\underline{\sf
K}_{d}-{\sf\Omega}_{X_{\mu}X_{\nu}}\partial_{p_{\nu}}\underline{\sf
K}_{d}\right)+O(\hslash^{2}).$ (115b)

This defines a velocity $\mathsf{v}_{\mu}\equiv\partial_{t}{\sf X}_{\mu}$ and
a force $\mathsf{f}_{\mu}\equiv\partial_{t}{\sf p}_{\mu}$.

We can similarly inspect the real space and momentum space components of
$\mathcal{J}_{\alpha}$. Taking $\alpha=X_{\mu}$, we obtain the quantity
describing the flow of energy in the $X_{\mu}$ direction. This is just the
momentum-resolved energy current. We have

$\displaystyle\mathcal{J}_{X_{\mu}}$ $\displaystyle=\frac{\upsilon}{2}{\rm
Tr}\left[\,\mathfrak{J}\,\partial_{t}{\sf X}_{\mu}\,\underline{\sf
F}_{d}+\hslash\left(\partial_{X_{\nu}}\left(\mathsf{M}_{p_{\mu}p_{\nu}}\underline{\sf
F}_{d}\right)-\partial_{p_{\nu}}\left(\mathsf{M}_{p_{\mu}X_{\nu}}\underline{\sf
F}_{d}\right)\right)\right]+O(\upsilon\hslash^{2})$
$\displaystyle=\frac{\upsilon}{2}{\rm
Tr}\left[\,\mathfrak{J}\,\partial_{t}{\sf X}_{\mu}\,\underline{\sf
F}_{d}+\hslash\epsilon^{\mu\nu\lambda}\left(\partial_{X_{\nu}}\mathfrak{M}_{pp}^{\lambda}-\partial_{p_{\nu}}\mathfrak{M}_{pX}^{\lambda}\right)\right]+O(\upsilon\hslash^{2}),$
(116)

where $\epsilon^{\mu\nu\lambda}$ (note here we use superscripts and not
subscripts for the indices) is the usual 3d Levi-Civita tensor, and

$\mathfrak{M}_{pp}^{\lambda}=\frac{1}{2}\epsilon^{\mu\nu\lambda}\mathsf{M}_{p_{\mu}p_{\nu}}\underline{\sf
F}_{d},\qquad\mathfrak{M}_{pX}^{\lambda}=\frac{1}{2}\epsilon^{\mu\nu\lambda}\mathsf{M}_{p_{\mu}X_{\nu}}\underline{\sf
F}_{d},$ (117)

because $\mathsf{M}_{\alpha\beta}=-\mathsf{M}_{\beta\alpha}$. The total energy
current at position $X_{\mu}$ is obtained by integrating
$\mathcal{J}_{X_{\mu}}(X_{\mu},p_{\mu})$ over $p$. Then clearly the
$\partial_{p_{\nu}}\mathfrak{M}_{pX}^{\lambda}$ term drops. Moreover, if one
computes the flux of $\mathcal{J}_{X_{\mu}}$ through an open 2d surface at
whose boundary $\mathfrak{M}_{pp}^{\lambda}$ vanishes, the
$\partial_{X_{\nu}}\mathfrak{M}_{pp}^{\lambda}$ term drops.

Taking now $\alpha=p_{\mu}$ describes the flow of energy in momentum space.
This is something like an “energy force” density. It is

$\displaystyle\mathcal{J}_{p_{\mu}}$ $\displaystyle=\frac{\upsilon}{2}{\rm
Tr}\left[\,\mathfrak{J}\,\partial_{t}{\sf p}_{\mu}\,\underline{\sf
F}_{d}-\hslash\left(\partial_{X_{\nu}}\left(\mathsf{M}_{X_{\mu}p_{\nu}}\underline{\sf
F}_{d}\right)-\partial_{p_{\nu}}\left(\mathsf{M}_{X_{\mu}X_{\nu}}\underline{\sf
F}_{d}\right)\right)\right]+O(\upsilon\hslash^{2})$
$\displaystyle=\frac{\upsilon}{2}{\rm
Tr}\left[\,\mathfrak{J}\,\partial_{t}{\sf p}_{\mu}\,\underline{\sf
F}_{d}-\hslash\epsilon^{\mu\nu\lambda}\left(\partial_{X_{\nu}}\mathfrak{M}_{Xp}^{\lambda}-\partial_{p_{\nu}}\mathfrak{M}_{XX}^{\lambda}\right)\right]+O(\upsilon\hslash^{2}),$
(118)

where

$\mathfrak{M}_{Xp}^{\lambda}=\frac{1}{2}\epsilon^{\mu\nu\lambda}\mathsf{M}_{X_{\mu}p_{\nu}}\underline{\sf
F}_{d},\qquad\mathfrak{M}_{XX}^{\lambda}=\frac{1}{2}\epsilon^{\mu\nu\lambda}\mathsf{M}_{X_{\mu}X_{\nu}}\underline{\sf
F}_{d}.$ (119)

##  Appendix F The case of a separable position and momentum dependence

Here we specialize our theory to the case where the dependence on $X$ and $p$
of ${\sf H}$ is separable, i.e. we can write ${\sf H}(X,p)=g(X){\sf H}^{\rm
h}(p)$, where $g$ is a c-number real function, and the superscript ${\rm h}$
stands for “homogeneous,” and we additionally require
${\sf\Gamma}(X,p)={\sf\Gamma}(p)={\sf\Gamma}^{\rm h}(p)$. In other words, all
the spatial dependence is encoded in the function $g(X)$.

###  F.1 Solving the inhomogeneous problem

In this section, we show how the solution to the inhomogeneous problem Eqs.
(31,32) can be deduced from that of the homogeneous one, Eqs. (62a,62b) in the
main text.

This can be worked out perturbatively in $\hslash$. To zeroth order, comparing
Eq. (62b) to Eq. (32) yields immediately

$\displaystyle\mathsf{S}_{0}$ $\displaystyle=\mathsf{S}^{\rm h}/\sqrt{g},$
(120)

whence

$\displaystyle\mathsf{\Lambda}_{\alpha}$
$\displaystyle=-\frac{\partial_{\alpha}g}{2g}+(\mathsf{S}^{\rm
h})^{-1}\partial_{\alpha}\mathsf{S}^{\rm h}+O(\hslash).$ (121)

Explicitly in terms of position and momentum components, and up to corrections
$O(\hslash)$, this reads $\mathsf{\Lambda}_{p_{\mu}}=\mathsf{\Lambda}^{\rm
h}_{p_{\mu}}$ so that $\mathsf{A}_{p_{\mu}}=\mathsf{A}^{\rm h}_{p_{\mu}}$, and
$\mathsf{\Lambda}_{X_{\mu}}=-\frac{1}{2}\psi_{\mu}$ and in turn
$\mathsf{A}_{X_{\mu}}=0$. Here we defined $\psi_{\mu}=\partial_{X_{\mu}}g/g$
which is sometimes called a “gravitational field” Luttinger (1964); Cooper
_et al._ (1997); Qin  _et al._ (2012); Matsumoto  _et al._ (2014).

The decomposition Eq. (B.2), ${\sf K}={\sf K}_{0}+\hslash{\sf
K}_{1}+O(\hslash^{2})$, becomes here

$\displaystyle\mathsf{K}_{0}$ $\displaystyle=g\mathsf{K}^{\rm h},$
$\displaystyle\mathsf{K}_{1}$
$\displaystyle=\hbar\,\frac{i}{2}\partial_{X_{\mu}}g\partial_{p_{\mu}}\mathsf{\Gamma}^{\rm
h}\mathsf{H}^{\rm h}.$ (122)

Thus $\mathsf{K}_{0,d}=g\mathsf{K}_{d}^{\rm h}$. Then one just need to plug
Eqs. (120, F.1) into Eq. (99). Because ${\sf S}_{0}^{-1}{\sf K}_{1}{\sf
S}_{0}$ is purely anti-hermitian, its real diagonal part vanishes, so the
first $O(\hslash)$ term in Eq. (99), ${\rm Re}\left[\left({\sf S}_{0}^{-1}{\sf
K}_{1}{\sf S}_{0}\right)^{({\rm d})}\right]$, is zero. Besides, because of the
$\alpha\leftrightarrow\beta$ antisymmetry, for the last term it is sufficient
to look at $\alpha=X_{\mu},\beta=p_{\mu}$, where clearly the commutator
vanishes because $\mathsf{\Lambda}_{X_{\mu}}$ is diagonal (actually $\propto
1$), so the third $O(\hslash)$ term, $\epsilon_{\alpha\beta}{\rm
Im}\\{{\sf\Lambda}_{\beta},[{\sf\Lambda}_{\alpha},{\sf K}_{0,d}]\\}^{({\rm
d})}$, is zero as well. Consequently, the only $O(\hslash^{1})$ contribution
in Eq. (99) comes from the second $O(\hslash^{1})$ term, and

$\displaystyle\mathsf{K}_{d}$
$\displaystyle=g\left(1+\hslash\psi_{\mu}\mathsf{A}_{p_{\mu}}^{\rm
h}\right)\mathsf{K}_{d}^{\rm h}+O(\hslash^{2}).$ (123)

This solves the problem to order $O(\hslash^{1})$.

###  F.2 Distribution function

The distribution $\mathsf{F}_{d}$ close to equilibrium is a function of the
energies ${\sf K}_{d}$ only,

$\displaystyle\mathsf{F}_{d}$ $\displaystyle=f({\sf K}_{d},T)$ (124)
$\displaystyle=f(g{\sf K}_{d}^{\rm h},T)+g\psi_{\mu}\mathsf{A}_{p_{\mu}}^{\rm
h}{\sf K}_{d}^{\rm h}f^{\prime}(g{\sf K}_{d}^{\rm h},T)+O(\hslash^{2}),$

where a function of a diagonal matrix is transparently defined (cf. Eq. (44)).
It should be understood that $f$ is a real function which depends solely on
_local equilibrium_ properties, such as the local temperature Cooper  _et al._
(1997); Qin  _et al._ (2011), and $f^{\prime}$ is the derivative of $f$ with
respect to its first variable. Because $f$ does not depend on the profile of
$g(X)$, it is the same function as in the homogeneous case,
$\mathsf{F}_{d}^{\rm h}=f({\sf K}_{d}^{\rm h},T)$. It is then not difficult to
show that

$\displaystyle f(\varepsilon,T)$ $\displaystyle=\varepsilon\left(n_{\rm
B}(\varepsilon,T)+\tfrac{1}{2}\right),$ (125)

where $n_{\rm B}(\varepsilon,T)=[\exp(\varepsilon/k_{\rm B}T)-1]^{-1}$ is the
Bose function, $T$ is the local temperature, and the extra $\varepsilon$
factor comes from the normalization of eigenvectors, Eq. (62b).

###  F.3 Energy current

We are now in a position to compute the energy current, Eq. (60). First
consider Eq. (61a). Because the only non-vanishing curvature here is
$\mathsf{\Omega}_{p_{\mu}p_{\nu}}$ since ${\sf A}_{\mu}=0$, the jacobian
factor reduces to $\mathfrak{J}=1+O(\hslash^{2})$. We specialize to
$\beta=p_{\mu}$, and only $\sigma=X_{\nu}$ in the second factor can contribute
(cf. Eq. (115a)), yielding

$\displaystyle\partial_{t}\mathsf{X}_{\mu}$
$\displaystyle=\partial_{p_{\mu}}\underline{\sf
K}_{d}-\hslash{\sf\Omega}_{p_{\mu}p_{\nu}}\partial_{X_{\nu}}\underline{\sf
K}_{d}+O(\hslash^{2})$ $\displaystyle=g\left(\partial_{p_{\mu}}{\sf
K}_{d}^{\rm h}-\hslash{\sf\Omega}^{\rm h}_{p_{\mu}p_{\nu}}\psi_{\nu}{\sf
K}_{d}^{\rm h}\right)+O(\hslash^{2}),$ (126)

where the $O(\hslash)$ terms generated by going from $\underline{\sf K}_{d}$
to ${\sf K}_{d}$ canceled against each other. Finally, plugging Eq. (124) into
the third factor and expanding, it is straightforward to show that the latter
reduces to $f(g{\sf K}_{d}^{\rm h},T)$. The first contribution to the current
is thus

$\displaystyle\mathcal{J}_{X_{\mu}}^{(1)}$
$\displaystyle=\frac{\upsilon}{2}{\rm Tr}\Big{[}g\left(\partial_{p_{\mu}}{\sf
K}_{d}^{\rm h}-\hslash{\sf\Omega}^{\rm h}_{p_{\mu}p_{\nu}}\psi_{\nu}{\sf
K}_{d}^{\rm h}\right)f(g{\sf K}_{d}^{\rm h},T)\Big{]}+O(\upsilon\hslash^{2}).$
(127)

Note that the first (order $O(\hslash^{0})$) contribution is a total momentum
derivative, therefore it always vanishes after momentum integration.

Now consider Eq. (61b). Because $\sf M_{\alpha\beta}$ is antisymmetric and
$\mathsf{\Lambda}_{X_{\mu}}$ is diagonal, the only (to $O(\hslash^{0})$)
nonvanishing component of $\sf M$ is

$\displaystyle{\sf M}_{p_{\mu}p_{\nu}}$ $\displaystyle=g\frac{1}{2}{\rm
Im}\left\\{\mathsf{\Lambda}_{p_{\nu}}^{\rm
h},\left[\mathsf{\Lambda}_{p_{\mu}}^{\rm h},{\sf K}_{d}^{\rm
h}\right]\right\\}^{({\rm d})}+O(\hslash),$ (128)

whence the second contribution to the current

$\displaystyle\mathcal{J}_{X_{\mu}}^{(2)}$
$\displaystyle=\frac{\upsilon\hslash}{2}\partial_{X_{\nu}}{\rm
Tr}\big{[}\mathsf{M}_{p_{\mu}p_{\nu}}{\sf
F}_{d}\big{]}+O(\upsilon\hslash^{2})$ (129)
$\displaystyle=\frac{\upsilon\hslash}{2}{\rm
Tr}\bigg{[}g\mathsf{M}_{p_{\mu}p_{\nu}}^{\rm h}\Big{(}f(g{\sf K}_{d}^{\rm
h},T)\left(\psi_{\nu}+\partial_{\nu}T/T\right)$ $\displaystyle+g{\sf
K}_{d}^{\rm h}f^{\prime}(g{\sf K}_{d}^{\rm
h},T)\left(\psi_{\nu}-\partial_{\nu}T/T\right)\Big{)}\bigg{]}+O(\upsilon\hslash^{2}).$

Eqs. (127,129) provide the total local, momentum-resolved energy current.

###  F.4 Discussion

Here we restrict ourselves to the equilibrium case, $\partial_{X_{\nu}}T=0$.
Then $\mathcal{J}_{X_{\mu}}^{(1)}$ besides $\mathcal{J}_{X_{\mu}}^{(2)}$ is
also a total $\partial_{X_{\nu}}$ derivative, so that the net energy current
$\boldsymbol{J}_{\rm tot}$ vanishes, as it should in equilibrium. Meanwhile,
neither the _local_ momentum-integrated currents

$\displaystyle{J}_{X_{\mu}}^{(1)}$ $\displaystyle=-\upsilon\hslash
g\;\psi_{\nu}\int_{p}{\rm Tr}_{+}\Big{[}{\sf\Omega}^{\rm
h}_{p_{\mu}p_{\nu}}{\sf K}_{d}^{\rm h}\,f(g{\sf K}_{d}^{\rm
h},T)\Big{]}+O(\upsilon\hslash^{2}),$ $\displaystyle{J}_{X_{\mu}}^{(2)}$
$\displaystyle=\upsilon\hslash\;\partial_{X_{\nu}}\int_{p}{\rm
Tr}_{+}\Big{[}g\mathsf{M}^{\rm h}_{p_{\mu}p_{\nu}}\,f(g{\sf K}_{d}^{\rm
h},T)\Big{]}+O(\upsilon\hslash^{2}),$ (130)

nor their sum $J_{X_{\mu}}(X)\equiv{J}_{X_{\mu}}^{(1)}+{J}_{X_{\mu}}^{(2)}$
need vanish even at equilibrium, because the energy magnetization is not zero
a priori. Note that in Eqs.(130) the trace $\rm Tr_{+}$ runs over positive
energy eigenstates only, following the argument developed in Sec. IV.1 and
App. B.3.

There is however a physical constraint that the energy magnetization should
obey the third law of thermodynamics, so that $J_{X_{\mu}}(X)$ should at least
vanish in the limit of zero temperature. Let us show that Eqs. (127,129)
indeed satisfy this property. We will resort to the fact that
$f(\epsilon,T)\overset{T\rightarrow 0}{=}\frac{1}{2}\epsilon$ for any
$\epsilon>0$, in other words only the zero-point fluctuations contribute to
the density matrix. Taking this limit within Eq.(130), expanding the
definition Eq.(18), using cyclicity of the trace and eventually identifying
$\mathsf{\Omega}_{p_{\mu}p_{\nu}}^{\rm h}={\rm
Im}[\mathsf{\Lambda}_{p_{\nu}}^{\rm h},\mathsf{\Lambda}_{p_{\mu}}^{\rm h}]$,
it is then straightforward to show that (in equilibrium)

$\displaystyle{J}_{X_{\mu}}^{(2)}\quad$ $\displaystyle\overset{T\rightarrow
0}{=}\quad\frac{\upsilon\hslash}{2}g^{2}\psi_{\nu}\int_{p}{\rm
Tr}_{+}\left(\mathsf{K}_{d}^{\rm h}\mathsf{\Omega}_{p_{\mu}p_{\nu}}^{\rm
h}\mathsf{K}_{d}^{\rm h}\right)+O(\upsilon\hslash^{2})$
$\displaystyle\overset{T\rightarrow
0}{=}\quad-{J}_{X_{\mu}}^{(1)}\;+O(\upsilon\hslash^{2}).$ (131)

This shows, to $O(\upsilon\hslash^{1})$, the cancellation of local
magnetization currents in the zero-temperature limit (this is in fact already
true for the momentum-resolved current). More generally, the above shows that
the zero-point fluctuations (the ‘$+\frac{1}{2}$’ in the definition of $f$)
cancel exactly between ${\mathcal{J}}_{X_{\mu}}^{(1)}$ and
${\mathcal{J}}_{X_{\mu}}^{(2)}$.

##  Appendix G Explicit spatial integration for a homogeneous system with a
boundary

We now perform the $\int\text{d}x$ (recall $x\equiv X_{x}$) integration in
Eqs. (63,64). We have

$\displaystyle\kappa_{xy}^{\rm tr}$
$\displaystyle=-\frac{1}{\partial_{x}T}{\rm
Tr}_{+}\int_{p}\lim_{L\rightarrow\infty}\frac{1}{L}\int\text{d}x\,{\sf\Omega}^{\rm
h}_{p_{y}p_{x}}\partial_{x}\left({\sf K}_{d}(x)\right){\sf
K}_{d}(x)\left(n_{\rm B}({\sf K}_{d}(x),T)\big{|}_{T(x)}+\tfrac{1}{2}\right),$
(132)

where we recall ${\sf K}_{d}(x)=g(x){\sf K}_{d}^{\rm h}$. Clearly the
integrand decays over distances $O(L^{0})$ outside of $I_{x}$, therefore one
does not change the result by multiplying the former by any cutoff function
$\Lambda(x)$ strictly equal to 1 wherever the integrand is nonzero, and
decaying slowly to zero at $x\rightarrow\pm\infty$. In addition we introduce
$1=\int\text{d}\varepsilon\,\delta(\varepsilon-{\sf K}_{d}(x))$ in order to
perform a change of variables from the spatial coordinate $x$ to energy
$\varepsilon$:

$\displaystyle\kappa_{xy}^{\rm tr}$
$\displaystyle=-\frac{1}{\partial_{x}T}{\rm Tr}_{+}\int_{p}{\sf\Omega}^{\rm
h}_{p_{y}p_{x}}\lim_{L\rightarrow\infty}\frac{1}{L}\int\text{d}x\,\partial_{x}\left({\sf
K}_{d}(x)\right)\int\text{d}\varepsilon\,\delta(\varepsilon-{\sf
K}_{d}(x))\,\varepsilon\,\left(n_{\rm
B}(\varepsilon,T)\big{|}_{T(x)}+\tfrac{1}{2}\right)\Lambda(x)$
$\displaystyle=\frac{1}{\partial_{x}T}{\rm Tr}_{+}\int_{p}{\sf\Omega}^{\rm
h}_{p_{y}p_{x}}\lim_{L\rightarrow\infty}\frac{1}{L}\int\text{d}x\int\text{d}\varepsilon\,\partial_{x}\Theta(\varepsilon-{\sf
K}_{d}(x))\,\varepsilon\,\left(n_{\rm
B}(\varepsilon,T)\big{|}_{T(x)}+\tfrac{1}{2}\right)\Lambda(x)$
$\displaystyle=-\frac{1}{\partial_{x}T}{\rm Tr}_{+}\int_{p}{\sf\Omega}^{\rm
h}_{p_{y}p_{x}}\lim_{L\rightarrow\infty}\frac{1}{L}\int\text{d}x\int\text{d}\varepsilon\,\Theta(\varepsilon-{\sf
K}_{d}(x))\,\varepsilon\,\partial_{x}n_{\rm B}(\varepsilon,T)\big{|}_{T(x)}.$
(133)

In going to the second line we used the chain rule, and in going to the third
line we integrated by parts – no boundary term is generated because of the
cutoff function $\Lambda(x)$, and the extra term proportional to
$\partial_{x}\Lambda$ vanishes like $O(1/L)$. We use finally that

$\displaystyle\partial_{x}n_{\rm B}(\varepsilon,T)\big{|}_{T(x)}$
$\displaystyle=-(\varepsilon/T(x))\partial_{x}T(x)\partial_{\varepsilon}n_{\rm
B}(\varepsilon,T)\big{|}_{T(x)},$ (134)

and do not include $T(x)$ position dependence anywhere else, and assume that
$\partial_{x}T(x)$ is uniform within $I_{x}$ and zero outside, i.e.
$\partial_{x}T(x)={\rm const}\neq 0$ for $x\in I_{x}$ and $\partial_{x}T(x)=0$
for $x\notin I_{x}$. This is enough for the linear response regime up to
$O(1/L)$ corrections thanks to the decay of $g(x)$ over distances $O(L^{0})$.
Since then $g(x)=1$ in the whole support of the integrand, one can replace
${\sf K}_{d}(x)\rightarrow{\sf K}_{d}$ as well as
$\frac{1}{L}\int_{I_{x}}\text{d}x\rightarrow 1$, which yields exactly Eq.
(65), where we also used ${\sf\Omega}^{\rm h}_{p_{x}p_{y}}=-{\sf\Omega}^{\rm
h}_{p_{y}p_{x}}$.

##  Appendix H The example of inhomogeneous elasticity with a time-reversal-
breaking term

###  H.1 General

The Lagrangian density Eq. (66) is equivalent to the Hamiltonian density

$\mathcal{H}_{\rm
ph}(r)=\frac{1}{2\rho(r)}\Pi_{i}\Pi_{i}+\frac{c_{ij\mu\nu}(r)}{2}\partial_{\mu}u_{i}\partial_{\nu}u_{j},$
(135)

upon introducing the conjugate lattice momentum

$\Pi_{i}=-\frac{i}{\hbar}\frac{\delta}{\delta u_{i}}-A_{i}[u],\quad
A_{i}[u]\equiv\eta_{ij\mu\nu}(r)\partial^{2}_{\mu\nu}u_{j}.$ (136)

Thus the viscosity term introduces an effective vector potential which couples
to the lattice momentum and entails a nontrivial
$[\Pi_{i}(r),\Pi_{j}(r^{\prime})]$ commutator. For convenience, we now define
the “reduced” momentum ${\pi}_{i}(r)=\frac{\Pi_{i}(r)}{\rho(r)}$. To connect
with the general formulation of the main text, we define the 6-dimensional
vector

$\Phi(r)=\begin{bmatrix}u_{1}(r),u_{2}(r),u_{3}(r),{\pi}_{1}(r),{\pi}_{2}(r),{\pi}_{3}(r)\end{bmatrix}^{\top},$
(137)

and by comparing Eq. (135) and the commutation relations of the $u_{i}(r)$ and
$\pi_{j}(r^{\prime})$ fields to Eqs. (1,2) in the main text, one can readily
identify $\mathsf{\Gamma}$ and $\mathsf{H}$. Performing the Wigner transform
to phase space coordinates, one then obtains the dynamical matrix
$\mathsf{K}=\hbar\,\mathsf{\Gamma}\star\mathsf{H}$:

${\sf\Gamma}=\frac{i}{\hbar\rho}\begin{bmatrix}0&1\\\
-1&[\mathsf{\Gamma}_{22}]\end{bmatrix},\;{\sf
H}=\begin{bmatrix}[\mathsf{H}_{11}]&0\\\
0&\rho\end{bmatrix},\;\mathsf{K}=i\begin{bmatrix}0&1\\\
[\mathsf{K}_{21}]&[\mathsf{K}_{22}]\end{bmatrix},$ (138)

where each block is a $3\times 3$ matrix indexed by $i,j=1..3$, and

$\displaystyle[\mathsf{H}_{11}]$
$\displaystyle=\frac{1}{4}\partial^{2}_{X_{\mu}X_{\nu}}c_{\mu\nu}+\frac{1}{\hbar^{2}}c_{\mu\nu}p_{\mu}p_{\nu},$
$\displaystyle[\mathsf{\Gamma}_{22}]$
$\displaystyle=\frac{1}{2\hbar\rho}\left(\partial^{2}_{X_{\mu}X_{\nu}}\eta_{\mu\nu}+2\eta_{\mu\nu}\left(\frac{\partial_{X_{\mu}}\rho\partial_{X_{\nu}}\rho}{\rho^{2}}-\frac{\partial^{2}_{X_{\mu}X_{\nu}}\rho}{\rho}-\frac{2}{\hbar^{2}}p_{\mu}p_{\nu}\right)\right),$
$\displaystyle[\mathsf{K}_{21}]$
$\displaystyle=-\frac{1}{\rho}\left(\frac{1}{4}\partial^{2}_{X_{\mu}X_{\nu}}c_{\mu\nu}+\frac{1}{\hbar^{2}}c_{\mu\nu}p_{\mu}p_{\nu}\right)-\upsilon^{2}\frac{c_{\mu\nu}}{4\rho}\left(\frac{\partial^{2}_{X_{\mu}X_{\nu}}\rho}{\rho}-2\frac{\partial_{X_{\mu}}\rho\partial_{X_{\nu}}\rho}{\rho^{2}}\right)+\frac{i\upsilon}{\hbar\rho}c_{\mu\nu}p_{\nu}\frac{\partial_{X_{\mu}}\rho}{\rho},$
$\displaystyle[\mathsf{K}_{22}]$
$\displaystyle=\frac{1}{2\hbar\rho}\left(\partial^{2}_{X_{\mu}X_{\nu}}\eta_{\mu\nu}+2\eta_{\mu\nu}\left(\frac{\partial_{X_{\mu}}\rho\partial_{X_{\nu}}\rho}{\rho^{2}}-\Big{(}1-\frac{\upsilon^{2}}{2}\Big{)}\frac{\partial^{2}_{X_{\mu}X_{\nu}}\rho}{\rho}-\frac{2}{\hbar^{2}}p_{\mu}p_{\nu}\right)\right)+\frac{i\upsilon}{\hbar\rho}\frac{2}{\hbar}\eta_{\mu\nu}p_{\nu}\frac{\partial_{X_{\mu}}\rho}{\rho}.$
(139)

Eqs. (138–139) involve parameters
$\rho,(c_{\mu\nu})_{ij},(\eta_{\mu\nu})_{ij}$ whose dependence on position $X$
can be chosen to be any arbitrary smooth functions of position $X$, so that
Eqs. (138–139) are the general expression, in phase space, of inhomogeneous
(linear) elasticity with a Hall viscosity term.

###  H.2 Specific

Now we focus on the special case of the simple spatial dependence from App. F.
For this application we need only consider the spatially independent version
of the above, Eq. (67) in the main text. We also choose the elasticity
parameters to have the specific form Eq. (68), corresponding to the Lagrangian
density

$\displaystyle\mathcal{L}_{\rm ph}$
$\displaystyle=\frac{\rho}{2}\left(\partial_{\tau}\boldsymbol{u}\right)^{2}+\frac{1}{2}\Big{(}-c_{1}\boldsymbol{u}\cdot(\nabla^{2}\boldsymbol{u})+c_{2}(\boldsymbol{\nabla}_{\perp}\cdot\boldsymbol{u}_{\perp})^{2}$
$\displaystyle+c_{3}(\partial_{z}u^{z})^{2}+c_{4}(\boldsymbol{\nabla}_{\perp}u^{z})^{2}+c_{5}(\partial_{z}\boldsymbol{u}_{\perp})^{2}\Big{)}$
$\displaystyle+i\left(\partial_{\tau}\boldsymbol{u}\times(\eta_{1}\nabla_{\perp}^{2}+\eta_{2}\partial^{2}_{zz})\boldsymbol{u}\right)\cdot\mathbf{\hat{z}},$
(140)

where $\boldsymbol{\nabla}_{\perp}=(\partial_{x},\partial_{y},0)$ and
$\boldsymbol{u}_{\perp}=(u^{x},u^{y},0)$.

This problem can be solved analytically for arbitrary values of
$\eta_{1},\eta_{2},c_{a},a=1..5$. In the following, we expand the solution in
powers of the viscosity coefficients $\eta_{1,2}$, which we assume to be small
with respect to $c_{1,2}\sqrt{\rho}$ – this assumption is largely valid in
known relevant cases Ye  _et al._ . We furthermore take $c_{5}=0$, which has
little physical consequence but makes analytical expressions considerably
shorter. We also define $\tilde{c}_{a}=c_{a}/\rho$ for $a=1..5$ and
$\tilde{\eta}_{b}=\eta_{b}/\rho$ for $b=1,2$.

We find the energy bands

$\displaystyle{\sf K}^{\rm h}_{d}$ $\displaystyle={\rm
diag}\left(\varepsilon_{1},-\varepsilon_{1},\varepsilon_{2},-\varepsilon_{2},\varepsilon_{3},-\varepsilon_{3}\right),$
(141) $\displaystyle\varepsilon_{1}$
$\displaystyle=[\tilde{c}_{1}\boldsymbol{p}^{2}]^{\frac{1}{2}},\;\quad\varepsilon_{2}=[\tilde{c}_{1}\boldsymbol{p}^{2}+\tilde{c}_{2}p_{\perp}^{2}]^{\frac{1}{2}},$
$\displaystyle\varepsilon_{3}$
$\displaystyle=[\tilde{c}_{1}\boldsymbol{p}^{2}+\tilde{c}_{3}p_{z}^{2}+\tilde{c}_{4}p_{\perp}^{2}]^{\frac{1}{2}},$

where $p_{\perp}^{2}=p_{x}^{2}+p_{y}^{2}$, and the “normalized” eigenvectors
(normalized according to Eq.(32)) $\Psi_{i},\Psi_{i}^{*}$ such that

$\displaystyle\mathsf{S}^{\rm h}$
$\displaystyle=\begin{bmatrix}\Psi_{1}\;\vline\;\Psi_{1}^{*}\;\vline\;\Psi_{2}\;\vline\;\Psi_{2}^{*}\;\vline\;\Psi_{3}\;\vline\;\Psi_{3}^{*}\end{bmatrix}.$
(142)

In writing the explicit expressions for the eigenvectors, we use the shorthand
notations $\zeta_{\boldsymbol{p}}^{(\rm
I)}\equiv\tilde{\eta}_{1}p_{\perp}^{2}+\tilde{\eta}_{2}\boldsymbol{p}^{2}$,
$\zeta_{\boldsymbol{p}}^{(\rm II)}\equiv
p_{\perp}^{2}(\tilde{\eta}_{1}p_{\perp}^{2}-\tilde{\eta}_{2}p_{z}^{2})$ and
$\zeta_{\boldsymbol{p}}^{(\rm
III)}\equiv\tilde{\eta}_{1}p_{\perp}^{4}-2\tilde{\eta}_{2}p_{z}^{4}-\tilde{\eta}_{2}p_{z}^{2}p_{\perp}^{2}$,
so that a choice of properly normalized eigenvectors can be

$\displaystyle\Psi_{1}$
$\displaystyle=\frac{1}{\sqrt{2\rho}}\left[1+(p_{y}/p_{x})^{2}\right]^{-\frac{1}{2}}\left(\frac{i}{\varepsilon_{1}}\frac{p_{y}}{p_{x}}+\frac{2\zeta_{\boldsymbol{p}}^{(\rm
I)}}{\tilde{c}_{2}p_{x}^{2}},-\frac{i}{\varepsilon_{1}},0,\frac{p_{y}}{p_{x}}-i\varepsilon_{1}\frac{2\zeta_{\boldsymbol{p}}^{(\rm
I)}}{\tilde{c}_{2}p_{x}^{2}},-1,0\right)^{\top}\quad+O(\eta^{2}),$
$\displaystyle\Psi_{2}$
$\displaystyle=\frac{1}{\sqrt{2\rho}}\left[1+(p_{x}/p_{y})^{2}\right]^{-\frac{1}{2}}\left(\frac{i}{\varepsilon_{2}}\frac{p_{x}}{p_{y}}+\frac{2\zeta_{\boldsymbol{p}}^{(\rm
I)}}{\tilde{c}_{2}p_{y}^{2}},\frac{i}{\varepsilon_{2}},0,\frac{p_{x}}{p_{y}}-i\varepsilon_{2}\frac{2\zeta_{\boldsymbol{p}}^{(\rm
I)}}{\tilde{c}_{2}p_{y}^{2}},1,0\right)^{\top}\quad+O(\eta^{2}),$ (143)

and
$\Psi_{3}=\frac{1}{\sqrt{2\rho}}\left(0,0,\frac{i}{\varepsilon_{3}},0,0,1\right)^{\top}$.

Here we provide $\Psi_{1}$ and $\Psi_{2}$ to linear order in $\eta_{1,2}$
only. This is sufficient to compute the leading contribution to the
magnetization ${\sf M}_{p_{x}p_{y}}$ and Berry curvature
$\mathsf{\Omega}_{p_{x}p_{y}}$, which are odd under time reversal, and thus
odd in powers of $\eta_{1,2}$. In particular $\Psi_{3}$ does not depend on
$\eta_{1,2}$, therefore the magnetization and Berry curvature are zero in the
third band. In the other two bands we find, explicitly,

$\displaystyle{\sf M}_{p_{x}p_{y}}^{(1)}$
$\displaystyle=\frac{1}{\varepsilon_{1}^{2}p_{\perp}^{2}}\tilde{c}_{1}\zeta_{\boldsymbol{p}}^{(\rm
III)}\quad+O(\eta^{3}),$ $\displaystyle{\sf M}_{p_{x}p_{y}}^{(2)}$
$\displaystyle=\frac{1}{\varepsilon_{2}^{2}p_{\perp}^{2}}\left(\tilde{c}_{2}\zeta_{\boldsymbol{p}}^{(\rm
II)}+\tilde{c}_{1}\zeta_{\boldsymbol{p}}^{(\rm III)}\right)\quad+O(\eta^{3}),$
(144) $\displaystyle{\mathsf{\Omega}}_{p_{x}p_{y}}^{(1)}$
$\displaystyle=\frac{1}{\tilde{c}_{2}p_{\perp}^{4}}\frac{\tilde{c}_{1}}{\varepsilon_{1}^{3}}\left(\tilde{c}_{2}p_{\perp}^{2}\zeta_{\boldsymbol{p}}^{(\rm
II)}+2\tilde{c}_{2}\eta_{1}p_{\perp}^{4}p_{z}^{2}+4\tilde{c}_{1}\boldsymbol{p}^{2}\zeta_{\boldsymbol{p}}^{(\rm
III)}\right)\quad+O(\eta^{3}),$
$\displaystyle{\mathsf{\Omega}}_{p_{x}p_{y}}^{(2)}$
$\displaystyle=-\frac{1}{\tilde{c}_{2}p_{\perp}^{4}}\frac{1}{\varepsilon_{2}^{3}}\left(3\tilde{c}_{2}^{2}p_{\perp}^{2}\zeta_{\boldsymbol{p}}^{(\rm
II)}+4\tilde{c}_{1}^{2}\boldsymbol{p}^{2}\zeta_{\boldsymbol{p}}^{(\rm
III)}+\tilde{c}_{1}\tilde{c}_{2}p_{\perp}^{2}\left[2p_{z}^{2}(\tilde{\eta}_{1}p_{\perp}^{2}+\tilde{\eta}_{2}p_{z}^{2})+7\zeta_{\boldsymbol{p}}^{(\rm
III)}\right]\right)\quad+O(\eta^{3}).$ (145)

We make a final remark about momentum integration. The integrand in Eqs.
(127,129) does not vanish at large energies due to the zero-point contribution
to the boson density (see details in Appendix F.4 and discussion in Sec.
III.6.2). While this does not have physical consequences for the total current
$\mathcal{J}_{X_{y}}=\mathcal{J}_{X_{y}}^{(1)}+\mathcal{J}_{X_{y}}^{(2)}$, it
makes both integrals diverge independently in a continuum model. Thus, when
evaluating momentum integrals we restore a conventional Brillouin zone ${\rm
BZ}\equiv[-\pi/\xi,\pi/\xi]^{3}$.