# Supplemental Materials for
”Nonreciprocal electron hydrodynamics under magnetic fields:
applications to nonreciprocal surface magnetoplasmons”

Ryotaro Sano<EMAIL_ADDRESS>Riki Toshio  Norio Kawakami
Department of Physics, Kyoto University, Kyoto 606-8502, Japan

##  I Formulation

In this section, we outline how to derive the magnetohydrodynamic equations
for noncentrosymmetric metals where TRS is broken by an external magnetic
field. We start from the Boltzmann equation,

$\partial_{t}f^{\alpha}+\dot{{\bf\it r}}^{\alpha}\dotproduct\gradient_{{\bf\it
r}}f^{\alpha}+\dot{{\bf\it p}}^{\alpha}\dotproduct\gradient_{{\bf\it
p}}f^{\alpha}=\mathcal{C}[f^{\alpha}],$ (S1)

where $\mathcal{C}$ is the collision integral due to the scattering processes,
given by the sum of the momentum-conserving collision integral
$\mathcal{C}^{\mathrm{mc}}$ and the momentum-relaxing collision integral
$\mathcal{C}^{\mathrm{mr}}$:
$C[f^{\alpha}]=C^{\mathrm{mc}}[f^{\alpha}]+C^{\mathrm{mr}}[f^{\alpha}]$. The
semi-classical equations of motion read [1]

$\displaystyle\dot{{\bf\it r}}^{\alpha}=\tilde{{\bf\it
v}}^{\alpha}-\dot{{\bf\it k}}^{\alpha}\times{\bf\it\Omega}^{\alpha},$ (S2)
$\displaystyle\dot{{\bf\it p}}^{\alpha}=-e{\bf\it E}-e\dot{{\bf\it
r}}^{\alpha}\times{\bf\it B},$

where the electric field ${\bf\it E}({\bf\it r},t)$ depends on position and
time, while the magnetic field ${\bf\it B}(t)$ is assumed to depend only on
time for simplicity. A more detailed explanation of other physical variables
is given in the main text.

One can rewrite Eq.(S2) in a more convenient manner as

$\displaystyle D^{\alpha}\dot{{\bf\it r}}^{\alpha}$
$\displaystyle=\tilde{{\bf\it v}}^{\alpha}+\frac{e}{\hbar}{\bf\it
E}\times{\bf\it\Omega}^{\alpha}+\frac{e}{\hbar}(\tilde{{\bf\it
v}}^{\alpha}\dotproduct{\bf\it\Omega}^{\alpha}){\bf\it B},$ (S3)
$\displaystyle D^{\alpha}\dot{{\bf\it p}}^{\alpha}$ $\displaystyle=-e{\bf\it
E}-e\tilde{{\bf\it v}}^{\alpha}\times{\bf\it B}-\frac{e^{2}}{\hbar}({\bf\it
E}\dotproduct{\bf\it B}){\bf\it\Omega}^{\alpha}.$

Following the standard approach [2, 3], the continuity equations for the
particle density and the momentum are obtained in the relaxation-time
approximation for momentum-relaxing scattering processes as follows:

$\displaystyle\partial_{t}n_{{\bf\it B}}+\divergence{\bf\it J}^{n}_{{\bf\it
B}}=0,$ (S4) $\displaystyle\partial_{t}P_{{\bf\it
B},i}+\partial_{j}\Pi_{{\bf\it B},ij}=F_{i}+\Gamma_{i},$ (S5)

where $n_{{\bf\it B}}$ and ${\bf\it P}_{{\bf\it B}}$ are the particle density
and the momentum of electrons under magnetic fields respectively, ${\bf\it F}$
is the driving force due to external fields and ${\bf\it\Gamma}$ is the
momentum-relaxing force. In Eqs.(S4) and (S5), ${\bf\it J}^{n}_{{\bf\it B}}$
and $\Pi_{{\bf\it B},ij}$ are the fluxes of particle density and the momentum.
In the following, we will show the detailed derivation of Eqs.(S4) and (S5).

##  II Particle conservation law

First, multiplying the Boltzmann equation (S1) with phase-space density of
states $D^{\alpha}({\bf\it p})$ and integrating over the momentum space, we
obtain the following equation:

$\displaystyle\sum_{\alpha}\int[\differential{{\bf\it
p}}]D^{\alpha}\partial_{t}f^{\alpha}+\sum_{\alpha}\int[\differential{{\bf\it
p}}]D^{\alpha}\dot{{\bf\it r}}^{\alpha}\dotproduct\gradient_{{\bf\it
r}}f^{\alpha}+\sum_{\alpha}\int[\differential{{\bf\it
p}}]D^{\alpha}\dot{{\bf\it p}}^{\alpha}\dotproduct\gradient_{{\bf\it
p}}f^{\alpha}=\sum_{\alpha}\int[\differential{{\bf\it
p}}]D^{\alpha}\mathcal{C}[f^{\alpha}].$ (S6)

In the following, we calculate each term respectively. Integrating by parts
the first term on the left-hand side, we obtain the expression of the particle
density $n_{{\bf\it B}}$,

$\displaystyle\sum_{\alpha}\int[\differential{{\bf\it
p}}]D^{\alpha}\partial_{t}f^{\alpha}=\partial_{t}\sum_{\alpha}\int[\differential{{\bf\it
p}}]D^{\alpha}f^{\alpha}-\sum_{\alpha}\int[\differential{{\bf\it
p}}](\partial_{t}D^{\alpha})f^{\alpha}=\partial_{t}n_{{\bf\it
B}}-\frac{e}{\hbar}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\partialderivative{{\bf\it
B}}{t}\dotproduct{\bf\it\Omega}^{\alpha}f^{\alpha}.$ (S7)

Here, $n_{{\bf\it B}}\equiv\sum_{\alpha}\int[\differential{{\bf\it
p}}]D^{\alpha}f^{\alpha}$ is the particle density under magnetic fields.
Performing similar procedures for the second and third terms on the left-hand
side of Eq.(S6), we obtain

$\displaystyle\sum_{\alpha}\int[\differential{{\bf\it
p}}]D^{\alpha}\dot{{\bf\it r}}^{\alpha}\cdot\gradient_{{\bf\it r}}f^{\alpha}$
$\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it
p}}]\left(\tilde{v}_{i}^{\alpha}+\frac{e}{\hbar}\epsilon_{ijk}E_{j}\Omega_{k}^{\alpha}+\frac{e}{\hbar}(\tilde{v}_{j}\Omega_{j}^{\alpha})B_{i}\right)\partial_{i}f^{\alpha}$
$\displaystyle=\partial_{i}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\left(\tilde{v}_{i}^{\alpha}+\frac{e}{\hbar}\epsilon_{ijk}E_{j}\Omega_{k}^{\alpha}+\frac{e}{\hbar}(\tilde{v}_{j}\Omega_{j}^{\alpha})B_{i}\right)f^{\alpha}-\frac{e}{\hbar}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\epsilon_{ijk}(\partial_{i}E_{j})\Omega_{k}^{\alpha}f^{\alpha}$
$\displaystyle=\divergence{\bf\it
J}^{N}+\frac{e}{\hbar}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\partialderivative{{\bf\it
B}}{t}\dotproduct{\bf\it\Omega}^{\alpha}f^{\alpha},$ (S8)

and

$\displaystyle\sum_{\alpha}\int[\differential{{\bf\it
p}}]D^{\alpha}\dot{{\bf\it p}}^{\alpha}\dotproduct\gradient_{{\bf\it
p}}f^{\alpha}$ $\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it
p}}]\left(-eE_{i}-e\epsilon_{ijk}\tilde{v}^{\alpha}_{j}B_{k}-\frac{e^{2}}{\hbar}(E_{j}B_{j})\Omega^{\alpha}_{i}\right)\partial_{p_{i}}f^{\alpha}$
$\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it
p}}]\partial_{p_{i}}\left[\left(-eE_{i}-e\epsilon_{ijk}\tilde{v}^{\alpha}_{j}B_{k}-\frac{e^{2}}{\hbar}(E_{j}B_{j})\Omega^{\alpha}_{i}\right)f^{\alpha}\right]$
$\displaystyle\quad-\sum_{\alpha}\int[\differential{{\bf\it
p}}]\left(-e\epsilon_{ijk}(\partial_{p_{i}}\tilde{v}^{\alpha}_{j})B_{k}-\frac{e^{2}}{\hbar}(E_{j}B_{j})(\partial_{p_{i}}\Omega^{\alpha}_{i})\right)f^{\alpha}$
$\displaystyle=0.$ (S9)

Here, we have used the Maxwell equations $\curl{\bf\it E}=-\partial_{t}{\bf\it
B}$ and $\divergence{\bf\it B}=0$ in the last step of Eq.(S8) and the identity
$\epsilon_{ijk}\partial_{p_{i}}\tilde{v}_{j}^{\alpha}=\epsilon_{ijk}\partial_{p_{i}}\partial_{p_{j}}\tilde{\epsilon}^{\alpha}=0$
and
$\partial_{p_{k}}\Omega^{\alpha}_{i}=\epsilon_{ijk}\partial_{p_{i}}\partial_{p_{j}}A^{\alpha}_{k}=0$
in the last step of Eq.(S9). Here and hereafter, the Einstein summation
convention is implied for repeated indices.

For the physical reasons that the total number of particles does not change
due to scatterings, both momentum-conserving one and momentum-relaxing one, we
conclude that the right-hand side of Eq.(S6) should vanish:

$\int[\differential{{\bf\it
p}}]D^{\alpha}\mathcal{C}^{\mathrm{mc}}[f^{\alpha}]=0,\quad\int[\differential{{\bf\it
p}}]D^{\alpha}\mathcal{C}^{\mathrm{mr}}[f^{\alpha}]=0,\quad\mathrm{i.e.}\quad\sum_{\alpha}\int[\differential{{\bf\it
p}}]D^{\alpha}\mathcal{C}[f^{\alpha}]=0.$ (S10)

Summarizing the above equations (S6)-(S10), we obtain the particle
conservation law

$\partial_{t}n_{{\bf\it B}}+\divergence{\bf\it J}^{n}_{{\bf\it B}}=0,$ (S11)

owing to the cancelation of $\frac{e}{\hbar}\int[\differential{{\bf\it
p}}]\partialderivative{{\bf\it B}}{t}\dotproduct{\bf\it\Omega}^{\alpha}f$ in
Eqs.(S7) and (S8).

##  III Momentum conservation law

The momentum conserving law is obtained by multiplying the Boltzmann equation
(S1) with $p_{i}D^{\alpha}({\bf\it p})$ and integrating over the momentum
space:

$\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}D^{\alpha}\partial_{t}f^{\alpha}+\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}D^{\alpha}\dot{{\bf\it r}}^{\alpha}\dotproduct\gradient_{{\bf\it
r}}f^{\alpha}+\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}D^{\alpha}\dot{{\bf\it p}}^{\alpha}\dotproduct\gradient_{{\bf\it
p}}f^{\alpha}=\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}D^{\alpha}\mathcal{C}[f^{\alpha}].$ (S12)

Each term is calculated by taking similar procudures with particle
conservation law:

$\displaystyle\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}D^{\alpha}\partial_{t}f^{\alpha}$
$\displaystyle=\partial_{t}\sum_{\alpha}\int[dp]p_{i}D^{\alpha}f^{\alpha}-\sum_{\alpha}\int[dp]p_{i}(\partial_{t}D^{\alpha})f^{\alpha}$
$\displaystyle=\partial_{t}P_{{\bf\it
B},i}-\frac{e}{\hbar}\sum_{\alpha}\int[dp]p_{i}\partialderivative{{\bf\it
B}}{t}\dotproduct{\bf\it\Omega}^{\alpha}f^{\alpha},$ (S13)
$\displaystyle\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}D^{\alpha}\dot{{\bf\it r}}\dotproduct\gradient_{{\bf\it
r}}f^{\alpha}$ $\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}\left[\tilde{v}_{j}^{\alpha}+\frac{e}{\hbar}\epsilon_{jkl}E_{k}\Omega_{l}^{\alpha}+\frac{e}{\hbar}(\tilde{v}_{k}^{\alpha}\Omega_{k}^{\alpha})B_{j}\right]\partial_{j}f^{\alpha}$
$\displaystyle=\partial_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}\left[\tilde{v}_{j}^{\alpha}+\frac{e}{\hbar}\epsilon_{jkl}E_{k}\Omega_{l}^{\alpha}+\frac{e}{\hbar}(\tilde{v}_{k}^{\alpha}\Omega_{k}^{\alpha})B_{j}\right]f^{\alpha}-\frac{e}{\hbar}\epsilon_{jkl}(\partial_{j}E_{k})\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}\Omega_{l}^{\alpha}f^{\alpha}$
$\displaystyle=\partial_{j}\Pi_{{\bf\it
B},ij}+\frac{e}{\hbar}\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}\partialderivative{{\bf\it
B}}{t}\dotproduct{\bf\it\Omega}^{\alpha}f^{\alpha},$ (S14)
$\displaystyle\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}D^{\alpha}\dot{{\bf\it p}}^{\alpha}\dotproduct\gradient_{{\bf\it
p}}f^{\alpha}$ $\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}\left(-eE_{j}-e\epsilon_{jkl}\tilde{v}^{\alpha}_{k}B_{l}-\frac{e^{2}}{\hbar}(E_{k}B_{k})\Omega^{\alpha}_{j}\right)\partial_{p_{j}}f^{\alpha}$
$\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it
p}}]\partial_{p_{j}}\left[p_{i}\left\\{-eE_{j}-e\epsilon_{jkl}\tilde{v}^{\alpha}_{k}B_{l}-\frac{e^{2}}{\hbar}(E_{k}B_{k})\Omega^{\alpha}_{j}\right\\}f^{\alpha}\right]$
$\displaystyle\qquad-\sum_{\alpha}\int[\differential{{\bf\it
p}}](\partial_{p_{j}}p_{i})\left[-eE_{j}-e\epsilon_{jkl}\tilde{v}^{\alpha}_{k}B_{l}-\frac{e^{2}}{\hbar}(E_{k}B_{k})\Omega^{\alpha}_{j}\right]f^{\alpha}$
$\displaystyle\qquad-\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}\partial_{p_{j}}\left[-eE_{j}-e\epsilon_{jkl}\tilde{v}^{\alpha}_{k}B_{l}-\frac{e^{2}}{\hbar}(E_{k}B_{k})\Omega^{\alpha}_{j}\right]f^{\alpha}$
$\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it
p}}]\left[eE_{i}+e\epsilon_{ikl}\tilde{v}_{k}^{\alpha}B_{l}+\frac{e^{2}}{\hbar}(E_{k}B_{k})\Omega^{\alpha}_{i}\right]f^{\alpha}.$
(S15)

Here, we have defined the total momentum $P_{{\bf\it B},i}$ and the momentum
flux $\Pi_{{\bf\it B},ij}$ as follows:

$\displaystyle P_{{\bf\it B},i}$
$\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}D^{\alpha}f^{\alpha},$ (S16) $\displaystyle\Pi_{{\bf\it B},ij}$
$\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}\left[\tilde{v}_{j}^{\alpha}+\frac{e}{\hbar}\epsilon_{jkl}E_{k}\Omega_{l}^{\alpha}+\frac{e}{\hbar}(\tilde{v}_{k}^{\alpha}\Omega_{k}^{\alpha})B_{j}\right]f^{\alpha}.$
(S17)

After all, they satisfy the momentum conservation law,

$\partial_{t}P_{{\bf\it B},i}+\partial_{j}\Pi_{{\bf\it
B},ij}=F_{i}+\Gamma_{i}$ (S18)

with the driving force $F_{i}$ and the momentum-relaxing force $\Gamma_{i}$
defined as follows respectively:

$\displaystyle F_{i}$ $\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it
p}}]\left[eE_{i}+e\epsilon_{ikl}\tilde{v}_{k}^{\alpha}B_{l}+\frac{e^{2}}{\hbar}(E_{k}B_{k})\Omega^{\alpha}_{i}\right]f^{\alpha},$
(S19) $\displaystyle\Gamma_{i}$
$\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}D^{\alpha}\mathcal{C}[f^{\alpha}]=\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}D^{\alpha}\mathcal{C}^{\mathrm{mr}}[f^{\alpha}].$ (S20)

Here, we have used the fact that the integration of the product of the
momentum-conserving collision integral $\mathcal{C}^{\mathrm{mc}}[f^{\alpha}]$
and a fuction of momentum should vanish: $\int[\differential{{\bf\it
p}}]p_{i}D^{\alpha}({\bf\it p})\mathcal{C}^{\mathrm{mc}}[f^{\alpha}]=0$.

##  IV Hydrodynamic variables

In order to construct the effective hydrodynamic theory, which is correct up
to the second order in electric fields and the first order in magnetic fields,
we introduce several concepts. In the hydrodynamic regime, the system reaches
the local equilibrium via normal electron-electron scatterings which conserve
the total momentum of electron systems. For this reason, we can assume that
the distribution functions are described with the Lagrange multiplier ${\bf\it
u}$ as follows:

$f_{\mathrm{local}}^{\alpha}({\bf\it r},{\bf\it
p},t)=\frac{1}{e^{\beta(\tilde{\epsilon}^{\alpha}({\bf\it p})-{\bf\it
u}\dotproduct{\bf\it p}-\mu)}+1},$ (S21)

which is referred to as the local equilibrium distribution function. Note that
the inverse temperature $\beta$, as well as the chemical potential $\mu$ and
the drift velocity ${\bf\it u}$, generally depend on position and time.

From now on, we assume that the magnetic field is sufficiently weak:
$f_{\mathrm{local}}^{\alpha}=f_{{\bf\it u}}^{\alpha}({\bf\it p})-({\bf\it
B}\dotproduct{\bf\it m}^{\alpha}({\bf\it p}))\partial_{\epsilon}f_{{\bf\it
u}}^{\alpha}({\bf\it p})+O({\bf\it B}^{2})$. Here we have introduced the
notation

$f_{{\bf\it u}}^{\alpha}({\bf\it
p})=\frac{1}{e^{\beta(\epsilon_{0}^{\alpha}({\bf\it p})-{\bf\it
u}\dotproduct{\bf\it p}-\mu)}+1}.$ (S22)

We carry out the calculation under an assumption that the band energy has an
isotropic parabolic dispersion with the same effective mass $m$ around some
valleys: $\epsilon^{\alpha}_{0}=p^{2}/2m$, where ${\bf\it p}$ is defined as a
deviation from the valley. For example, when considering graphene with
inversion breaking or ML-TMDs, the centers of valleys correspond to $K$ and
$K^{\prime}$ points in the Brillouin zone and the above condition is
satisfied.

In the following, we assume that the underlying effective theory is invariant
under Galilean transformation:

$f_{{\bf\it u}}^{\alpha}({\bf\it p}+m{\bf\it u})=f_{0}(\epsilon_{0})$ (S23)

for a free-like dispersion $\epsilon_{0}({\bf\it p})={\bf\it p}^{2}/2m$ [4].
Here, $f_{0}(\epsilon)\equiv[e^{\beta(\epsilon-\mu_{0})}+1]^{-1}$ is the
conventional Fermi-Dirac distribution function. Eq.(S23) is obtained by
transforming the energy and the chemical potential from the frame of reference
$\mathcal{S}$ in which the fluid moves with velocity ${\bf\it u}$ to the frame
$\mathcal{S}_{0}$ in which it is at rest: $\mu=\mu_{0}-m{\bf\it u}^{2}/2$.

###  IV.1 I. particle conservation law

We are now ready to express two conservation laws in terms of hydrodynamic
variables. First, the particle density $n_{{\bf\it B}}$ can be expressed as
follows:

$\displaystyle n_{{\bf\it B}}$
$\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it p}}]D^{\alpha}({\bf\it
p})f^{\alpha}_{\mathrm{local}}$
$\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it p}}]D^{\alpha}({\bf\it
p})\biggl{[}f_{{\bf\it u}}^{\alpha}({\bf\it p})-{\bf\it B}\dotproduct{\bf\it
m}^{\alpha}({\bf\it p})\partial_{\epsilon}f_{{\bf\it u}}^{\alpha}({\bf\it
p})\biggr{]}+O({\bf\it B}^{2})$
$\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it p}}]D^{\alpha}({\bf\it
p}+m{\bf\it u})\biggl{[}f_{{\bf\it u}}^{\alpha}({\bf\it p}+m{\bf\it
u})-{\bf\it B}\dotproduct{\bf\it m}^{\alpha}({\bf\it p}+m{\bf\it
u})\partial_{\epsilon}f_{{\bf\it u}}^{\alpha}({\bf\it p}+m{\bf\it
u})\biggr{]}+O({\bf\it B}^{2})$
$\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it
p}}]\left(1+\frac{e}{\hbar}{\bf\it B}\cdot{\bf\it\Omega}^{\alpha}({\bf\it
p}+m{\bf\it u})\right)f^{\alpha}_{0}(\epsilon_{0})-{\bf\it
B}\dotproduct\sum_{\alpha}\int[\differential{{\bf\it p}}]{\bf\it
m}^{\alpha}({\bf\it p}+m{\bf\it
u})\partial_{\epsilon}f^{\alpha}_{0}(\epsilon_{0})+O({\bf\it B}^{2})$
$\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it
p}}]f_{0}(\epsilon_{0})+\frac{e}{\hbar}B_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\frac{\partial\Omega^{\alpha}_{j}}{\partial
p_{i}}mu_{i}f_{0}(\epsilon_{0})-B_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\frac{\partial m^{\alpha}_{j}}{\partial
p_{i}}mu_{i}\partial_{\epsilon}f_{0}(\epsilon_{0})+O({\bf\it B}^{2},{\bf\it
B}{\bf\it u}^{2})$
$\displaystyle=n_{0}+\frac{me}{\hbar}u_{i}B_{j}(D_{ij}+N_{ij})+O({\bf\it
B}^{2},{\bf\it B}{\bf\it u}^{2}).$ (S24)

Here, we have defined the particle number density in the absence of magnetic
field

$n_{0}\equiv\sum_{\alpha}\int[\differential{{\bf\it p}}]f_{{\bf\it u}}({\bf\it
p})=\sum_{\alpha}\int[\differential{{\bf\it p}}]f_{{\bf\it u}}({\bf\it
p}+m{\bf\it u})=\sum_{\alpha}\int[\differential{{\bf\it
p}}]f_{0}(\epsilon_{0}),$ (S25)

and one can observe that new geometrical tensors such as Berry curvature
dipole $\hat{D}$ [5] and $\hat{N}$ appear:

$D_{ij}\equiv\sum_{\alpha}\int[\differential{{\bf\it
p}}]\frac{\partial\Omega^{\alpha}_{j}}{\partial
p_{i}}f_{0}(\epsilon_{0}),\quad
N_{ij}\equiv-\frac{\hbar}{e}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\frac{\partial m_{j}^{\alpha}}{\partial
p_{i}}\partial_{\epsilon}f_{0}(\epsilon_{0}).$ (S26)

Surprisingly, Eq.(S24) shows that an interplay between the magnetic field
${\bf\it B}$ and the drift velocity ${\bf\it u}$ changes the particle number
density due to the lack of inversion symmetry.

The same calculations can be performed for the particle number flux ${\bf\it
J}^{n}_{{\bf\it B}}$ :

$\displaystyle{\bf\it J}^{n}_{{\bf\it
B}}=\sum_{\alpha}\int[\differential{{\bf\it p}}]\left[\tilde{{\bf\it
v}}^{\alpha}+\frac{e}{\hbar}{\bf\it
E}\times{\bf\it\Omega}^{\alpha}+\frac{e}{\hbar}(\tilde{{\bf\it
v}}^{\alpha}\dotproduct{\bf\it\Omega}^{\alpha}){\bf\it
B}\right]f_{\mathrm{local}}^{\alpha}.$ (S27)

Now we consider the $i$-component of Eq.(S27). Each term can be calculated as
follows.

$\displaystyle\sum_{\alpha}\int[\differential{{\bf\it
p}}]\tilde{v}^{\alpha}_{i}f_{\mathrm{local}}^{\alpha}$
$\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it
p}}]\tilde{v}^{\alpha}_{i}({\bf\it p})\biggl{[}f_{{\bf\it u}}^{\alpha}({\bf\it
p})-{\bf\it B}\dotproduct{\bf\it m}^{\alpha}({\bf\it
p})\partial_{\epsilon}f_{{\bf\it u}}^{\alpha}({\bf\it p})\biggr{]}+O({\bf\it
B}^{2})$ $\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it
p}}]\tilde{v}^{\alpha}_{i}({\bf\it p}+m{\bf\it u})\biggl{[}f_{{\bf\it
u}}^{\alpha}({\bf\it p}+m{\bf\it u})-{\bf\it B}\dotproduct{\bf\it
m}^{\alpha}({\bf\it p}+m{\bf\it u})\partial_{\epsilon}f_{{\bf\it
u}}^{\alpha}({\bf\it p}+m{\bf\it u})\biggr{]}+O({\bf\it B}^{2})$
$\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it
p}}](v^{\alpha}_{i}({\bf\it p}+m{\bf\it
u})-B_{j}\partial_{p_{i}}m^{\alpha}_{j}({\bf\it p}+m{\bf\it
u}))f_{0}(\epsilon_{0})$ $\displaystyle\quad-
B_{j}\sum_{\alpha}\int[\differential{\mathbf{p}}]v_{i}^{\alpha}({\bf\it
p}+m{\bf\it u})m^{\alpha}_{j}({\bf\it p}+m{\bf\it
u})\partial_{\epsilon}f_{0}(\epsilon_{0})+O({\bf\it B}^{2})$
$\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it
p}}]\left[\frac{p_{i}+mu_{i}}{m}-B_{j}\partial_{p_{i}}\left(m^{\alpha}_{j}({\bf\it
p})+\partialderivative{m^{\alpha}_{j}}{p_{k}}mu_{k}\right)\right]f_{0}(\epsilon_{0})$
$\displaystyle\quad-B_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\frac{p_{i}+mu_{i}}{m}\left(m_{j}^{\alpha}({\bf\it
p})+\partialderivative{m_{j}^{\alpha}}{p_{k}}mu_{k}\right)\partial_{\epsilon}f_{0}(\epsilon_{0})+O({\bf\it
B}^{2},{\bf\it B}{\bf\it u}^{2})$
$\displaystyle=u_{i}\sum_{\alpha}\int[\differential{{\bf\it
p}}]f_{0}(\epsilon_{0})-B_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\partialderivative{m_{j}^{\alpha}}{p_{i}}f_{0}(\epsilon_{0})-B_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\frac{p_{i}}{m}m_{j}^{\alpha}({\bf\it
p})\partial_{\epsilon}f_{0}(\epsilon_{0})+O({\bf\it B}^{2},{\bf\it B}{\bf\it
u}^{2})$
$\displaystyle=n_{0}u_{i}-B_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\partialderivative{m_{j}^{\alpha}}{p_{i}}f_{0}(\epsilon_{0})-B_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}]m_{j}^{\alpha}\partial_{p_{i}}f_{0}(\epsilon_{0})+O({\bf\it B}^{2},{\bf\it
B}{\bf\it u}^{2})$
$\displaystyle=n_{0}u_{i}-B_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\partialderivative{m_{j}^{\alpha}}{p_{i}}f_{0}(\epsilon_{0})+B_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\partialderivative{m_{j}^{\alpha}}{p_{i}}f_{0}(\epsilon_{0})+O({\bf\it
B}^{2},{\bf\it B}{\bf\it u}^{2})$ $\displaystyle=n_{0}u_{i}+O({\bf\it
B}^{2},{\bf\it B}{\bf\it u}^{2}).$ (S28)
$\displaystyle\sum_{\alpha}\int[\differential{{\bf\it
p}}]\epsilon_{ijk}E_{j}\Omega^{\alpha}_{k}f_{\mathrm{local}}^{\alpha}$
$\displaystyle=\epsilon_{ijk}E_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\Omega^{\alpha}_{k}({\bf\it p})[f_{{\bf\it u}}^{\alpha}({\bf\it
p})-B_{l}m^{\alpha}_{l}({\bf\it p})\partial_{\epsilon}f_{{\bf\it
u}}^{\alpha}({\bf\it p})]+O({\bf\it B}^{2})$
$\displaystyle=\epsilon_{ijk}E_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\Omega^{\alpha}_{k}({\bf\it p}+m{\bf\it u})[f_{{\bf\it
u}}^{\alpha}({\bf\it p}+m{\bf\it u})-B_{l}m_{l}^{\alpha}({\bf\it p}+m{\bf\it
u})\partial_{\epsilon}f_{{\bf\it u}}^{\alpha}({\bf\it p}+m{\bf\it
u})]+O({\bf\it B}^{2})$
$\displaystyle=\epsilon_{ijk}E_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\left(\Omega^{\alpha}_{k}({\bf\it
p})+\partialderivative{\Omega^{\alpha}_{k}}{p_{m}}mu_{m}\right)\left[f_{0}(\epsilon_{0})-B_{l}\left(m_{l}^{\alpha}({\bf\it
p})+\partialderivative{m_{l}^{\alpha}}{p_{n}}mu_{n}\right)\partial_{\epsilon}f_{0}(\epsilon_{0})\right]+O({\bf\it
B}^{2},{\bf\it B}{\bf\it u}^{2})$
$\displaystyle=\epsilon_{ijk}E_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\partialderivative{\Omega^{\alpha}_{k}}{p_{l}}mu_{l}f_{0}(\epsilon_{0})-\epsilon_{ijk}E_{j}B_{l}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\Omega^{\alpha}_{k}m^{\alpha}_{l}\partial_{\epsilon}f_{0}(\epsilon_{0})+O({\bf\it
B}^{2},{\bf\it B}{\bf\it u}^{2})$
$\displaystyle=\epsilon_{ijk}mE_{j}(D_{lk}u_{l}+Y_{lk}B_{l}).$ (S29)
$\displaystyle\sum_{\alpha}\int[\differential{{\bf\it
p}}]\tilde{v}^{\alpha}_{j}\Omega^{\alpha}_{j}B_{i}f_{\mathrm{local}}^{\alpha}$
$\displaystyle=B_{i}\sum_{\alpha}\int[\differential{{\bf\it
p}}]v^{\alpha}_{j}({\bf\it p})\Omega^{\alpha}_{j}({\bf\it p})f_{{\bf\it
u}}^{\alpha}({\bf\it p})+O({\bf\it B}^{2})$
$\displaystyle=B_{i}\sum_{\alpha}\int[\differential{{\bf\it
p}}]v^{\alpha}_{j}({\bf\it p}+m{\bf\it u})\Omega^{\alpha}_{j}({\bf\it
p}+m{\bf\it u})f_{{\bf\it u}}^{\alpha}({\bf\it p}+m{\bf\it u})+O({\bf\it
B}^{2})$ $\displaystyle=B_{i}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\frac{p_{j}+mu_{j}}{m}\left(\Omega^{\alpha}_{j}({\bf\it
p})+\partialderivative{\Omega_{j}^{\alpha}}{p_{k}}mu_{k}\right)f_{0}(\epsilon_{0})+O({\bf\it
B}^{2},{\bf\it B}{\bf\it u}^{2})$
$\displaystyle=B_{i}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\frac{p_{j}}{m}\Omega^{\alpha}_{j}f_{0}(\epsilon_{0})+O({\bf\it
B}^{2},{\bf\it B}{\bf\it u}^{2})$
$\displaystyle=B_{i}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\Omega_{j}^{\alpha}\partial_{p_{j}}f_{0}(\epsilon_{0})+O({\bf\it
B}^{2},{\bf\it B}{\bf\it u}^{2})$
$\displaystyle=-B_{i}\sum_{\alpha}\int[\differential{{\bf\it
p}}](\partial_{p_{j}}\Omega_{j}^{\alpha})f_{0}({\bf\it p})+O({\bf\it
B}^{2},{\bf\it B}{\bf\it u}^{2})$ $\displaystyle=O({\bf\it B}^{2},{\bf\it
B}{\bf\it u}^{2}).$ (S30)

In the above calculation, we have used the relation
$\partial_{p_{i}}f_{0}(\epsilon_{0}^{\alpha})=\partialderivative{\epsilon_{0}^{\alpha}}{p_{i}}\partial_{\epsilon}f_{0}(\epsilon_{0}^{\alpha})=\frac{p_{i}}{m}\partial_{\epsilon}f_{0}(\epsilon_{0}^{\alpha})$
and defined the new geometrical tensor

$Y_{lk}=-\frac{1}{m}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\Omega_{k}^{\alpha}m_{l}^{\alpha}\partial_{\epsilon}f_{0}(\epsilon_{0}).$
(S31)

In the end, we can obtain the particle number flux ${\bf\it J}^{n}_{{\bf\it
B}}$ as follows,

$\displaystyle{\bf\it J}^{n}_{{\bf\it B}}=n_{0}{\bf\it
u}+\frac{me}{\hbar}{\bf\it E}\times({}^{t}\hat{D}{\bf\it
u}+{}^{t}\hat{Y}{\bf\it B})+O({\bf\it B}^{2},{\bf\it B}{\bf\it u}^{2}).$ (S32)

Summarizing the above equations, the particle conservation law Eq. (S4) is
expressed as

$\displaystyle\partialderivative{t}\left[n_{0}+\frac{me}{\hbar}u_{i}B_{j}(D_{ij}+N_{ij})\right]+\divergence\left[n_{0}{\bf\it
u}+\frac{me}{\hbar}{\bf\it E}\times({}^{t}\hat{D}{\bf\it
u}+{}^{t}\hat{Y}{\bf\it B})\right]=0$ (S33)

up to the order of ${\bf\it B}^{2}$, ${\bf\it B}{\bf\it u}^{2}$.

###  IV.2 II. momentum conservation law

First we calculate the total momentum:

$\displaystyle P_{{\bf\it B},i}$
$\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}D^{\alpha}({\bf\it p})f^{\alpha}_{\mathrm{local}}$
$\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}D^{\alpha}({\bf\it p})[f^{\alpha}_{{\bf\it u}}({\bf\it p})-{\bf\it
B}\dotproduct{\bf\it m}^{\alpha}({\bf\it p})\partial_{\epsilon}f_{{\bf\it
u}}^{\alpha}({\bf\it p})]+O({\bf\it B}^{2})$
$\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it
p}}](p_{i}+mu_{i})D^{\alpha}({\bf\it p}+m{\bf\it u})[f^{\alpha}_{{\bf\it
u}}({\bf\it p}+m{\bf\it u})-{\bf\it B}\dotproduct{\bf\it m}^{\alpha}({\bf\it
p}+m{\bf\it u})\partial_{\epsilon}f_{{\bf\it u}}({\bf\it p}+m{\bf\it
u})]+O({\bf\it B}^{2})$ $\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it
p}}](p_{i}+mu_{i})\left(1+\frac{e}{\hbar}{\bf\it
B}\dotproduct{\bf\it\Omega}^{\alpha}({\bf\it p}+m{\bf\it
u})\right)f_{0}(\epsilon_{0})-{\bf\it
B}\dotproduct\sum_{\alpha}\int[\differential{{\bf\it p}}](p_{i}+mu_{i}){\bf\it
m}^{\alpha}({\bf\it p}+m{\bf\it
u})\partial_{\epsilon}f_{0}(\epsilon_{0})+O({\bf\it B}^{2})$
$\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it
p}}](p_{i}+mu_{i})\left[1+\frac{e}{\hbar}B_{j}\left(\Omega^{\alpha}_{j}({\bf\it
p})+\partialderivative{\Omega^{\alpha}_{j}}{p_{k}}mu_{k}\right)\right]f_{0}({\bf\it
p})$ $\displaystyle\quad-B_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}](p_{i}+mu_{i})\left(m_{j}^{\alpha}({\bf\it
p})+\partialderivative{m_{j}^{\alpha}}{p_{k}}mu_{k}\right)\partial_{\epsilon}f_{0}(\epsilon_{0})+O({\bf\it
B}^{2},{\bf\it B}{\bf\it u}^{2})$
$\displaystyle=mu_{i}\sum_{\alpha}\int[\differential{{\bf\it
p}}]f_{0}(\epsilon_{0})+\frac{e}{\hbar}B_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}\Omega_{j}^{\alpha}f_{0}(\epsilon_{0})-B_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}m_{j}^{\alpha}\partial_{\epsilon}f_{0}(\epsilon_{0})+O({\bf\it
B}^{2},{\bf\it B}{\bf\it u}^{2})$
$\displaystyle=mu_{i}n_{0}+\frac{e}{\hbar}B_{j}(C_{ij}+M_{ij}).$ (S34)

In the third term, we have used the relation for parabolic dispersion
$m\partial_{p_{i}}f_{0}(\epsilon_{0}^{\alpha})=p_{i}\partial_{\epsilon}f_{0}(\epsilon_{0}^{\alpha})$,

$\displaystyle-\int[\differential{{\bf\it
p}}]p_{i}m_{j}^{\alpha}\partial_{\epsilon}f_{0}(\epsilon_{0})=-m\int[\differential{{\bf\it
p}}]m_{j}^{\alpha}\partial_{p_{i}}f_{0}(\epsilon_{0})=m\int[\differential{{\bf\it
p}}]\partialderivative{m_{j}^{\alpha}}{p_{i}}f_{0}(\epsilon_{0}).$ (S35)

Here, two new geometrical tensors appear in the expression of the total
momentum:

$\displaystyle C_{ij}=\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}\Omega_{j}^{\alpha}f_{0}(\epsilon_{0}),\quad
M_{ij}=\frac{m\hbar}{e}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\frac{\partial m_{j}^{\alpha}}{\partial p_{i}}f_{0}(\epsilon_{0}).$ (S36)

Next, we would like to calculate the momentum flux

$\displaystyle\Pi_{{\bf\it B},ij}=\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}\left[\tilde{v}_{j}^{\alpha}+\frac{e}{\hbar}\epsilon_{jkl}E_{k}\Omega_{l}^{\alpha}+\frac{e}{\hbar}(\tilde{v}_{k}^{\alpha}\Omega_{k}^{\alpha})B_{j}\right]f^{\alpha}_{\mathrm{local}}.$
(S37)

The first term of $\Pi_{{\bf\it B},ij}$ is calculated as

$\displaystyle\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}\tilde{v}_{j}^{\alpha}({\bf\it p})f^{\alpha}_{\mathrm{local}}$
$\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}\tilde{v}_{j}^{\alpha}({\bf\it p})[f_{{\bf\it u}}^{\alpha}({\bf\it
p})-{\bf\it B}\dotproduct{\bf\it m}^{\alpha}({\bf\it
p})\partial_{\epsilon}f_{{\bf\it u}}^{\alpha}({\bf\it p})]+O({\bf\it B}^{2})$
$\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it
p}}](p_{i}+mu_{i})\tilde{v}_{j}^{\alpha}({\bf\it p}+m{\bf\it u})[f_{{\bf\it
u}}^{\alpha}({\bf\it p}+m{\bf\it u})-{\bf\it B}\dotproduct{\bf\it
m}^{\alpha}({\bf\it p}+m{\bf\it u})\partial_{\epsilon}f_{{\bf\it
u}}^{\alpha}({\bf\it p}+m{\bf\it u})]+O({\bf\it B}^{2})$
$\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it
p}}](p_{i}+mu_{i})\left[\frac{p_{j}+mu_{j}}{m}-B_{k}\partial_{p_{j}}\left(m_{k}^{\alpha}({\bf\it
p})+\partialderivative{m_{k}^{\alpha}}{p_{l}}mu_{l}\right)\right]f_{0}(\epsilon_{0})$
$\displaystyle\quad-B_{k}\sum_{\alpha}\int[\differential{{\bf\it
p}}](p_{i}+mu_{i})\frac{p_{j}+mu_{j}}{m}\left(m_{k}^{\alpha}({\bf\it
p})+\partialderivative{m_{k}^{\alpha}}{p_{l}}mu_{l}\right)\partial_{\epsilon}f_{0}(\epsilon_{0})+O({\bf\it
B}^{2},{\bf\it B}{\bf\it u}^{2})$
$\displaystyle=\frac{1}{m}\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}p_{j}f_{0}(\epsilon_{0})+mu_{i}u_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}]f_{0}(\epsilon_{0})$ $\displaystyle\quad-
B_{k}mu_{i}\int[\differential{{\bf\it
p}}]\partialderivative{m_{k}^{\alpha}}{p_{j}}f_{0}(\epsilon_{0})-B_{k}\int[\differential{{\bf\it
p}}]p_{i}\frac{\partial^{2}m_{k}^{\alpha}}{\partial p_{j}\partial
p_{l}}mu_{l}f_{0}(\epsilon_{0})$ $\displaystyle\quad-
B_{k}\sum_{\alpha}\int[\differential{{\bf\it
p}}](p_{i}u_{j}+u_{i}p_{j})m_{k}^{\alpha}\partial_{\epsilon}f_{0}(\epsilon_{0})-B_{k}\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}p_{j}\partialderivative{m^{\alpha}_{k}}{p_{l}}u_{l}\partial_{\epsilon}f_{0}(\epsilon_{0})+O({\bf\it
B}^{2},{\bf\it B}{\bf\it u}^{2}).$ (S38)

The first term of Eq. (S38) becomes [2]

$\frac{1}{m}\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}p_{j}f_{0}(\epsilon_{0})=p\delta_{ij}.$ (S39)

The fourth term of Eq.(S38) can be calculated as

$\displaystyle-B_{k}\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}\frac{\partial^{2}m_{k}^{\alpha}}{\partial p_{j}\partial
p_{l}}mu_{l}f_{0}(\epsilon_{0})$
$\displaystyle=-B_{k}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\partial_{p_{j}}\left(p_{i}\partialderivative{m_{k}^{\alpha}}{p_{l}}mu_{l}f_{0}(\epsilon_{0})\right)+B_{k}\sum_{\alpha}\int[\differential{{\bf\it
p}}](\partial_{p_{j}}p_{i})\partialderivative{m_{k}^{\alpha}}{p_{l}}mu_{l}f_{0}(\epsilon_{0})$
$\displaystyle\quad+B_{k}\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}\partialderivative{m_{k}^{\alpha}}{p_{l}}mu_{l}\partial_{p_{j}}f_{0}(\epsilon_{0})$
$\displaystyle=mB_{k}u_{l}\delta_{ij}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\partialderivative{m_{k}^{\alpha}}{p_{l}}f_{0}(\epsilon_{0})+B_{k}\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}\partialderivative{m_{k}^{\alpha}}{p_{l}}u_{l}p_{j}\partial_{\epsilon}f_{0}(\epsilon_{0})$
$\displaystyle=\frac{e}{\hbar}B_{k}u_{l}\delta_{ij}M_{lk}+B_{k}\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}p_{j}\partialderivative{m_{k}^{\alpha}}{p_{l}}u_{l}\partial_{\epsilon}f_{0}(\epsilon).$
(S40)

The sixth term of Eq.(S38) can be calculated as

$\displaystyle-B_{k}\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}u_{j}m_{k}^{\alpha}\partial_{\epsilon}f_{0}(\epsilon_{0})+(i\leftrightarrow
j)$ $\displaystyle=-mB_{k}u_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}]m_{k}^{\alpha}\partial_{p_{i}}f_{0}(\epsilon_{0})+(i\leftrightarrow j)$
$\displaystyle=mB_{k}u_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\partialderivative{m_{k}^{\alpha}}{p_{i}}f_{0}(\epsilon_{0})+(i\leftrightarrow
j)$ $\displaystyle=\frac{e}{\hbar}B_{k}u_{j}M_{ik}+(i\leftrightarrow j).$
(S41)

Finally, we obtain

$\displaystyle\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}\tilde{v}_{j}^{\alpha}f_{\mathrm{local}}^{\alpha}=mn_{0}u_{i}u_{j}+p\delta_{ij}+\frac{e}{\hbar}B_{k}(u_{l}M_{lk}\delta_{ij}+u_{j}M_{ik})+O({\bf\it
B}^{2},{\bf\it B}{\bf\it u}^{2}).$ (S42)

The rest term of $\Pi_{{\bf\it B},ij}$ can be calculated in a similar way.

$\displaystyle\epsilon_{jkl}E_{k}\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}\Omega_{l}^{\alpha}f_{\mathrm{local}}^{\alpha}$
$\displaystyle=\epsilon_{jkl}E_{k}\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}\Omega_{l}^{\alpha}({\bf\it p})[f_{{\bf\it u}}^{\alpha}({\bf\it
p})-{\bf\it B}\dotproduct{\bf\it m}^{\alpha}({\bf\it
p})\partial_{\epsilon}f_{{\bf\it u}}^{\alpha}({\bf\it p})]+O({\bf\it B}^{2})$
$\displaystyle=\epsilon_{jkl}E_{k}\sum_{\alpha}\int[\differential{{\bf\it
p}}](p_{i}+mu_{i})\Omega_{l}^{\alpha}({\bf\it p}+m{\bf\it u})[f_{{\bf\it
u}}^{\alpha}({\bf\it p}+m{\bf\it u})-{\bf\it B}\dotproduct{\bf\it
m}^{\alpha}({\bf\it p}+m{\bf\it u})\partial_{\epsilon}f_{{\bf\it
u}}^{\alpha}({\bf\it p}+m{\bf\it u})]+O({\bf\it B}^{2})$
$\displaystyle=\epsilon_{jkl}E_{k}\sum_{\alpha}\int[\differential{{\bf\it
p}}](p_{i}+mu_{i})\Omega_{l}^{\alpha}({\bf\it p}+m{\bf\it
u})f_{0}(\epsilon_{0})$
$\displaystyle\qquad-\epsilon_{jkl}E_{k}B_{n}\sum_{\alpha}\int[\differential{{\bf\it
p}}](p_{i}+mu_{i})\Omega_{l}^{\alpha}({\bf\it p}+m{\bf\it
u})m_{n}^{\alpha}({\bf\it p}+m{\bf\it
u})\partial_{\epsilon}f_{0}(\epsilon_{0})+O({\bf\it B}^{2})$
$\displaystyle=\epsilon_{jkl}E_{k}\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}\Omega_{l}^{\alpha}f_{0}(\epsilon_{0})+O({\bf\it B}^{2},{\bf\it
B}{\bf\it u}^{2},{\bf\it u}^{3})$
$\displaystyle=\epsilon_{jkl}E_{k}C_{il}+O({\bf\it B}^{2},{\bf\it B}{\bf\it
u}^{2},{\bf\it u}^{3}),$ (S43) $\displaystyle
B_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}\tilde{v}_{k}^{\alpha}\Omega_{k}^{\alpha}f_{\mathrm{local}}^{\alpha}$
$\displaystyle=B_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}v_{k}^{\alpha}({\bf\it p})\Omega_{k}^{\alpha}({\bf\it p})f_{{\bf\it
u}}^{\alpha}({\bf\it p})+O({\bf\it B}^{2})$
$\displaystyle=B_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}](p_{i}+mu_{i})v_{k}^{\alpha}({\bf\it p}+m{\bf\it
u})\Omega_{k}^{\alpha}({\bf\it p}+m{\bf\it u})f_{{\bf\it u}}^{\alpha}({\bf\it
p}+m{\bf\it u})+O({\bf\it B}^{2})$
$\displaystyle=B_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}](p_{i}+mu_{i})\frac{p_{k}+mu_{k}}{m}\Omega_{k}^{\alpha}({\bf\it
p}+m{\bf\it u})f_{0}(\epsilon_{0})+O({\bf\it B}^{2})$
$\displaystyle=B_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}p_{k}\partialderivative{\Omega^{\alpha}_{k}}{p_{l}}u_{l}f_{0}(\epsilon_{0})+B_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}](p_{i}u_{k}+u_{i}p_{k})\Omega_{k}^{\alpha}f_{0}(\epsilon_{0})+O({\bf\it
B}^{2},{\bf\it B}{\bf\it u}^{2}).$ (S44)

To obtain a simpler form, we perform partial integration of the first term on
the right-hand side of Eq. (S44),

$\displaystyle u_{l}\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}p_{k}\partialderivative{\Omega^{\alpha}_{k}}{p_{l}}f_{0}(\epsilon_{0})$
$\displaystyle=u_{l}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\partial_{p_{l}}[\Omega^{\alpha}_{k}p_{i}p_{k}f_{0}(\epsilon_{0})]-u_{l}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\Omega^{\alpha}_{k}\partial_{p_{l}}[p_{i}p_{k}f_{0}(\epsilon_{0})]$
$\displaystyle=-u_{i}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\Omega^{\alpha}_{k}p_{k}f_{0}(\epsilon_{0})-u_{k}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\Omega^{\alpha}_{k}p_{i}f_{0}(\epsilon_{0})-u_{l}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\Omega^{\alpha}_{k}p_{i}p_{k}\partial_{p_{l}}f_{0}(\epsilon_{0})$
$\displaystyle=u_{i}C_{kk}+u_{k}C_{ik}-u_{l}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\Omega^{\alpha}_{k}p_{i}p_{k}\partial_{p_{l}}f_{0}(\epsilon_{0}),$ (S45)

and Eq. (S44) reduces to

$\displaystyle B_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}\tilde{v}_{k}^{\alpha}\Omega_{k}^{\alpha}f_{\mathrm{local}}^{\alpha}$
$\displaystyle=-B_{j}u_{l}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\Omega^{\alpha}_{k}p_{i}p_{k}\partial_{p_{l}}f_{0}(\epsilon_{0})+O({\bf\it
B}^{2},{\bf\it B}{\bf\it u}^{2})$ $\displaystyle=B_{j}u_{l}O_{il}+O({\bf\it
B}^{2},{\bf\it B}{\bf\it u}^{2}).$ (S46)

Here we have defined the new geometrical tensor

$O_{il}=-\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}p_{k}\Omega^{\alpha}_{k}\partial_{p_{l}}f_{0}(\epsilon_{0}).$ (S47)

After all, we obtain the following equation for the momenum flux:

$\displaystyle\Pi_{{\bf\it B},ij}$
$\displaystyle=mn_{0}u_{i}u_{j}+p\delta_{ij}+\frac{e}{\hbar}\epsilon_{jkl}E_{k}C_{il}+\frac{e}{\hbar}B_{k}\left(u_{l}M_{lk}\delta_{ij}+u_{j}M_{ik}\right)+\frac{e}{\hbar}B_{j}u_{l}O_{il}+O({\bf\it
B}^{2},{\bf\it B}{\bf\it u}^{2}).$ (S48)

Driving force is defined as

$\displaystyle F_{i}$ $\displaystyle=-\sum_{\alpha}\int[\differential{{\bf\it
p}}]\left[eE_{i}+e\epsilon_{ikl}\tilde{v}_{k}^{\alpha}B_{l}+\frac{e^{2}}{\hbar}(E_{k}B_{k})\Omega^{\alpha}_{i}\right]f^{\alpha}_{\mathrm{local}}.$
(S49)

The $i$-component of each term can be calculated as follows:

$\displaystyle eE_{i}\sum_{\alpha}\int[\differential{{\bf\it
p}}]f_{\mathrm{local}}^{\alpha}$
$\displaystyle=eE_{i}\sum_{\alpha}\int[\differential{{\bf\it p}}][f_{{\bf\it
u}}^{\alpha}({\bf\it p})-{\bf\it B}\dotproduct{\bf\it m}^{\alpha}({\bf\it
p})\partial_{\epsilon}f_{{\bf\it u}}^{\alpha}({\bf\it p})]+O({\bf\it B}^{2})$
$\displaystyle=eE_{i}\sum_{\alpha}\int[\differential{{\bf\it p}}][f_{{\bf\it
u}}^{\alpha}({\bf\it p}+m{\bf\it u})-{\bf\it B}\dotproduct{\bf\it
m}^{\alpha}({\bf\it p}+m{\bf\it u})\partial_{\epsilon}f_{{\bf\it
u}}^{\alpha}({\bf\it p}+m{\bf\it u})]+O({\bf\it B}^{2})$
$\displaystyle=eE_{i}\sum_{\alpha}\int[\differential{{\bf\it
p}}]f_{0}(\epsilon_{0})-eE_{i}B_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\left(m_{j}^{\alpha}({\bf\it
p})+\partialderivative{m_{j}^{\alpha}}{p_{k}}mu_{k}\right)\partial_{\epsilon}f_{0}(\epsilon_{0})+O({\bf\it
B}^{2},{\bf\it B}{\bf\it u}^{2})$ $\displaystyle=eE_{i}n_{0}+O({\bf\it
B}^{2},{\bf\it B}{\bf\it u}^{2}),$ (S50) $\displaystyle
e\epsilon_{ikl}B_{l}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\tilde{v}_{k}^{\alpha}f_{\mathrm{local}}^{\alpha}$
$\displaystyle=e\epsilon_{ikl}B_{l}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\tilde{v}_{k}^{\alpha}({\bf\it p})f_{{\bf\it u}}^{\alpha}({\bf\it
p})+O({\bf\it B}^{2})$
$\displaystyle=e\epsilon_{ikl}B_{l}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\tilde{v}_{k}^{\alpha}({\bf\it p}+m{\bf\it u})f_{{\bf\it
u}}^{\alpha}({\bf\it p}+m{\bf\it u})+O({\bf\it B}^{2})$
$\displaystyle=e\epsilon_{ikl}B_{l}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\frac{p_{k}+mu_{k}}{m}f_{0}(\epsilon_{0})+O({\bf\it B}^{2},{\bf\it
B}{\bf\it u}^{2})$ $\displaystyle=e\epsilon_{ikl}B_{l}u_{k}n_{0}+O({\bf\it
B}^{2},{\bf\it B}{\bf\it u}^{2}),$ (S51) $\displaystyle
E_{k}B_{k}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\Omega_{i}^{\alpha}f_{\mathrm{local}}^{\alpha}$
$\displaystyle=E_{k}B_{k}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\Omega_{i}^{\alpha}({\bf\it p})f_{{\bf\it u}}^{\alpha}({\bf\it
p})+O({\bf\it B}^{2})$
$\displaystyle=E_{k}B_{k}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\Omega_{i}^{\alpha}({\bf\it p}+m{\bf\it u})f_{{\bf\it u}}^{\alpha}({\bf\it
p}+m{\bf\it u})+O({\bf\it B}^{2})$
$\displaystyle=E_{k}B_{k}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\left(\Omega_{i}^{\alpha}({\bf\it
p})+\partialderivative{\Omega_{i}^{\alpha}}{p_{j}}mu_{j}\right)f_{0}(\epsilon_{0})+O({\bf\it
B^{2}},{\bf\it B}{\bf\it u}^{2})$ $\displaystyle=O({\bf\it B^{2}},{\bf\it
B}{\bf\it u}^{2}).$ (S52)

After all, the driving force can be calculated as

$F_{i}=-en_{0}(E_{i}+\epsilon_{ikl}B_{l}u_{k})+O({\bf\it B}^{2},{\bf\it
B}{\bf\it u}^{2}),\quad\mathrm{i.e.}\quad{\bf\it F}=-en_{0}({\bf\it E}+{\bf\it
u}\times{\bf\it B})+O({\bf\it B}^{2},{\bf\it B}{\bf\it u}^{2}).$ (S53)

One recognizes that the driving force $F_{i}$ corresponds to the conventional
Lorentz force for a charged particle. At this stage, we assume the relaxation
time approximation for momentum-relaxating scatterings:

$C^{\mathrm{mr}}[f^{\alpha}]=-\frac{f^{\alpha}-f_{0}(\tilde{\epsilon}^{\alpha})}{\tau},$
(S54)

and we can evaluate the relaxing force as

$\Gamma_{i}=\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}D^{\alpha}({\bf\it
p})C^{\mathrm{mr}}[f^{\alpha}]=-\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}D^{\alpha}({\bf\it
p})\frac{f^{\alpha}-f_{0}(\tilde{\epsilon}^{\alpha})}{\tau}.$ (S55)

The first term on the right-hand side of Eq. (S55) can be calculated as

$\displaystyle\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}D^{\alpha}({\bf\it p})f_{\mathrm{local}}^{\alpha}$
$\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}D^{\alpha}({\bf\it p})f_{{\bf\it u}}^{\alpha}({\bf\it p})$
$\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it
p}}](p_{i}+mu_{i})\left(1+\frac{e}{\hbar}{\bf\it
B}\dotproduct{\bf\it\Omega}^{\alpha}({\bf\it p}+m{\bf\it u})\right)[f_{{\bf\it
u}}^{\alpha}({\bf\it p}+m{\bf\it u})-{\bf\it B}\dotproduct{\bf\it
m}^{\alpha}({\bf\it p}+{\bf\it u})\partial_{\epsilon}f_{{\bf\it
u}}^{\alpha}({\bf\it p}+m{\bf\it u})]$ $\displaystyle\quad+O({\bf\it B}^{2})$
$\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it
p}}](p_{i}+mu_{i})\left[1+\frac{e}{\hbar}B_{j}\left(\Omega^{\alpha}_{j}({\bf\it
p})+\partialderivative{\Omega_{j}^{\alpha}}{p_{k}}mu_{k}\right)\right]f_{0}(\epsilon_{0})$
$\displaystyle\quad-B_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}](p_{i}+mu_{i})\left[m^{\alpha}_{j}({\bf\it
p})+\partialderivative{m^{\alpha}_{j}}{p_{k}}mu_{k}\right]\partial_{\epsilon}f_{0}(\epsilon_{0})+O({\bf\it
B}^{2},{\bf\it B}{\bf\it u}^{2})$
$\displaystyle=mu_{i}\sum_{\alpha}\int[\differential{{\bf\it
p}}]f_{0}(\epsilon_{0})+\frac{e}{\hbar}B_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}\Omega_{j}^{\alpha}f_{0}(\epsilon_{0})-B_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}m^{\alpha}_{j}({\bf\it
p})\partial_{\epsilon}f_{0}(\epsilon_{0})+O({\bf\it B}^{2},{\bf\it B}{\bf\it
u}^{2})$
$\displaystyle=mu_{i}n_{0}+\frac{e}{\hbar}B_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}\Omega_{j}^{\alpha}f_{0}(\epsilon_{0})-B_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}m^{\alpha}_{j}({\bf\it
p})\partial_{\epsilon}f_{0}(\epsilon_{0})+O({\bf\it B}^{2},{\bf\it B}{\bf\it
u}^{2}).$ (S56)

On the other hand, the second term on the right-hand side of Eq. (S55) can be
calculated as

$\displaystyle\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}D^{\alpha}({\bf\it p})f_{0}(\tilde{\epsilon}^{\alpha})$
$\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}D^{\alpha}({\bf\it p})[f_{0}(\epsilon_{0})-{\bf\it
B}\dotproduct{\bf\it m}^{\alpha}({\bf\it
p})\partial_{\epsilon}f_{0}(\epsilon_{0})]+O({\bf\it B}^{2})$
$\displaystyle=\frac{e}{\hbar}B_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}\Omega^{\alpha}_{j}f_{0}(\epsilon_{0})-B_{j}\sum_{\alpha}\int[\differential{{\bf\it
p}}]p_{i}m^{\alpha}_{j}({\bf\it
p})\partial_{\epsilon}f_{0}(\epsilon_{0})+O({\bf\it B}^{2}).$ (S57)

Consequently, the momentum-relaxing force $\Gamma_{i}$ can be expressed in
terms of hydrodynamic variables as follows:

$\Gamma_{i}=-\frac{mu_{i}n_{0}}{\tau}.$ (S58)

##  V Derivation of Generalized Euler equation

Combining the particle conservation law Eq. (S4) and the momentum conservation
law Eq. (S5), we can derive the hydrodynamic equation which is a
generalization of conventional Euler equation correct up to ${\bf\it B}^{2}$,
${\bf\it B}{\bf\it u}^{2}$, ${\bf\it u}^{3}$. First, the particle conservation
law

$\partial_{t}n_{{\bf\it B}}+\divergence{\bf\it J}^{n}_{{\bf\it B}}=0$ (S59)

can be expressed by hydrodynamic variables of the order of $O({\bf\it
B}^{2},{\bf\it B}{\bf\it u}^{2},{\bf\it u}^{3})$,

$\partialderivative{t}\left[n_{0}+\frac{me}{\hbar}u_{i}B_{j}(D_{ij}+N_{ij})\right]+\partial_{i}\left[n_{0}u_{i}+\frac{me}{\hbar}\epsilon_{ijk}E_{j}(D_{lk}u_{l}+Y_{lk}B_{l})\right]=0.$
(S60)

In a similar way, the momentum conservation law

$\partial_{t}P_{{\bf\it B},i}+\partial_{j}\Pi_{{\bf\it
B},ij}=F_{i}+\Gamma_{i}$ (S61)

can be expressed in terms of hydrodynamic variables as follows:

$\displaystyle\partialderivative{t}\left[mn_{0}u_{i}+\frac{e}{\hbar}B_{j}(C_{ij}+M_{ij})\right]+\partial_{j}\left[mn_{0}u_{i}u_{j}+p\delta_{ij}+\frac{e}{\hbar}\epsilon_{jkl}E_{k}C_{il}+\frac{e}{\hbar}B_{k}(u_{l}M_{lk}\delta_{ij}+u_{j}M_{ik})+\frac{e}{\hbar}B_{j}u_{l}O_{il}\right]$
$\displaystyle=-en_{0}(E_{i}+\epsilon_{ikl}u_{k}B_{l})-\frac{mn_{0}u_{i}}{\tau}.$
(S62)

Substituting Eq. (S60) into Eq. (S62), we obtain the generalized Euler
equation

$\displaystyle mn_{0}\partialderivative{u_{i}}{t}+mn_{0}({\bf\it
u}\dotproduct\gradient)u_{i}+\partial_{i}p+en_{0}(E_{i}+\epsilon_{ikl}u_{k}B_{l})+\frac{e}{\hbar}M_{ij}\partialderivative{B_{j}}{t}+\frac{e}{\hbar}\epsilon_{jkl}E_{k}\partial_{j}C_{il}+\frac{e}{\hbar}B_{k}\partial_{j}(u_{l}M_{lk}\delta_{ij}+u_{j}M_{ik})$
$\displaystyle+\frac{e}{\hbar}B_{j}\partial_{j}(u_{l}O_{il})=-\frac{mn_{0}u_{i}}{\tau}.$
(S63)

If we omit the spatial derivatives of $T$ and $\mu$ for simplicity, we obtain
the Euler equation which is shown in the main text:

$\displaystyle mn_{0}\partialderivative{u_{i}}{t}+mn_{0}({\bf\it
u}\dotproduct\gradient)u_{i}+\partial_{i}p+en_{0}(E_{i}+\epsilon_{ikl}u_{k}B_{l})+\frac{e}{\hbar}M_{ij}\partialderivative{B_{j}}{t}+\frac{e}{\hbar}B_{k}(M_{lk}\partial_{i}u_{l}+M_{ik}\partial_{j}u_{j})+\frac{e}{\hbar}B_{j}O_{il}\partial_{j}u_{l}$
$\displaystyle=-\frac{mn_{0}u_{i}}{\tau}.$ (S64)

##  VI local current

To relate the hydrodynamic theory with an observable current in optical
experiments, we further need to introduce the so-called local current [6] when
we consider the magneto-optical responses,

$\displaystyle{\bf\it j}=-e\sum_{\alpha}\int[\differential{{\bf\it
p}}]D^{\alpha}\dot{{\bf\it
r}}^{\alpha}f^{\alpha}+\curl\sum_{\alpha}\int[\differential{{\bf\it
p}}]{\bf\it m}^{\alpha}f^{\alpha}.$ (S65)

The detailed calculation of the second term of the right-hand side is given by

$\displaystyle\sum_{\alpha}\int[\differential{{\bf\it
p}}]m_{k}^{\alpha}f^{\alpha}_{\mathrm{local}}$
$\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it
p}}]m_{k}^{\alpha}({\bf\it p})[f^{\alpha}_{{\bf\it u}}({\bf\it p})-{\bf\it
B}\dotproduct{\bf\it m}^{\alpha}({\bf\it p})\partial_{\epsilon}f_{{\bf\it
u}}^{\alpha}({\bf\it p})]+O({\bf\it B}^{2})$
$\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it
p}}]m_{k}^{\alpha}({\bf\it p}+m{\bf\it u})[f^{\alpha}_{{\bf\it u}}({\bf\it
p}+m{\bf\it u})-{\bf\it B}\dotproduct{\bf\it m}^{\alpha}({\bf\it p}+m{\bf\it
u})\partial_{\epsilon}f_{{\bf\it u}}^{\alpha}({\bf\it p}+m{\bf\it
u})]+O({\bf\it B}^{2})$ $\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it
p}}]\left(m_{k}^{\alpha}({\bf\it
p})+\partialderivative{m_{k}^{\alpha}}{p_{l}}mu_{l}\right)f_{0}(\epsilon_{0})$
$\displaystyle\quad-B_{l}\sum_{\alpha}\int[\differential{{\bf\it
p}}]\left(m_{k}^{\alpha}({\bf\it
p})+\partialderivative{m_{k}^{\alpha}}{p_{m}}mu_{m}\right)\left(m_{l}^{\alpha}({\bf\it
p})+\partialderivative{m_{l}}{p_{n}}mu_{n}\right)\partial_{\epsilon}f_{0}(\epsilon_{0})+O({\bf\it
B}^{2},{\bf\it B}{\bf\it u}^{2})$
$\displaystyle=\sum_{\alpha}\int[\differential{{\bf\it
p}}]\partialderivative{m_{k}^{\alpha}}{p_{l}}mu_{l}f_{0}(\epsilon_{0})-B_{l}\int[\differential{{\bf\it
p}}]m_{k}^{\alpha}m_{l}^{\alpha}\partial_{\epsilon}f_{0}(\epsilon_{0})+O({\bf\it
B}^{2},{\bf\it B}{\bf\it u}^{2}).$ (S66)

Therefore, we obtain the local current as follows:

$\displaystyle{\bf\it j}=-en_{0}{\bf\it u}-\frac{me^{2}}{\hbar}{\bf\it
E}\times({}^{t}\hat{D}{\bf\it u}+{}^{t}\hat{Y}{\bf\it
B})+\frac{e}{\hbar}\gradient\times({}^{t}\hat{M}{\bf\it
u})-\gradient\times(\hat{W}{\bf\it B}),$ (S67)

with $W_{kl}=\int[\differential{{\bf\it
p}}]m_{k}^{\alpha}m_{l}^{\alpha}\partial_{\epsilon}f_{0}(\epsilon_{0})$. In
the main text, the local current is described as

${\bf\it j}=-en_{0}{\bf\it u}-\frac{me^{2}}{\hbar}{\bf\it
E}\times({}^{t}\hat{D}{\bf\it u}+{}^{t}\hat{Y}{\bf\it
B})+\frac{e}{\hbar}\gradient\times({}^{t}\hat{M}{\bf\it u}),$ (S68)

where the terms with spatial derivatives of $T$ and $\mu$ are omitted for
simplicity.

##  VII linearization of generalized Euler equation

Next we move on to the analyses of collective modes. In the study of
collective modes, the deviations of the local thermodynamic parameters from
their equilibrium values are small. Then, the use of linearized hydrodynamic
equations is sufficient. They are obtained by looking for a solution in the
form of plane waves, i.e. ${\bf\it E}=\tilde{{\bf\it E}}e^{i{\bf\it
k}\dotproduct{\bf\it r}-i\omega t}$ together with similar expressions for
other oscillating variables such as $n$, ${\bf\it u}$ and ${\bf\it B}$. Here
and hereafter, the subscript ${\bf\it B}$ is omitted for simplicity and the
subscript 0 describes the non-oscillating part of physical variables.

First, we linearize the equation of particle conservation law. Substituting
$n=n_{0}+\tilde{n}e^{i{\bf\it k}\dotproduct{\bf\it r}-i\omega t}$, ${\bf\it
u}=\tilde{{\bf\it u}}e^{i{\bf\it k}\dotproduct{\bf\it r}-i\omega t}$ and
${\bf\it B}={\bf\it B}_{0}+\tilde{{\bf\it B}}e^{-i\omega t}$ into the particle
conservation law

$\partialderivative{n}{t}+\divergence{\bf\it J}^{n}=0,$ (S69)

we obtain the following equation:

$\tilde{n}=\frac{{\bf\it k}\dotproduct\tilde{{\bf\it
u}}}{\omega}n_{0}+\frac{e}{\hbar}\frac{({\bf\it k}\times\tilde{{\bf\it
E}})_{k}}{\omega}Y_{lk}B_{0}^{l}.$ (S70)

Next, we linearize the hydrodynamic equation in a similar way and obtain the
linearized hydrodynamic equation:

$\displaystyle
mn_{0}(-i\omega\tilde{u}_{i})+\frac{K}{n_{0}}i{k}_{i}\tilde{n}+en_{0}(\tilde{{E}}_{i}+\epsilon_{ijk}\tilde{{u}}_{j}{B}^{k}_{0})-\frac{e}{\hbar}M_{ik}(i\epsilon_{klj}{k}_{l}\tilde{{E}}_{j})+\frac{e}{\hbar}B_{0}^{k}M_{lk}\tilde{u}_{l}(i{k}_{i})+\frac{e}{\hbar}(i{k}_{j}\tilde{{u}}_{j})M_{ik}{B}_{0}^{k}+\frac{e}{\hbar}(i{k}_{k}{B}^{k}_{0})O_{ij}\tilde{{u}}_{j}$
$\displaystyle=-\frac{mn_{0}}{\tau}\tilde{{u}}_{i},$ (S71)

where we use the relation between compressibility $K$ and particle density $n$

$K\equiv\frac{1}{\kappa}=-V\left(\partialderivative{p}{V}\right)=n\partialderivative{p}{n},\quad\mathrm{i.e.}\quad\gradient
p=\partialderivative{p}{n}\gradient n=\frac{K}{n}\gradient
n\simeq\frac{K}{n_{0}}i{\bf\it k}\tilde{n}.$ (S72)

Substituting $\tilde{n}$ of the form in Eq.(S70) into Eq.(S71), we obtain

$\displaystyle\left[\left(-i\omega+\frac{1}{\tau}\right)\delta_{ij}+\frac{e}{m}\epsilon_{ijk}B_{0}^{k}+\frac{e}{\hbar}\frac{iB_{0}^{k}}{mn_{0}}\left(M_{jk}k_{i}+M_{ik}k_{j}+O_{ij}k_{k}\right)+\frac{iK}{mn_{0}\omega}k_{i}k_{j}\right]\tilde{u}_{j}$
$\displaystyle=\left[-\frac{e}{m}\delta_{ij}+\frac{e}{\hbar}\frac{iM_{ik}}{mn_{0}}\epsilon_{jkl}k_{l}-\frac{e}{\hbar}\frac{iK}{n_{0}^{2}\omega}Y_{mk}B_{0}^{m}\epsilon_{jkl}k_{i}k_{l}\right]\tilde{E}_{j}.$
(S73)

Solving the above simultaneous equations for $\tilde{{\bf\it u}}$ and
substituting it into Eq. (S68), we obtain the optical conductivity
$\hat{\sigma}(\omega,{\bf\it k})$.

##  VIII nonreciprocity in surface magnetoplasmons

In order to study the nonreciprocal collective modes from hydrodynamic
equations, one notices that the inverse compressibility $K\equiv 1/\kappa$ is
not important for the nonreciprocity because it appears in linearized
hydrodynamic equations with $k_{i}k_{j}$, so we ignore the term proportional
to $K$ in the following analyses. According to symmetry consideration [7], the
geometrical tensor $\hat{M}$ vanishes in a pseudo-vector constrained in the 2D
plane : $M_{ij}=M_{i}\delta_{jz}$. Here we take the vector ${\bf\it M}$ along
$y$-axis without loss of generality. The $x$ and $y$ component of Eq. (S73)
are given by

$\left(\begin{array}[]{cc}1-i\omega\tau&\omega_{\mathrm{c}}\tau(1+i\zeta
q_{x})\\\ &\\\ -\omega_{\mathrm{c}}\tau(1-i\zeta
q_{x})&1-i\omega\tau+2i\omega_{\mathrm{c}}\tau\zeta
q_{y}\end{array}\right)\left(\begin{array}[]{c}\tilde{u}_{x}\\\ \\\
\tilde{u}_{y}\end{array}\right)=-\frac{e\tau}{m}\left(\begin{array}[]{cc}1&0\\\
&\\\ i\zeta q_{y}&1-i\zeta
q_{x}\end{array}\right)\left(\begin{array}[]{c}\tilde{E}_{x}\\\ \\\
\tilde{E}_{y}\end{array}\right),$ (S74)

respectively with $\zeta=M_{y}/\hbar n_{0}$. Solving this equation, we obtain

$\displaystyle\left(\begin{array}[]{c}\tilde{u}_{x}\\\ \\\
\tilde{u}_{y}\end{array}\right)$
$\displaystyle=-\frac{e\tau/m}{Z(\omega,{\bf\it
q})}\left(\begin{array}[]{cc}1-i\omega\tau+2i\omega_{\mathrm{c}}\tau\zeta
q_{y}-\omega_{\mathrm{c}}\tau(1+i\zeta q_{x})i\zeta
q_{y}&-\omega_{\mathrm{c}}\tau(1+\zeta^{2}q_{x}^{2})\\\ &\\\
\omega_{\mathrm{c}}\tau(1-i\zeta q_{x})+(1-i\omega\tau)i\zeta
q_{y}&(1-i\omega\tau)(1-i\zeta
q_{x})\end{array}\right)\left(\begin{array}[]{c}\tilde{E}_{x}\\\ \\\
\tilde{E}_{y}\end{array}\right).$ (S84)

Here, $Z(\omega,{\bf\it q})=(1-i\omega\tau+2i\omega_{\mathrm{c}}\tau\zeta
q_{y})(1-i\omega\tau)+(\omega_{\mathrm{c}}\tau)^{2}(1+\zeta^{2}q_{x}^{2})$.
Using the form of local current in terms of hydrodynamic valuables:

$\tilde{j}_{i}=-en_{0}\tilde{u}_{i}-\frac{me^{2}}{\hbar}\epsilon_{ijk}Y_{lk}B_{0}^{l}\tilde{E}_{j}+\frac{e}{\hbar}\epsilon_{ijk}iq_{j}M_{lk}\tilde{u}_{l},$
(S85)

we obtain the optical conductivity given by

$\displaystyle\sigma_{xx}$
$\displaystyle=\sigma_{\mathrm{D}}(1-i\omega\tau)(1+\zeta^{2}q_{y}^{2})/Z,$
(S86) $\displaystyle\sigma_{xy}$ $\displaystyle=-\sigma_{\mathrm{D}}(1-i\zeta
q_{x})[\omega_{\mathrm{c}}\tau(1+i\zeta q_{x})+i\zeta
q_{y}(1-i\omega\tau)]/Z-\frac{me^{2}}{\hbar}Y_{zz}B_{0},$ (S87)
$\displaystyle\sigma_{yx}$ $\displaystyle=\sigma_{\mathrm{D}}(1+i\zeta
q_{x})[\omega_{\mathrm{c}}\tau(1-i\zeta q_{x})+i\zeta
q_{y}(1-i\omega\tau)]/Z+\frac{me^{2}}{\hbar}Y_{zz}B_{0},$ (S88)
$\displaystyle\sigma_{yy}$
$\displaystyle=\sigma_{\mathrm{D}}(1-i\omega\tau)(1+\zeta^{2}q_{x}^{2})/Z.$
(S89)

Here, $\sigma_{\mathrm{D}}$ is the Drude weight and ${\bf\it q}=(q_{x},q_{y})$
is the in-plane wave vector of the incident light.

In the following analysis, we consider surface plasmons in 2D electron
hydrodynamic materials such as monolayer graphene. Surface plasmons originate
from the collective excitations of electrons coupled to the electromagnetic
field at an interface between a dielectric material and a metal. The equation
which determines the dispersion relation of surface plasmons is given by

$f(\omega,{\bf\it
q})=\left[\frac{\epsilon_{1}}{\kappa_{1}}+\frac{\epsilon_{2}}{\kappa_{2}}+\frac{i\sigma_{xx}}{\epsilon_{0}\omega}\right]\left[\kappa_{1}+\kappa_{2}-\frac{i\omega\sigma_{yy}}{\epsilon_{0}c^{2}}\right]-\frac{\sigma_{xy}\sigma_{yx}}{(\epsilon_{0}c^{2})}=0.$
(S90)

Here, $\kappa_{i}=\sqrt{q^{2}-\epsilon_{i}(\omega/c)^{2}}$ $(i=1,2)$,
$\epsilon_{i}$ are the relative dielectric constants in the medium above
$(i=1)$ and below $(i=2)$ the graphene, $\epsilon_{0}$ and $c$ are the
dielectric constant and the speed of light in free space respectively.

This equation $f(\omega,{\bf\it q})=0$ is difficult to solve analytically with
respect to $\omega=\omega({\bf\it q})$, however, we can know the
nonreciprocity of the surface magnetoplasmons without solving the equation.
Let us consider two different directions of propagation.

  * •

If we consider the surface plasmons propagating along the $x$-direction, the
in-plane wave vector is ${\bf\it q}=(q_{x},0)$ and Eq.(S90) becomes
$f(\omega,q_{x}^{2})=0$. This can be solved as $\omega=\omega(q_{x}^{2})$ and
we conclude

$\omega(q_{x})=\omega(-q_{x}).$ (S91)

  * •

If we consider the surface plasmons propagating along the $y$-direction, the
in-plane wave vector is ${\bf\it q}=(0,q_{y})$ and $f(\omega,{\bf\it q})$
generally satisfies

$f(\omega,q_{y})-f(\omega,-q_{y})\neq 0$ (S92)

due to the nonreciprocity of $\sigma_{ij}$, i.e.,
$\sigma_{ij}(\omega,q_{y})\neq\sigma_{ij}(\omega,-q_{y})$. Substituting
$\omega=\omega(-q_{y})$ in Eq.(S92), we obtain

$f(\omega(-q_{y}),q_{y})-f(\omega(-q_{y}),-q_{y})\neq 0.$ (S93)

Using the fact that $\omega=\omega(-q_{y})$ is the solution of
$f(\omega,-q_{y})=0$ and satisfies $f(\omega(-q_{y}),-q_{y})=0$, we obtain
from Eq. (S93)

$f(\omega(-q_{y}),q_{y})\neq 0.$ (S94)

However, $\omega(q_{y})$ satisfies $f(\omega(q_{y}),q_{y})=0$, thus we
conclude

$\omega(q_{y})\neq\omega(-q_{y}).$ (S95)

In summary, the surface magnetoplasmons acquire nonreciprocity only when they
propagate along ${\bf\it M}$. On the other hand, the surface magnetoplasmons
propagating perpendicular to ${\bf\it M}$ are reciprocal.

## References

  * Xiao  _et al._ [2010] D. Xiao, M.-C. Chang, and Q. Niu, Rev. Mod. Phys. 82, 1959 (2010).
  * Lifshitz and Pitaevskii [1981] E. M. Lifshitz and L. P. Pitaevskii,  _Course of Theoretical Physics: Physical Kinetics_ (Pergamon, New York, 1981).
  * Lucas and Fong [2018] A. Lucas and K. C. Fong, Journal of Physics: Condensed Matter 30, 053001 (2018).
  * Tavakol and Kim [2021] O. Tavakol and Y. B. Kim, Phys. Rev. Research 3, 013290 (2021).
  * Sodemann and Fu [2015] I. Sodemann and L. Fu, Phys. Rev. Lett. 115, 216806 (2015).
  * Xiao  _et al._ [2006] D. Xiao, Y. Yao, Z. Fang, and Q. Niu, Phys. Rev. Lett. 97, 026603 (2006).
  * Toshio  _et al._ [2020] R. Toshio, K. Takasan, and N. Kawakami, Phys. Rev. Research 2, 032021 (2020).


# Efficient and Accurate Automatic Python Bindings with cppyy & Cling

Baidyanath Kundu$\ast$$\dagger$ Vassil Vassilev$\ast$$\dagger$ Wim
Lavrijsen$\S$ $\ast$European Council for Nuclear Research, Espl. des
Particules 1, 1211 Meyrin, Switzerland
$\dagger$Department of Physics, Princeton University, Princeton, New Jersey
08544, USA
$\S$LBNL, MS 50B-6222, 1 Cyclotron Road, Berkeley CA 94720-8147, USA
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The simplicity of Python and the power of C++ force stark choices on a
scientific software stack. There have been multiple developments to mitigate
language boundaries by implementing language bindings, but the impedance
mismatch between the static nature of C++ and the dynamic one of Python
hinders their implementation; examples include the use of user-defined Python
types with templated C++ and advanced memory management.

The development of the C++ interpreter Cling has changed the way we can think
of language bindings as it provides an incremental compilation infrastructure
available at runtime. That is, Python can interrogate C++ on demand, and
bindings can be lazily constructed at runtime. This automatic binding
provision requires no direct support from library authors and offers better
performance than alternative solutions, such as PyBind11. ROOT pioneered this
approach with PyROOT, which was later enhanced with its successor, cppyy.
However, until now, cppyy relied on the reflection layer of ROOT, which is
limited in terms of provided features and performance.

This paper presents the next step for language interoperability with cppyy,
enabling research into uniform cross-language execution environments and
boosting optimization opportunities across language boundaries. We illustrate
the use of advanced C++ in Numba-accelerated Python through cppyy. We outline
a path forward for re-engineering parts of cppyy to use upstream LLVM
components to improve performance and sustainability. We demonstrate cppyy
purely based on a C++ reflection library, InterOp, which offers
interoperability primitives based on Cling and Clang-Repl.

##  1 Introduction

The C++ programming language was adopted in the mid-‘90s by the high energy
physics (HEP) community to provide a more modern yet performant programming
environment. At the time, HEP had some of the world’s largest data sets (and
processing rates). Python usage in HEP started in the mid-2000s, initially
only for interactive access to experiments’ frameworks and analysis.
Meanwhile, the rest of the world caught up with and even surpassed HEP data
processing needs, and a complete ecosystem has grown up around Python to
support data science. Much of the modern software for HEP physics analysis,
such as Awkward Array [1], Uproot [2], and Coffea [3], is grounded in this
Python ecosystem. C++ remains favored because of its performance, access to
accelerators in heterogeneous computing environments, and having foundational
libraries in HEP, such as ROOT [4], GEANT4 [5], and most experiments’
processing and production software frameworks. With both languages playing
important roles, Python-C++ integration has become instrumental. However, more
advanced integration is necessary for high-performance codes and codes that
run in heterogeneous environments. For example, a data processing task
performed by a C++ framework running a user-provided Python function or kernel
incurs significant overhead because of the large number of language crossings.
Similarly, if the task is run as part of a GPU workflow, it will incur a
performance penalty due to offloading since the Python code can only run on
the CPU.

Numba [6], a just-in-time compiler (JIT) for Python code, addresses this
problem: it is capable of compiling Python code, targeting either the CPU or
GPU, and provides callable interfaces to use the JITed closures from low-level
libraries. Thus, Numba greatly improves the performance of Python code not
only by lowering it to machine code but also by removing unnecessary crossings
of language (or device) boundaries in inner loops. However, Numba has two
shortcomings: preparing external code, such as functions written in C, to be
usable for the Numba JIT is an involved, manual process; and Numba does not
support C++ directly.

Cppyy, an automatic runtime bindings generator [7] based on the C++
interpreter Cling [8], can address both these shortcomings. One goal of our
work is to bring advanced C++, in particular highly optimized numerics
libraries, into Numba-accelerated Python in a fully automatic way, with
minimal language crossing overhead. In this paper, we present a fully generic
prototype that provides, lazily and at runtime, Numba extensions for bound C++
through cppyy. Our prototype supports overloading, C++ templates, data member
access, and instance method calls in Numba-accelerated Python. We outline a
path forward for re-engineering cppyy’s backend to use LLVM components more
directly, thus improving performance and sustainability. We demonstrate cppyy
on InterOp, a new library that implements interoperability primitives based on
Cling and Clang-Repl.

This paper is organized as follows: Section 2 motivates our technology and
design choices; Section 3 provides a bird’s eye overview of the implementation
design; Section 4 shows timing results with the current implementation;
Section 5 discusses the enhancements brought about by a new reflection
library, InterOp, which will allow for even better performance and greater
flexibility when retargeting codes; Section 6 discusses related work; finally,
Section 7 summarizes our findings.

##  2 Motivation

The typical approach to speed up a Python application has been to write the
performance-critical portions of it in a lower-level language, usually C, and
then access that code in the Python interpreter through an extension module. A
much more elegant solution is to lower the original Python code to native
through the use of a JIT so that the developer can stay within a single
environment, Python, and only needs to write and debug the code once. In HEP,
it is critical that such a JIT works well with bound C++ code and that JITed
Python code is usable from C++, since any full HEP software stack relies
heavily on both languages. For this paper, we use Numba as the Python JIT and
integrate C++ into it through cppyy.

Numba has an extension application programming interface (API) to add
functions and types for the JIT to recognize, annotate, and lower to LLVM IR
(intermediate representation; which in turn gets assembled to native code).
Manually writing such extensions is tedious and time-consuming, so automation
is critically important, especially to meet the goal of being able to use the
JIT transparently. Integrating C++ with Numba would additionally make it easy
to use Python kernels in C++ code without losing performance, thus reaping the
benefits of continued use of a large installed base of C++ codes. To enable
this, we use cppyy to generate first-class Python objects for bound C++, which
enables Numba to query for reflection information through regular Python
introspection and for cppyy to provide the necessary Numba extensions lazily
at runtime. The result is a generic, compact, and efficient implementation.

Modern C++ is especially amenable to automation: successive C++ standards have
moved the language away from low-level code to more expressive abstractions
for the purposes of code quality and optimizing compilers. In modern C++, many
previously difficult to analyze corner cases are now more clear; compare,
e.g., ownership rules for raw v.s. smart pointers. Furthermore, cppyy is based
on Cling, the C++ interpreter, to have C++ match Python’s dynamism,
interactivity, and runtime behavior. Cling itself is based on Clang/LLVM,
which allows it to stay up to date even as the C++ language continues to
evolve rapidly. The use of Cling also means that unresolved corner cases can
be handled at runtime in either C++ or Python as appropriate, without
necessitating an intermediate bridging or mapping language.

##  3 Implementation

We introduce C++ into Numba through a new reflection interface on top of
cppyy, which resolves overloaded functions and instantiates C++ templates
during Numba’s type annotation step; and which provides the low-level
description for lowering to LLVM IR for subsequent optimization and conversion
of assembly to machine code. This process is completely automatic and
transparent, including integration with other manually provided Numba
extensions. Performance is on par, with only a moderate additional warm-up
cost. Python programmers can thus continue to develop and debug their codes
fully in Python while simply switching on the Numba JIT for selected
performance-critical closures.

Python is dynamically typed. It achieves this by wrapping objects into so-
called PyObjects via a method known as boxing and retrieves their original
value via a method known as unboxing. Since Numba needs to eliminate these
Python boxing and unboxing steps, it requires accurate typing information
ahead of time. To provide detailed, exact and low-level C++ type information
to Numba, which is needed in addition to the type information already
available through Python’s introspection, we introduce a new Reflection API to
cppyy.

The reflection API follows the request/reply pattern and adds a new function
to all cppyy-bound objects called

    __cpp_reflex__}. The function has a constant signature for all cppyy object types, but its functionality depends on which one it is called on. It takes two arguments, the type of reflection information and its format, and returns the information requested. The request can range from asking if the object represents a namespace or class to requesting the object’s C++ type. The format argument specifies whether we want the information in a string format, a C++ type, or in the manner that the API thinks is the most optimal.
    \begin

figure*[h]  The interaction diagram for Numba, numba_ext and cppyy Fig. 3
shows the layered approach taken by Numba to separate concerns. Numba analyzes
a Python closure as directed, and upon encountering cppyy types, it queries
cppyy’s pre-registered numba_ext module for exact type information. This
module, in turn, queries cppyy’s new reflection API if it encounters a type
that has not been seen before. The module then uses Numba’s extension and core
modules to generate the necessary typing classes, to be returned to Numba
during this typing step; and lowering methods, which Numba will expect in its
next step, to write the LLVM IR. Each core language construct (free functions,
namespaces, classes, methods, data members, etc.) has its own implementation.
For instance, to support C++ free functions, we extend the Numba Callable type
and register it as the type provider for cppyy function objects. When Numba
encounters a C++ function bound by cppyy, it queries this registered type for
the function signature (i.e., the return and argument types), given the actual
argument types traced so far. Cppyy’s numba_ext uses the argument types to
select an appropriate overload, aided by the new reflection API, and wraps
this overload in a unique typing class. Finally, numba_ext registers a
lowering method uniquely generated for the address of the chosen C++ function
overload, which will be called during the Numba lowering step. This process
provides Numba with all the information it needs to convert the function call
to LLVM IR.

##  4 Results

We benchmark the running time of Numba JITed functions with cppyy objects
against their Python counterparts to obtain the time taken by Numba to JIT the
function (Numba JIT time), the time taken by cppyy to create the typing info
and possibly perform lookups and instantiate templates (cppyy JIT time), the
time taken to run the function after it has been JITed (Hot run time), and the
time taken to run the equivalent Python function. The results are obtained on
a 3.1GHz Intel NUC Core i7-8809G CPU with 32G RAM.

⬇

\parcppyy.cppdef(r”””

template<class T>

T add42(T t) {

return T(t+42);

}

”””)

\par

C++ code in cppyy

⬇

|\phantom{empty space}|

def go_slow(a):

trace = 0.0

for i in range(a.shape[0]):

trace +=

|\colorbox{higlightcolor}{cppyy.gbl.add42(a[i, i]) + }|

|\colorbox{higlightcolor}{cppyy.gbl.add42(int(a[i, i]))}|

return a + trace

Python/C++ binding with cppyy

⬇

@numba.jit(nopython=True)

def go_fast(a):

trace = 0.0

for i in range(a.shape[0]):

trace +=

|\colorbox{higlightcolor}{cppyy.gbl.add42(a[i, i]) + }|

|\colorbox{higlightcolor}{cppyy.gbl.add42(int(a[i, i]))}|

return a + trace

Numba-acceleration of cppyy

Figure 1: Benchmark fixture for ‘Templated free functions’ case. On the left
side, the C++ templated function is declared in cppyy. In the center, we use a
Python kernel to run this C++ function. On the right, we use the same kernel
but add the Numba JIT decorator to accelerate it and compare the timing of the
two kernels. The other cases use the same setup, with only the highlighted
code being replaced by the case being benchmarked.

To evaluate the speedup obtained by Numba JITing of cppyy objects, we use the
fixture in Fig. 1. For each benchmark case in Tab. 1, a Numpy array of size
$100\times 100$ was passed to the function; times indicated are averages of
3000 runs. Numba JITed functions achieve a minimum speedup of 2.3 times in the
case of methods and a maximum speedup of nearly 21 times in the case of
templated free functions. Tab. 1 summarizes timings obtained on our chosen
benchmarks.

Benchmark Case |  | Cppyy JIT
---
time (s)
| Numba JIT
---
time (s)
| Hot run
---
time (s)
| Python run
---
time (s)
Speedup
Function w/o args | 1.72e-03 | 3.33e-01 | 3.58e-06 | 1.73e-05 | $4.84\times$
Overloaded functions | 1.05e-03 | 1.35e-01 | 4.51e-06 | 3.47e-05 | $7.70\times$
Templated free functions | 8.92e-04 | 1.45e-01 | 3.48e-06 | 7.18e-05 | $20.66\times$
Class data members | 1.43e-06 | 1.33e-01 | 5.87e-06 | 1.82e-05 | $3.10\times$
Class methods | 2.16e-03 | 1.39e-01 | 6.06e-06 | 1.43e-05 | $2.36\times$
Table 1: Performance comparison. The columns cppyy JIT time and Numba JIT time
represent respectively the time taken by cppyy to create the Python objects
for C++ entities and time taken for Numba to compile the Python function. Both
do not execute the compiled functions. Column Hot run time shows the time
taken to run the Numba JITed function, whereas Python run time shows the time
taken to run the uncompiled Python function with cppyy code. Column Speedup
compares the Hot run time to Python run time.

##  5 Further Improvements

A promising area of research is incorporating the LLVM IR of the C++ code
directly into the Numba-generated LLVM IR, in effect inlining the C++
functions into Python. This would open up more optimization opportunities for
the compiler and simplify support for alternate targets, such as GPUs. In
practice, this requires IR to be sent from Cling to Numba or vice versa. One
of the challenges here is that Numba uses the llvmlite [9] package to
synthesize LLVM IR, while Cling operates directly on the underlying IR. Since
LLVM IR has no guarantees of backward or forward compatibility, this could
make a rather fragile dependency. Furthermore, the current design of Cling
does not provide the right abstractions and does not have the necessary low-
level infrastructure to either extract or inject IR into the interactive C++
environment. We therefore created a new, specialized library called InterOp to
mitigate these problems. This library serves as a high-performance interface
between C++ and Python, abstracting out details from the underlying compiler
API. The InterOp library is built on top of Clang-Repl, a new capability
available in LLVM, and has a reflection-focused API surface. InterOp will
follow LLVM’s release cycles, enabling its client to quickly adopt new
features.

Figure 2: Time taken and memory used during class template instantiation. On
the left, we compare template instantiations with `w`here the number of
template instantiations done by the C++ interpreter increases with the number
of template arguments. On the right, we compare instantiating nested
templates, for example, `,` where cppyy has to instantiate each nesting
individually from the innermost to the outermost class template. These are
common features of high-performance, templated numerics libraries that utilize
template expressions.

The design of InterOp brings significant improvements to cppyy in the time
taken and memory used for template instantiations. It shows promising
preliminary results thanks to the new way of modeling them. In Fig. 5 (left),
we show improvements in unpacking templates where we recursively instantiate a
textbook example of a multitype array with an underlying implementation based
on `std::tuple}s. The InterOp version is about $40\%$ faster and $4.5\%$ more
memory efficient. In Fig.~\ref`fig:InterOpBenchmarks (right), we show that
instantiating deeply nested templates scales better. The initial speedup of
6.2 times decreases to 3.8 times at about 4 levels of nesting but improves
rapidly after that. In addition, we reduce the total memory used. In addition
to the improved performance, the engineering merits of using Clang-Repl and
InterOp in cppyy are: (a) adopting new features faster, including potential
CUDA support; (b) better runtime and memory performance due to using opaque
data structures directly from the underlying compiler rather than their string
representations; and (c) a wider range of supported stock LLVM versions.

##  6 Related Work

Automatically improving the performance of Python code, i.e. without having to
rewrite it, has been a research and/or engineering topic for many projects.
Chosen Methodologies include just-in-time compilation [10][11], transpilation
followed by ahead-of-time compilation [12], or by using Python as a domain-
specific language [13][14]. Just-in-time compilation provides the most
pythonic experience, however, since it can be used interactively, and is very
relevant today as it is the most suitable in heterogeneous computing
environments. The PyPy project implements an alternative Python interpreter
that uses a tracing JIT to improve performance. PyPy is fully compatible with
the standard CPython implementation and is able to provide significant
speedups, especially for numeric code. Cppyy is integrated at the interpreter
level into PyPy, which allows C++ code to be called directly from the JITed
traces without the need for function wrappers, and which enables the inlining
of access to C++ data. PyPy’s tracing JIT is, contrary to all other options,
not based on LLVM. This reduces dependencies and overheads, but it does mean
that the JIT can target fewer platforms and offers fewer optimizations. A weak
point of PyPy is support for existing CPython extension modules, access to
which, because of different memory models, can be rather slow. Furthermore,
because the memory models need to be mapped, it is not as lenient as CPython
when it comes to reference counting bugs, leading to incompatibilities. JAX is
an alternative implementation of Numpy on top of XLA, a domain-specific
compiler to optimize linear algebra for machine learning. It can JIT compile
any Numpy-based code, targeting CPUs, GPUs, and TPUs. It takes advantage of
the runtime knowledge of actual data types and sizes used, to specialize the
optimization. JAX can be extended, through the underlying XLA, with user-
defined functions and types. However, there is, at the time of writing, no
officially supported extension interface and suggested implementations use
undocumented APIs that are unsupported and subject to change. We find that
both PyPy and JAX are good alternative candidates for our stated goals, but
that Numba provides a better platform: where JAX supports GPUs well and where
PyPy supports extension with user-defined functions and types, Numba does
both. We will continue to track the developments of these alternative
technologies, however, and may revisit our choice in the future if warranted.

##  7 Conclusion

Seamless integration between programming languages is crucial in many
scientific domains, including high-energy physics. Python and C++ are widely
used in this field, and integrating the two languages is important for
efficient and accessible workflows. The just-in-time compiler Numba is capable
of compiling Python code to native but does not support C++. To address this
issue, we introduce C++ into Numba through a new reflection interface on top
of cppyy, which is an automatic, runtime bindings generator based on Cling,
the C++ interpreter. This approach provides a fully automatic and transparent
process, including integration with other Numba extensions, with performance
on par with traditional methods. This technology allows HEP programmers to
develop and debug their codes in Python, utilizing C++ libraries, while
switching on the Numba JIT for selected performance-critical closures. We
demonstrate 2-20 times speedup when using Numba to accelerate cppyy through
our extension. We outlined further gains using the Clang-Repl component of
LLVM and the newly developed library InterOp. The preliminary results from
using new infrastructure showed 1.4 to 144 times faster handling of templated
code in two different scenarios in cppyy, which will indirectly improve the
Numba-accelerated Python.

##  8 Acknowledgement

This work is supported by the National Science Foundation under Grant
OAC-1931408 and under Cooperative Agreement OAC-1836650. This material is
based upon work supported by the U.S. Department of Energy, Office of Science,
Advanced Scientific Computing Research (ASCR), SC Program 31.

## References

## References

  * [1] Pivarski J, Escott C, Smith N, Hedges M, Proffitt M, Escott C, Nandi J, Rembser J, bfis, benkrikler, Gray L, Davis D, Rodrigues E, Schreiner H, Gizdov K, Nollde and Fackeldey P 2021 scikit-hep/awkward-0.x: 0.15.5 URL https://doi.org/10.5281/zenodo.4520562

* [2] Pivarski J, Das P, Burr C, Smirnov D, Feickert M, Gal T, Kreczko L, Smith N, Biederbeck N, Shadura O, Proffitt M, benkrikler, Dembinski H, Schreiner H, Rembser J, R M, Gu C, Edoardo, Rodrigues E, JMSchoeffmann, Rübenach J, Koch L, Peresano M, Eich N, Turra R and bfis 2021 scikit-hep/uproot3: 3.14.4 URL https://doi.org/10.5281/zenodo.4537826
* [3] Gray L, Smith N, Tovar B, Novak A, Chakraborty J, Fackeldey P, Hartmann N, Watts G, Thain D, Stark G, BenGalewsky, Rübenach J, Fischer B, Taylor D, MoAly98, Kondratyev D, Gessinger P, Chen Y M E, Pata J, Woodard A, Potrebko A, Albert A, R M, slehti, Surma Z, Pedro K, dnoonan08 and Hennessee A 2023 Coffeateam/coffea: Release v0.7.21 URL https://doi.org/10.5281/zenodo.7733568
* [4] Brun R and Rademakers F 1997 Nucl. Instrum. Meth. A 389 81–86

* [5] Agostinelli S et al. (GEANT4) 2003 Nucl. Instrum. Meth. A 506 250–303
* [6] Lam S K, stuartarchibald, Pitrou A, Florisson M, Seibert S, Markall G, esc, Anderson T A, rjenc29, Leobas G, luk-f a, Collison M, Bourque J, Meurer A, Oliphant T E, Riasanovsky N, Kaustubh, Wang M, densmirn, njwhite, Pronovost E, Totoni E, Wieser E, Seefeld S, Grecco H, Peterson P, Virshup I, G M, Turner-Trauring I and Bourbeau J 2022 numba/numba: Version 0.56.4 URL https://doi.org/10.5281/zenodo.7289231
* [7] Lavrijsen W T L P and Dutta A 2016 Proceedings of the 6th Workshop on Python for High-Performance and Scientific Computing PyHPC ’16 (IEEE Press) p 27–35 ISBN 9781509052202
* [8] Vassilev V, Canal P, Naumann A, Moneta L and Russo P 2012 vol 396 (IOP Publishing) p 052071 URL https://iopscience.iop.org/article/10.1088/1742-6596/396/5/052071/pdf
* [9] llvmlite - A lightweight LLVM python binding for writing JIT compilers URL https://llvmlite.readthedocs.io/en/latest/
* [10] PyPy - A fast, compliant alternative implementation of Python URL https://www.pypy.org/
* [11] Bradbury J, Frostig R, Hawkins P, Johnson M J, Leary C, Maclaurin D, Necula G, Paszke A, VanderPlas J, Wanderman-Milne S and Zhang Q 2018 JAX: composable transformations of Python+NumPy programs URL http://github.com/google/jax
* [12] Shed Skin - a restricted-Python-to-C++ compiler URL https://github.com/shedskin/shedskin
* [13] Selective Embedded Just-In-Time Specialization URL https://sejits.eecs.berkeley.edu/
* [14] Arkouda - Interactive Data Analytics at Supercomputing Scale a Python API powered by Chapel URL https://bears-r-us.github.io/arkouda/