# Semiclassical equations of motion for disordered conductors: extrinsic
interband velocity, corrected collision integral and spin-orbit torques

Rhonald Burgos Atencia  Facultad de Ingenierías, Departamento de Ciencias
Básicas, Universidad del Sinú, Cra.1w No.38-153, 4536534, Montería, Córdoba
230002, Colombia  Qian Niu  Department of Physics, The University of Texas at
Austin, Austin TX 78712  Dimitrie Culcer  School of Physics, The University of
New South Wales, Sydney 2052, Australia

(August 27, 2024)

###### Abstract

The semiclassical equations of motion are widely used to describe carrier
transport in conducting materials. Nevertheless, the substantial challenge of
incorporating disorder systematically into the semiclassical model persists,
leading to quantitative inaccuracies and occasionally erroneous predictions
for the expectation values of physical observables. To address this issue, in
the present work we provide a general prescription for reformulating the
semiclassical equations of motion for carriers in disordered conductors by
taking the quantum mechanical density matrix as the starting point. We focus
on the case when only external electric fields are present, without magnetic
fields, and the disorder potential is spin-independent. The density matrix
approach allows averaging over impurity configurations, and the trace of the
velocity operator with the disorder-averaged density matrix can be
reinterpreted as the semiclassical velocity weighted by the Boltzmann
distribution function. Through this rationale the well-known intrinsic group
and anomalous velocities are trivially recovered, while we demonstrate the
existence of an extrinsic interband velocity, namely a disorder correction to
the semiclassical velocity of Bloch electrons, mediated by the interband
matrix elements of the Berry connection. A similar correction is present in
the non-equilibrium expectation value of the spin operator, contributing to
spin-orbit torques. To obtain agreement with diagrammatic approaches the
scattering term in the Boltzmann equation must be corrected to first order in
the applied electric field, and the Boltzmann equation itself must be solved
up to sub-leading order in the disorder potential. Our prescription ensures
all vertex corrections present in diagrammatic treatments are taken into
account, and to illustrate this we discuss model cases in topological
insulators, including the anomalous Hall effect as well as spin-orbit torques.

##  I Introduction

Carrier transport in extended conductors is an inherently semiclassical
phenomenon, requiring an effective single-particle description as well as
averaging over real space and momentum space degrees of freedom. The
semiclassical modelSundaram and Niu (1999), which takes these ingredients as
natural building blocks, has been a staple of transport theory for the best
part of a centuryAshcroft and Mermin (1976). In recent years, aside from its
well-established application to the anomalous Hall effect and dynamics in
magnetic systems Jungwirth et al. (2002); Matsumoto and Murakami (2011); Cheng
et al. (2016); Gao and Xiao (2019); Zhang et al. (2019), it has frequently
been used to describe transport in electric and magnetic fields in systems
with non-trivial topological textures Shindou and Imura (2005); Dumitrescu et
al. (2012); Son and Spivak (2013); Song et al. (2015); Yang et al. (2015);
Lensky et al. (2015); Chang and Yang (2015); Gorbar et al. (2018); Misaki et
al. (2018); Araki (2018); Alexandradinata and Glazman (2018); Das and Agarwal
(2019); Nandy et al. (2019); Hou et al. (2019); Zhu et al. (2021); Yokoyama
(2021); Wang et al. (2021); Das and Agarwal (2021); Chaudhary et al. (2021);
Bhowal and Vignale (2021), including non-linear electromagnetic responses
Moore and Orenstein (2010); Sodemann and Fu (2015); Morimoto et al. (2016);
Zeng et al. (2019); Golub et al. (2020) and has recently found substantial
applications in computational approaches to non-equilibrium physics Wang et
al. (2006, 2007); Gradhand et al. (2012); He et al. (2012); Chen et al.
(2013); Bianco et al. (2014); Chen et al. (2014); Olsen and Souza (2015); Feng
et al. (2016); Dai et al. (2017); Martiny et al. (2019); Wuttke et al. (2019);
Du et al. (2020); He et al. (2020); He and Law (2021). In addition to its
profound physical insight, broad applicability, and relative simplicity, the
semiclassical method naturally accounts for topological effects, and enables a
clear identification of Fermi surface and Fermi sea contributions to
transportChang and Niu (2008); Culcer et al. (2017); Stedman et al. (2019);
Stedman and Woods (2020).

The central idea of the semiclassical model is the separation between the
dynamics of individual carriers and the carrier distribution. Carrier dynamics
between collisions are described by the semiclassical equations of motion,
which do not incorporate disorder, while collisions are taken into account
through the Boltzmann equation, and affect solely the distribution function,
inducing changes in the occupation of quantum statesAshcroft and Mermin
(1976). The semiclassical velocity, originally assumed to be simply the band
group velocity, is now known to incorporate a transversal anomalous component
linear in the driving electric field and proportional to the Berry curvature
${\bm{\Omega}}_{m}$ of a given band $m$. Sundaram and Niu (1999); Xiao et al.
(2010) Although written in terms of the curvature for a single band this
anomalous velocity includes inter-band coherence effects, and is associated
with band mixing by an electric field.Culcer et al. (2017) The anomalous
velocity lies at the heart of the quantum Hall effect and of the intrinsic
contribution to the anomalous Hall effect, together with its quantized
counterpart. In recent years, however, it has been realised that disorder
itself leads to band-mixing effects which are not captured by the Boltzmann
equation, and are challenging to include in the wave packet description, since
averaging over disorder configurations cannot be done at the level of the wave
function. It is well established that a naive application of the semiclassical
model to the anomalous and spin-Hall effects in disordered systems makes
inaccurate predictionsInoue et al. (2004, 2006). Indeed, the role of disorder
in the anomalous Hall effect,Karplus and Luttinger (1954); Luttinger (1958);
Smit (1958); Berger (1970); Nozières and Lewiner (1973) and its relationship
to semiclassical dynamics, remains an intensely researched topicWang et al.
(2007); Kovalev et al. (2010); Yang et al. (2011); Bianco et al. (2014); Ado
et al. (2015, 2016, 2017); Rauch et al. (2018); Keser et al. (2019); Shen et
al. (2020).

Nevertheless, since transport is fundamentally semiclassical, all transport-
related quantities must be expressible in semiclassical terms. The assumptions
behind wave packet dynamics and diagrammatic approaches are the same: external
fields are treated classically and are assumed to be slowly-varying in space,
a separation is made between scattering processes and the dynamics between
scattering events, and the calculation is performed in the regime
$\epsilon_{F}\tau/\hbar\gg 1$, where $\epsilon_{F}$ is the Fermi energy and
$\tau$ the momentum relaxation time. Recent work has investigated strategies
for incorporating the findings of diagrammatic linear response theory into
semiclassical dynamics. Xiao and Niu showed that agreement is obtained with
diagrammatic approaches if all semiclassical quantities are dressed by
disorder, the cost being that the introduction of a disorder-dressed Berry
curvature Xiao and Niu (2017). Sinitsyn et al introduced a spin-dependent
coordinate shift into the position operator and identified a side-jump
velocity.Sinitsyn et al. (2005, 2007); Sinitsyn (2007) These quantities,
however, are difficult to work with, since a spin-dependent position operator
introduces complications of its own.

In light of the above, in this paper we formulate the semiclassical equations
of motion for Bloch electrons so as to include disorder using the density
matrix formalism, Culcer et al. (2017) in the process making a connection to
Green’s functions approaches. A strong motivation for our work is the recent
surge in applications of the semiclassical model in computational studies of
transport in topological materials. Wang et al. (2006, 2007); Gradhand et al.
(2012); He et al. (2012); Chen et al. (2013); Bianco et al. (2014); Chen et
al. (2014); Olsen and Souza (2015); Feng et al. (2016); Dai et al. (2017);
Martiny et al. (2019); Wuttke et al. (2019); Du et al. (2020); He et al.
(2020); He and Law (2021) Our primary aim is to provide a straightforward
method to incorporate disorder into such computational strategies once a model
of disorder is chosen. Taking the density matrix as the starting pointVasko
and Raichev (2005) allows one to average over disorder configurations,
something that cannot be done using a wave function. We demonstrate that
disorder affects not only the state occupation but also the semiclassical
equations of motion, and that it generates a correction to the velocity that
accounts for band mixing mediated by the Berry connection and disorder. This
approach enables one to distinguish disorder effects on the distribution
function from disorder effects on carrier dynamics, yet entails a change in
one’s point of view so as to regard the semiclassical equations as describing
carrier propagation averaged over many disorder scattering events. The carrier
undergoes transitions between bands as it scatters, and its trajectory can be
determined by averaging over impurity configurations. Whereas the equation of
motion for the wave vector ${\bm{k}}$ follows trivially from the time
derivative of the momentum operator, our central result is the revised
semiclassical equation of motion for the position of a carrier in band $m$,
with dispersion $\epsilon_{m}$, propagating under the action of an electric
field ${\bm{E}}$ in the presence of disorder

$\begin{array}[]{rl}\displaystyle\dot{\bm{r}}_{m}=&\displaystyle\frac{1}{\hbar}\frac{\partial\epsilon_{m}}{\partial\bm{k}}+\frac{e}{\hbar}\,{\bm{E}}\times{\bm{\Omega}}_{m}+{\bm{\beta}}^{m}_{\bm{k}}.\end{array}$
(1)

We identify a new contribution to the velocity, which we term the extrinsic
interband velocity ${\bm{\beta}}^{m}_{\bm{k}}$, defined as:

${\bm{\beta}}^{m}_{\bm{k}}=\frac{1}{2\pi\hbar}\int^{\infty}_{-\infty}d\epsilon\langle[U,G^{A}_{0}(\epsilon)[U,\bm{\mathcal{R}}^{\prime}]G^{R}_{0}(\epsilon)]\rangle^{mm}_{\bm{k}},$
(2)

where $\bm{\mathcal{R}}^{\prime}$ represents the inter-band Berry connection
with matrix elements $\bm{\mathcal{R}}^{{}^{\prime}mn}_{{\bm{k}}}=i\langle
u^{m}_{\bm{k}}|\nabla_{\bm{k}}u^{n}_{\bm{k}}\rangle_{m\neq n}$, and
$|u^{m}_{\bm{k}}\rangle$ is the lattice-periodic part of the Bloch wave
function. The extrinsic interband velocity ${\bm{\beta}}^{m}_{\bm{k}}$ is
proportional to the disorder strength, which is typically quantified by the
impurity density $n_{i}$, scattering potential strength $u^{2}_{0}$, or
alternatively $1/\tau$, where $\tau$ is the characteristic scattering time.
Formally, ${\bm{\beta}}^{m}_{\bm{k}}$ is similar to the customary scattering
term in the Born approximation, except the distribution function is replaced
by the band off-diagonal elements of the Berry connection. Since only the band
off-diagonal elements of the Berry connection appear
${\bm{\beta}}^{m}_{\bm{k}}$ is by construction gauge covariant. In general
${\bm{\beta}}^{m}_{\bm{k}}$, being independent of applied fields, can be
thought of as a disorder-dependent correction to the semiclassical band
velocity, or a random inter-band walk on the Fermi surface. We find that
${\bm{\beta}}^{m}_{\bm{k}}$ is nonzero in systems in which time reversal
symmetry is broken by e.g. a magnetization. Whereas
${\bm{\beta}}^{m}_{\bm{k}}$ is similar to the side jump as defined in Ref.
[Sinitsyn et al., 2006a], unlike Ref. Sinitsyn et al., 2006a, the present
formalism does not employ coordinate shifts, so that the formal and physical
position operators coincide. Furthermore, unlike Ref. [Xiao and Niu, 2017],
the Berry curvature is the same as in the clean system, rather than being
dressed by disorder. Importantly, we show that the scattering term in the
Boltzmann equation, needed to determine the effective distribution function,
acquires a correction to first order in the electric field, which is
equivalent to a gradient expansion in the electrostatic potential. In
addition, the Boltzmann equation needs to be solved up to the sub-leading
order in the impurity strength, in order to incorporate processes customarily
termed skew scattering and side jump.Smit (1958); Berger (1970); Nagaosa et
al. (2010) The method we present here also enables us to calculate spin
densities using the semiclassical model and obtain accurate results for spin-
orbit torques. In order to accomplish this the bare spin expectation value
needs to be supplemented with an electric field contribution, and we find an
analogous quantity to ${\bm{\beta}}^{m}_{\bm{k}}$ in the spin expectation
value.

More generally, we present a prescription for mapping steady-state expectation
values onto the semiclassical model by expressing traces purely in terms of
the band diagonal elements of the density matrix. Since in linear response
theory all expectation values are traced back to the equilibrium density
matrix, which is band-diagonal, they can all be recast in terms of
semiclassical quantities. The band-diagonal elements of the density matrix
represent the Boltzmann distribution, which can be evaluated from the much
simpler Boltzmann equation. In fact we will argue briefly in the latter part
of this work that linear response theories can be thought of as a family tree
with its roots in the quantum Liouville equation: the Kubo approach is the
integral formulation, the quantum kinetic, or quantum Boltzmann, approach is
the integro-differential formulation, while the semiclassical model is an
offshoot of the latter, which arises as a result of an additional separation
between the carrier dynamics and distribution. Sharing a common origin, these
methods yield equivalent results, and, in particular, vertex corrections
present in diagrammatic approaches have straightforward equivalents in the
semiclassical language. The blueprint presented in this work can be used in
the future to incorporate electron-electron interactions into computational
approaches in a mean field picture.

The outline of this paper is as follows. In Sec. (II) we introduce the
Hamiltonian and the model of disorder. In Sec.(III) we review linear response
theory based on the density matrix and introduce the electric-field correction
to the collision term. Next in Sec.(IV) we outline the general methodology for
deriving the semiclassical equations of motion from the quantum kinetic
equation, and discuss also disorder effects on spin expectation values. In
Sec. (V) we analysed the relation among different linear response
methodologies commonly used to calculate transport coefficients. In Sec.(VI)
we discuss at length two model examples, the anomalous Hall effect and spin-
orbit torques in magnetic topological insulators. We end with a summary and
conclusions.

##  II Model Hamiltonian

We consider a Hamiltonian of the form:

$H=H_{0}+V(\bm{r})+U(\bm{r}),$ (3)

with $H_{0}$ the low energy effective band Hamiltonian, in principle assumed
to include the Zeeman interaction with an external magnetic field, $V(\bm{r})$
is the electrostatic potential, and $U=U(\bm{r})$ represents the disorder
scattering potential. We emphasize the Hamiltonian is always Hermitian. Non-
Hermitian systems were considered in Ref. [Silberstein et al., 2020].

We work in the crystal momentum representation
$|m,\bm{k}\rangle=e^{i\bm{k}\cdot\bm{r}}|u^{m}_{\bm{k}}\rangle$. The matrix
elements of a scalar disorder potential $U(\bm{r})$ are given by the equation

$U^{mm^{\prime}}_{\bm{k}\bm{k}^{\prime}}=\langle
u^{m}_{\bm{k}}|u^{m^{\prime}}_{\bm{k}^{\prime}}\rangle\mathcal{U}_{\bm{q}},$
(4)

where we have defined the Fourier transform of the spatial function in $d$
dimensions

$\mathcal{U}_{\bm{q}}=\int d^{d}\bm{r}U(\bm{r})e^{-i\bm{q}\cdot\bm{r}},$ (5)

with $\bm{q}=\bm{k}-\bm{k}^{\prime}$. The impurity average is defined by

$\langle
U^{mm^{\prime}}_{\bm{k}\bm{k}^{\prime}}U^{nn^{\prime}}_{\bm{k}^{\prime}\bm{k}}\rangle=\langle
u^{m}_{\bm{k}}|u^{m^{\prime}}_{\bm{k}^{\prime}}\rangle\langle
u^{n^{\prime}}_{\bm{k}^{\prime}}|u^{n}_{\bm{k}}\rangle\langle\mathcal{U}_{\bm{q}}\mathcal{U}_{-\bm{q}}\rangle,$
(6)

with

$\langle\mathcal{U}_{\bm{q}}\mathcal{U}_{-\bm{q}}\rangle=\int\int
d^{d}\bm{r}d^{d}\bm{r}^{\prime}\langle U(\bm{r})U(\bm{r}^{\prime})\rangle
e^{-i\bm{q}\cdot(\bm{r}-\bm{r}^{\prime})},$ (7)

where $\langle\cdots\rangle$ refers to an average over impurity
configurations. For concreteness we will use a model of disorder whose spatial
correlations function is defined as:

$\displaystyle\langle U(\bm{r})\rangle$ $\displaystyle=0$ (8)
$\displaystyle\langle U(\bm{r})U(\bm{r}^{\prime})\rangle$
$\displaystyle=u^{2}_{0}\delta(\bm{r}-\bm{r}^{\prime})$ (9)

Then, it follows that

$\langle
U^{mm^{\prime}}_{\bm{k}\bm{k}^{\prime}}U^{nn^{\prime}}_{\bm{k}^{\prime}\bm{k}}\rangle=u^{2}_{0}\langle
u^{m}_{\bm{k}}|u^{m^{\prime}}_{\bm{k}^{\prime}}\rangle\langle
u^{n^{\prime}}_{\bm{k}^{\prime}}|u^{n}_{\bm{k}}\rangle,$ (10)

where $u^{2}_{0}$ is a parameter that takes into account the strength of the
disorder potential.

##  III Quantum Kinetic Equation

In this section we give a brief presentation of the quantum kinetic equation
in a somewhat different language than that used in Ref. [Culcer et al., 2017].
We note that similar density-matrix based approaches have been used recently
to describe carrier dynamics in the semiclassical regimeStedman et al. (2019);
Stedman and Woods (2020). The starting point is the quantum Liouville equation
for the single particle density operator $\rho$, namely:

$\frac{\partial\rho}{\partial t}+\frac{i}{\hbar}[H,\rho]=0.$ (11)

For the sake of convenience, we introduce at this stage the free retarded
Green function

$G^{R}_{0}(t)=-i\theta(t)e^{-itH_{0}/\hbar},$ (12)

In the frequency domain

$\displaystyle G^{R}_{0}(\epsilon)$
$\displaystyle=-\frac{i}{\hbar}\int^{\infty}_{0}dte^{-iH_{0}t/\hbar}e^{i\epsilon
t/\hbar}e^{-\eta t},$ (13)

where we introduced the factor $e^{-\eta t}$ to ensure convergence. The
advanced Green function follows by Hermitian conjugation.

###  III.1 Kinetic equation in equilibrium

For the sake of simplicity, let us for now ignore the effect of the driving
electric field in the kinetic equation. Using a decomposition of the density
matrix as $\rho=\langle\rho\rangle+g_{0}$ in the quantum Liouville
equationCulcer et al. (2017), we get for the disorder averaged part in
equilibrium:

$\frac{\partial\langle\rho\rangle}{\partial
t}+\frac{i}{\hbar}[H_{0},\langle\rho\rangle]+\frac{i}{\hbar}\langle[U,g_{0}]\rangle=0.$
(14)

while for $g_{0}$ we get the equation:

$\frac{\partial g_{0}}{\partial
t}+\frac{i}{\hbar}[H_{0},g_{0}]+\frac{i}{\hbar}[U,g_{0}]-\frac{i}{\hbar}\langle[U,g_{0}]\rangle=-\frac{i}{\hbar}[U,\langle\rho\rangle].$
(15)

In order to solve the kinetic equation for $\langle\rho\rangle$, we first
solve Eq.(15) for $g_{0}$ and then we use it in Eq.(14). In the first Born
approximation Culcer et al. (2017) we neglect the last two terms on the left
hand side of eq.(15). We are left with

$\frac{\partial g_{0}}{\partial
t}+\frac{i}{\hbar}[H_{0},g_{0}]=-\frac{i}{\hbar}[U,\langle\rho\rangle].$ (16)

Solving for $g_{0}$

$\displaystyle g_{0}$
$\displaystyle=-\frac{i}{\hbar}\int^{\infty}_{0}dt^{\prime}[e^{-iH_{0}t^{\prime}/\hbar}Ue^{iH_{0}t^{\prime}/\hbar},\langle\rho(t)\rangle].$
(17)

In terms of Green’s functions $g_{0}$ can be expressed as

$\displaystyle g_{0}$ $\displaystyle=\frac{1}{2\pi
i}\int^{\infty}_{0}d\epsilon[G^{R}_{0}(\epsilon)UG^{A}_{0}(\epsilon),\langle\rho(t)\rangle].$
(18)

This solution is substituted into Eq. (14). We arrive at the equation

$\frac{\partial\langle\rho\rangle}{\partial
t}+\frac{i}{\hbar}[H_{0},\langle\rho\rangle]+J(\langle\rho\rangle)=0,$ (19)

with the collision integral $J(\langle\rho\rangle)$ defined as:

$J(\langle\rho\rangle)=\frac{i}{\hbar}\langle[U,g_{0}]\rangle.$ (20)

###  III.2 Adding an electric field

Let us now consider the effect of the driving electrostatic potential up to
linear order. For simplicity we take this potential to have the form
$V(\bm{r})=e\bm{E}\cdot\bm{r}$, implying a uniform electric field, which
corresponds to the overwhelming majority of experimental setups. The case of
inhomogeneous systems, including systems in inhomogeneous electric fields,
entails additional subtleties which we postpone for later
considerationMarrazzo and Resta (2017); Lapa and Hughes (2019); Tu et al.
(2020); Kozii et al. (2021). Adding an electric field to the Hamiltonian
implies a correction to the function $g$, which can then be written as
$g=g_{0}+g_{E}$, where $g_{0}$ was found in the previous section, and

$\frac{\partial g_{E}}{\partial
t}+\frac{i}{\hbar}[H_{0},g_{E}]=-\frac{i}{\hbar}[V,g_{0}].$ (21)

The notation $g_{E}$ reflects the fact that eventually it is the electric
field that appears in the final expressions, rather than the electrostatic
potential. For $g_{E}$ we find explicitly

$\displaystyle g_{E}$
$\displaystyle=-\frac{i}{\hbar}\int^{\infty}_{0}dt^{\prime\prime}e^{-iH_{0}t^{\prime\prime}/\hbar}[V,g_{0}(t-t^{\prime\prime})]e^{iH_{0}t^{\prime\prime}/\hbar}.$
(22)

The function $g_{E}$ is off-diagonal in the momentum as well as in the band
index. We solve Eq.(22) by introducing Markovian approximation which reads
$g_{0}(t-t^{\prime\prime})\approx g_{0}(t)$, whereupon in the frequency domain
we obtain

$\displaystyle g_{E}$ $\displaystyle=\frac{1}{2\pi
i}\int^{\infty}_{0}d\epsilon G^{R}_{0}(\epsilon)[V,g_{0}]G^{A}_{0}(\epsilon),$
(23)

and in the commutator, we should use Eq.(18) as a functional of the
equilibrium distribution function $f_{0}(\epsilon)$ to fulfill linear
response.

The kinetic equation for the disorder averaged density matrix
$\langle\rho\rangle$ is now modified to

$\frac{\partial\langle\rho\rangle}{\partial
t}+\frac{i}{\hbar}[H_{0},\langle\rho\rangle]+J_{0}(\langle\rho\rangle)=-\frac{i}{\hbar}[V,\langle\rho\rangle]-J_{E}(\langle\rho\rangle),$
(24)

with the collision integral $J_{E}(\langle\rho\rangle)$ defined as:

$\displaystyle J_{E}(\langle\rho\rangle)$
$\displaystyle=\frac{i}{\hbar}\langle[U,g_{E}]\rangle.$ (25)

As we show below, the electric field correction to the collision integral in
Eq. (25) with the off-diagonal density function as given in Eq. (23) will
provide results in agreement with previous calculations based on diagrammatic
perturbation theoryAdo et al. (2015); Ndiaye et al. (2017). We consider such
an agreement as a positive test of the Markovian approximation. We note also
that Eq. (23) was used in a different but equivalent form in a previous paper
Culcer et al. (2010) in order to calculate side jump effects in a system with
extrinsic spin-orbit coupling. In this paper, we will focus on systems with
intrinsic spin-obit coupling.

###  III.3 Kinetic equation and linear response

When Eq. (24) is expressed in the crystal momentum representation we obtain
the quantum kinetic equationSekine et al. (2017)

$\frac{\partial f_{\bm{k}}}{\partial
t}+\frac{i}{\hbar}\,[H_{0{\bm{k}}},f_{\bm{k}}]+J_{0}(f_{\bm{k}})=\frac{e{\bm{E}}}{\hbar}\cdot\frac{Df_{\bm{k}}}{D{\bm{k}}}-J_{E}(f_{\bm{k}}).$
(26)

We have written the matrix elements of $\langle\rho\rangle$ in this
representation as $f_{\bm{k}}$. We refer to $f_{\bm{k}}$ henceforth as the
density matrix, noting that it has matrix elements connecting different bands,
although the band index $n$ has not been written explicitly. The covariant
derivative $\frac{Df_{\bm{k}}}{D{\bm{k}}}=\frac{\partial
f_{\bm{k}}}{\partial{\bm{k}}}-i[\bm{\mathcal{R}}_{\bm{k}},f_{\bm{k}}]$.

To solve Eq. 26, the density matrix is separated into a band diagonal and a
band off-diagonal part, namely, we write
$f_{\bm{k}}=n_{{\bm{k}}}+S_{{\bm{k}}}$. The band diagonal term $n_{{\bm{k}}}$
represents the fraction of carriers in a specific band and is essentially the
solution of the ordinary Boltzmann equation, while $S_{{\bm{k}}}$ contains the
effect of inter-band coherence, or band mixing. All our effort in recovering
the semiclassical theory consists of eliminating $S_{{\bm{k}}}$. The effective
Boltzmann equation that we shall derive is simply what is obtained for
$n_{{\bm{k}}}$ once all references to $S_{{\bm{k}}}$ have been eliminated.
Fortunately, as we recapitulate below, the solution for $S_{{\bm{k}}}$ in an
electric field is relatively simple, making it straightforward to express
expectation values in terms of $n_{{\bm{k}}}$ alone.

The equilibrium density matrix is band diagonal, its elements represented by
the Fermi-Dirac distribution for each band $n_{FD}(\epsilon^{m}_{\bm{k}})$. In
an electric field one may expand to linear order
$f_{\bm{k}}^{mn}=n_{FD}(\epsilon^{m}_{\bm{k}})\delta_{mn}+f_{E{\bm{k}}}^{mn}$,
with corresponding expressions for $n_{E{\bm{k}}}$ and $S_{E{\bm{k}}}$. The
kinetic equation is split into two coupled equations for $n_{E{\bm{k}}}$ and
$S_{E{\bm{k}}}$, whose solution, based on an expansion in the small parameter
$\hbar/(\epsilon_{F}\tau)$, is explained in detail in Ref. [Culcer et al.,
2017]. It was shown that $n_{E{\bm{k}}}$ starts at order $-1$ in this small
parameter, since it is proportional to the scattering time $\tau$, while
$S_{E{\bm{k}}}$ starts at order $0$. Consequently, the sub-leading correction
to $n_{E{\bm{k}}}$, referred to as $n_{E{\bm{k}}}^{(0)}$, is also required.

To leading order in $\hbar/(\epsilon_{F}\tau)$, the diagonal part reads:

$[J_{0}(n_{E}^{(-1)})]^{m}_{{\bm{k}}}=\frac{e\bm{E}}{\hbar}\cdot\frac{\partial
n_{FD}(\epsilon^{m}_{\bm{k}})}{\partial\bm{k}},$ (27)

where the Born approximation collision integral is

$\displaystyle[J_{0}(n_{E})]^{m}_{{\bm{k}}}$
$\displaystyle=\frac{2\pi}{\hbar}\sum_{m^{\prime},\bm{k}^{\prime}}\langle
U^{mm^{\prime}}_{\bm{k}\bm{k}^{\prime}}U^{m^{\prime}m}_{\bm{k}^{\prime}\bm{k}}\rangle$
(28)
$\displaystyle\times\left(n^{m}_{E\bm{k}}-n^{m^{\prime}}_{E\bm{k}^{\prime}}\right)\delta(\epsilon^{m}_{\bm{k}}-\epsilon^{m^{\prime}}_{\bm{k}^{\prime}}).$

The solution of Eq. (27) is in general rather complicatedAllen (1978). For a
system with isotropic dispersion it reduces to the simple form

$n_{E{\bm{k}}}^{m(-1)}=\tau^{m}_{p}\frac{e\bm{E}}{\hbar}\cdot\frac{\partial\epsilon^{m}_{\bm{k}}}{\partial\bm{k}}\frac{\partial
n_{FD}(\epsilon^{m}_{\bm{k}})}{\partial\epsilon^{m}_{\bm{k}}},$ (29)

where the transport time $\tau^{m}_{p}$ is defined as:

$\frac{1}{\tau^{m}_{p}}=\frac{2\pi}{\hbar}\sum_{m^{\prime},\bm{k}^{\prime}}\langle
U^{mm^{\prime}}_{\bm{k}\bm{k}^{\prime}}U^{m^{\prime}m}_{\bm{k}^{\prime}\bm{k}}\rangle[1-\cos(\theta_{\bm{k}^{\prime}}-\theta_{\bm{k}})]\delta(\epsilon^{m}_{\bm{k}}-\epsilon^{m^{\prime}}_{\bm{k}^{\prime}}).$
(30)

The solution for $S^{(0)}_{E\bm{k}}$ takes the simple formCulcer et al. (2017)

$S^{(0)mm^{\prime}}_{E{\bm{k}}}=\frac{\hbar(D+D^{\prime})^{mm^{\prime}}_{E\bm{k}}}{i(\epsilon^{m}_{\bm{k}}-\epsilon^{m^{\prime}}_{\bm{k}}-i\eta)}$
(31)

with the intrinsic and anomalous driving terms, Culcer et al. (2017)

$\displaystyle D^{mm^{\prime}}_{E\bm{k}}$
$\displaystyle=\frac{ie}{\hbar}\bm{E}\cdot\bm{\mathcal{R}}^{mm^{\prime}}_{\bm{k}}[n_{FD}(\epsilon^{m}_{\bm{k}})-n_{FD}(\epsilon^{m^{\prime}}_{\bm{k}})]$
(32) $\displaystyle D^{\prime mm^{\prime}}_{E\bm{k}}$
$\displaystyle=-\frac{\pi}{\hbar}\sum_{m^{\prime\prime},\bm{k}^{\prime}}\langle
U^{mm^{\prime\prime}}_{\bm{k}\bm{k}^{\prime}}U^{m^{\prime\prime}m^{\prime}}_{\bm{k}^{\prime}\bm{k}}\rangle\left\\{(n_{E{\bm{k}}}^{m^{\prime}}-n_{E{\bm{k}^{\prime}}}^{m^{\prime\prime}})\delta(\epsilon_{m^{\prime}\bm{k}}-\epsilon_{m^{\prime\prime}{\bm{k}^{\prime}}})\right.+\left.(n_{E{\bm{k}}}^{m}-n_{E{\bm{k}^{\prime}}}^{m^{\prime\prime}})\delta(\epsilon_{m^{\prime\prime}{\bm{k}^{\prime}}}-\epsilon_{m{\bm{k}}})\right\\}.$
(33)

Since $S_{E\bm{k}}$ starts at zeroth order in the parameter
$\hbar/(\epsilon_{F}\tau)$, we also require the sub-leading term
$n^{(0)}_{E\bm{k}}$, which is found from the equation

$[J_{0}(n^{(0)}_{E})]^{m}_{\bm{k}}=-[J_{sk}(n^{(-1)}_{E})]^{m}_{\bm{k}}-[J_{E}(n_{FD})]^{m}_{\bm{k}},$
(34)

where the right hand side acts as the driving term, whose constituents will be
explained shortly. Solving this equation will yield two different
contributions to the sub-leading diagonal density matrix
$n^{(0)}_{E{\bm{k}}}$, which we write as $n^{(0)}_{E{\bm{k}}}=n^{(\rm
sk)}_{E{\bm{k}}}+n^{(\rm sj)}_{E{\bm{k}}}$. Although both $n^{(\rm
sj)}_{E{\bm{k}}}$ and $n^{(\rm sk)}_{E{\bm{k}}}$ are of zeroth order in
$\hbar/(\epsilon_{F}\tau)$, they are parametrically different with respect to
magnetisation and Fermi energy, as we will see later on. We can solve for
these two terms separately as follows.

The contribution $n^{(\rm sk)}_{E{\bm{k}}}$ stems from $D^{\prime}$ and is
associated with skew scattering in the semiclassical theory. It is solved in
an analogous manner to Eq. (27), namely, the driving term is found by
substituting Eq. (31) into a collision integral of the form of Eq. (28),
obtaining

$J_{0}[n^{(\rm sk)}_{E}]=-[J_{sk}(n^{(-1)}_{E})]^{m}_{\bm{k}},$ (35)

which can be solved for $n^{(\rm sk)}_{E{\bm{k}}}$ using the standard
techniques of Boltzmann theoryAllen (1978). The driving term in this equation
can be written explicitly as a function of the leading-order density matrix
$n^{(-1)}_{E{\bm{k}}}$ as

$\displaystyle[J_{sk}(n^{(-1)}_{E})]^{m}_{\bm{k}}=\frac{2\pi^{2}}{\hbar}\sum_{m^{\prime}m^{\prime\prime}n{\bm{k}}^{\prime}{\bm{k}}^{\prime\prime}}$
$\displaystyle{\rm Im}\left[\frac{\langle
U^{mm^{\prime\prime}}_{\bm{k}\bm{k}^{\prime}}U^{m^{\prime}m}_{\bm{k}^{\prime}\bm{k}}\rangle\langle
U^{m^{\prime\prime}n}_{\bm{k}^{\prime}\bm{k}^{\prime\prime}}U^{nm^{\prime}}_{\bm{k}^{\prime\prime}\bm{k}^{\prime}}\rangle}{(\epsilon^{m^{\prime\prime}}_{\bm{k}^{\prime}}-\epsilon^{m^{\prime}}_{\bm{k}^{\prime}})}\right]$
$\displaystyle\left\\{(n_{E{\bm{k}^{\prime}}}^{m^{\prime}(-1)}-n_{E{\bm{k}^{\prime\prime}}}^{n(-1)})\delta(\epsilon^{m^{\prime}}_{\bm{k}^{\prime}}-\epsilon^{n}_{\bm{k}^{\prime\prime}})+(n_{E{\bm{k}^{\prime}}}^{m^{\prime\prime}(-1)}-n_{E{\bm{k}^{\prime\prime}}}^{n(-1)})\delta(\epsilon^{n}_{\bm{k}^{\prime\prime}}-\epsilon^{m^{\prime\prime}}_{\bm{k}^{\prime}})\right\\}\delta(\epsilon^{m^{\prime\prime}}_{\bm{k}^{\prime}}-\epsilon^{m}_{\bm{k}})$
(36)
$\displaystyle-\frac{2\pi^{2}}{\hbar}\sum_{m^{\prime}m^{\prime\prime}n{\bm{k}}^{\prime}{\bm{k}}^{\prime\prime}}$
$\displaystyle{\rm Im}\left[\frac{\langle
U^{m^{\prime\prime}m^{\prime}}_{\bm{k}\bm{k}^{\prime}}U^{m^{\prime}m}_{\bm{k}^{\prime}\bm{k}}\rangle\langle
U^{mn}_{\bm{k}\bm{k}^{\prime\prime}}U^{nm^{\prime\prime}}_{\bm{k}^{\prime\prime}\bm{k}}\rangle}{(\epsilon^{m}_{\bm{k}}-\epsilon^{m^{\prime\prime}}_{\bm{k}})}\right]$
$\displaystyle\left\\{(n_{E{\bm{k}}}^{m^{\prime\prime}(-1)}-n_{E{\bm{k}^{\prime\prime}}}^{n(-1)})\delta(\epsilon^{m^{\prime\prime}}_{\bm{k}}-\epsilon^{n}_{\bm{k}^{\prime\prime}})+(n_{E{\bm{k}}}^{m(-1)}-n_{E{\bm{k}^{\prime\prime}}}^{n(-1)})\delta(\epsilon^{n}_{\bm{k}^{\prime\prime}}-\epsilon^{m}_{\bm{k}})\right\\}\delta(\epsilon^{m^{\prime}}_{\bm{k}^{\prime}}-\epsilon^{m^{\prime\prime}}_{\bm{k}}),$

and we recall that $n_{E{\bm{k}}}^{(-1)}$ was found in Eq. (29).

The second contribution to the driving term in Eq. (34) is due to the electric
field correction of the collision integral $J_{E}$ acting on the equilibrium
distribution function. Since this contribution is associated with side jump
scattering in the semiclassical theory, we will refer to it as $n^{(\rm
sj)}_{E{\bm{k}}}$. To determine $n^{(\rm sj)}_{E{\bm{k}}}$ we need to solve
the equation

$J_{0}[n^{(\rm sj)}_{E}]^{m}_{\bm{k}}=-[J_{E}(n_{FD})]^{m}_{\bm{k}}.$ (37)

In the crystal momentum representation the electric field correction to the
collision integral takes the form

$\displaystyle[J_{E}(n_{FD})]^{m}_{\bm{k}}$
$\displaystyle=\frac{2\pi}{\hbar}\frac{\partial
n_{FD}(\epsilon^{m}_{\bm{k}})}{\partial\epsilon^{m}_{\bm{k}}}e\bm{E}\cdot\sum_{\bm{k}^{\prime}}\langle
U^{mm}_{\bm{k}\bm{k}^{\prime}}U^{mm}_{\bm{k}^{\prime}\bm{k}}\rangle\left[\bm{\mathcal{R}}^{mm}_{\bm{k}^{\prime}}-\bm{\mathcal{R}}^{mm}_{\bm{k}}\right]\delta(\epsilon^{m}_{\bm{k}^{\prime}}-\epsilon^{m}_{\bm{k}})$
$\displaystyle+\frac{2\pi}{\hbar}\frac{\partial
n_{FD}(\epsilon^{m}_{\bm{k}})}{\partial\epsilon^{m}_{\bm{k}}}e\bm{E}\cdot\sum_{m^{\prime}\bm{k}^{\prime}}{\rm
Im}\left\\{\left\langle\left[\left(\nabla_{\bm{k}}+\nabla_{\bm{k}^{\prime}}\right)U^{mm^{\prime}}_{\bm{k}\bm{k}^{\prime}}\right]U^{m^{\prime}m}_{\bm{k}^{\prime}\bm{k}}\right\rangle\right\\}\delta(\epsilon^{m^{\prime}}_{\bm{k}^{\prime}}-\epsilon^{m}_{\bm{k}}),$
(38)

where the derivatives act only on $U^{mm^{\prime}}_{\bm{k}\bm{k}^{\prime}}$.
This equation should be compared with the side jump velocity calculated from a
coordinate shift introduced in Ref.[Sinitsyn et al., 2006a]. The balance
between two collision integrals, as stated in Eq. (37) provides the necessary
information to calculate a new subleading density function $n^{(\rm sj)}_{E}$
that in the semiclassical language Sinitsyn et al. (2006b) is interpreted as
an anomalous distribution due to coordinate shift of the scattered particle
after many collisions.

##  IV Recovering the semiclassical theory

In this section we decompose the kinetic equation into a part representing
carrier dynamics and a part representing the distribution, which is found from
a modified Boltzmann equation. Since the equation of motion for the carrier
wave vector, yielding $\hbar\dot{\bm{k}}=-e{\bm{E}}$, follows immediately from
the operator commutator $[{\bm{p}},V({\bm{r}})]$, the bulk of our effort is
devoted to finding the disorder-averaged velocity, which will yield the time
evolution of the carrier position $\dot{\bm{r}}_{n}$. The prescription for
recovering the semiclassical theory from the quantum kinetic equation proceeds
as follows:

  * •

Determine the velocity expectation value as the operator trace Tr
$(\dot{\bm{r}}f)$, where $\bm{\dot{r}}=\frac{i}{\hbar}[H,\bm{r}]$ represents
the matrix elements of the velocity operator. In the crystal momentum
representation these are given by the covariant derivative
$\dot{\bm{r}}=\frac{1}{\hbar}\frac{DH}{D{\bm{k}}}$.

  * •

Reduce the trace to a form in which only band-diagonal elements of the density
matrix appear. These will contain either the equilibrium Fermi-Dirac
distribution $n_{0{\bm{k}}}$, or the correction to the band-diagonal part
$n_{E{\bm{k}}}$, which we recall has three constituents:
$n_{E{\bm{k}}}=n_{E{\bm{k}}}^{(-1)}+n_{E{\bm{k}}}^{(sk)}+n_{E{\bm{k}}}^{(sj)}$.

  * •

The result follows a natural separation into a contribution associated with
the equation of motion $\dot{\bm{r}}_{n}$ and one associated with the
Boltzmann equation.

  * •

For the spin density, we follow similar steps, namely, we take the trace of
the spin operator in the Bloch basis with the averaged density matrix. We will
also find an extrinsic spin matrix element that accounts for spin rotations
during scattering events.

The Hamiltonian is $H=H_{0}+V({\bm{r}})+U({\bm{r}})$, and since the last two
terms commute with the position operator they do not contribute to the
velocity operator. The band Hamiltonian yields

$\displaystyle{\rm Tr}(\dot{\bm{r}}f)\rightarrow{\rm
Tr}\left\\{\frac{i}{\hbar}[H_{0},\bm{r}]f\right\\}$
$\displaystyle=\frac{1}{\hbar}\sum_{m^{\prime},m,\bm{k}}\left[\frac{\partial\epsilon^{m^{\prime}}_{\bm{k}}}{\partial\bm{k}}\delta_{m,m^{\prime}}+i(\epsilon^{m}_{\bm{k}}-\epsilon^{m^{\prime}}_{\bm{k}})\bm{\mathcal{R}}^{mm^{\prime}}_{\bm{k}}\right]f^{m^{\prime}m}_{\bm{k}},$
(39)
$\displaystyle=\sum_{m,\bm{k}}\bm{v}^{m}_{\bm{k}}[n_{E{\bm{k}}}^{(-1)m}+n^{(0)mm}_{\bm{E}\bm{k}}]+\frac{i}{\hbar}\sum_{m^{\prime},m,\bm{k}}(\epsilon^{m}_{\bm{k}}-\epsilon^{m^{\prime}}_{\bm{k}})\bm{\mathcal{R}}^{mm^{\prime}}_{\bm{k}}S^{(0)m^{\prime}m}_{E{\bm{k}}}.$
(40)

The Berry connection $\bm{\mathcal{R}}^{mm^{\prime}}_{\bm{k}}=i\langle
u^{m}_{\bm{k}}|\nabla_{\bm{k}}u^{m^{\prime}}_{\bm{k}}\rangle$. The first term
gives the usual group velocity
$\bm{v}^{m}_{\bm{k}}=\nabla_{\bm{k}}\epsilon^{m}_{\bm{k}}/\hbar$ which is
diagonal, while the second term gives a contribution due to band mixing and is
purely off diagonal. We will concentrate on the second factor or band mixing
velocity. The off-diagonal density matrix is composed of two terms: an
intrinsic one and an extrinsic one. Let us first consider the intrinsic one.
It is

$\displaystyle{\rm Tr}\\{\dot{\bm{r}}f\\}^{\rm int}$
$\displaystyle=-\sum_{m^{\prime},m,\bm{k}}\bm{\mathcal{R}}^{mm^{\prime}}_{\bm{k}}D^{m^{\prime}m}_{E\bm{k}}.$
(41)

After replacing the driving term $D^{m^{\prime}m}_{E\bm{k}}$ by exchanging
$m\rightarrow m^{\prime}$ in the first term, and summing over intermediate
states, the intrinsic contribution can be written as the average of the
transverse velocity

$\displaystyle{\rm Tr}\\{\dot{\bm{r}}f\\}^{\rm int}$
$\displaystyle=\frac{e}{\hbar}\sum_{m,\bm{k}}\bm{E}\times{\bm{\Omega}^{m}_{\bm{k}}}n_{FD}(\epsilon^{m}_{\bm{k}})$
(42)

with the Berry curvature

$\Omega^{m}_{k,z}=i\left[\left\langle\frac{\partial u^{m}_{\bm{k}}}{\partial
k_{x}}\Big{|}\frac{\partial u^{m}_{\bm{k}}}{\partial
k_{y}}\right\rangle-\left\langle\frac{\partial u^{m}_{\bm{k}}}{\partial
k_{y}}\Big{|}\frac{\partial u^{m}_{\bm{k}}}{\partial
k_{x}}\right\rangle\right].$ (43)

The extrinsic contribution reads

$\displaystyle{\rm Tr}\\{\dot{\bm{r}}f\\}^{\rm ext}$
$\displaystyle=\sum_{n,m,\bm{k}}\bm{\mathcal{R}}^{nm}_{\bm{k}}[J(f)]^{mn}_{\bm{k}}.$
(44)

After some algebra it can be written as

$\displaystyle{\rm Tr}\\{\dot{\bm{r}}f\\}^{\rm ext}=\langle\bm{\beta}\rangle$
$\displaystyle=\sum_{m,\bm{k}}n_{E{\bm{k}}}^{(-1)m}\bm{\beta}^{m}_{\bm{k}},$
(45)

with the function $\bm{\beta}^{m}_{\bm{k}}$ formally defined as

$\displaystyle\bm{\beta}^{m}_{\bm{k}}$
$\displaystyle=\frac{1}{2\pi\hbar}\int^{\infty}_{-\infty}d\epsilon\langle[U,G^{A}_{0}(\epsilon)[U,\bm{\mathcal{R}^{\prime}}]G^{R}_{0}(\epsilon)]\rangle^{mm}_{\bm{k}},$
(46)

where the prime in $\bm{\mathcal{R}^{\prime}}$ indicates that only the band
off-diagonal matrix elements of the Berry connection enter. In Eq. (45) we
have written directly the electric-field dependent correction to the
distribution function, since Eq. (44) makes it obvious that this contribution
vanishes when $f$ is replaced by the equilibrium distribution $n_{FD}$. This
is because, for scalar scattering as studied in this work, the equilibrium
distribution causes the entire collision integral to vanish. For computational
evaluations it will be useful to list the explicit equation for
$\bm{\beta}^{m}_{\bm{k}}$:

$\displaystyle\bm{\beta}^{m}_{\bm{k}}$
$\displaystyle=\frac{\pi}{\hbar}\sum_{n,m^{\prime},\bm{k}^{\prime}}\Big{\\{}\left[\bm{\mathcal{R}}^{{}^{\prime}mn}_{\bm{k}}\langle
U^{nm^{\prime}}_{\bm{k},\bm{k}^{\prime}}U^{m^{\prime}m}_{\bm{k}^{\prime},\bm{k}}\rangle+\langle
U^{mm^{\prime}}_{\bm{k},\bm{k}^{\prime}}U^{m^{\prime}n}_{\bm{k}^{\prime},\bm{k}}\rangle\bm{\mathcal{R}}^{{}^{\prime}nm}_{\bm{k}}\right]\delta(\epsilon^{m}_{\bm{k}}-\epsilon^{m^{\prime}}_{\bm{k}^{\prime}})$
(47) $\displaystyle-\left[\langle
U^{mn}_{\bm{k},\bm{k}^{\prime}}U^{m^{\prime}m}_{\bm{k}^{\prime},\bm{k}}\rangle\bm{\mathcal{R}}^{{}^{\prime}nm^{\prime}}_{\bm{k}^{\prime}}+\langle
U^{mm^{\prime}}_{\bm{k},\bm{k}^{\prime}}U^{nm}_{\bm{k}^{\prime},\bm{k}}\rangle\bm{\mathcal{R}}^{{}^{\prime}m^{\prime}n}_{\bm{k}^{\prime}}\right]\delta(\epsilon^{n}_{\bm{k}^{\prime}}-\epsilon^{m}_{\bm{k}})\Big{\\}}$

Note that $\bm{\beta}^{m}_{\bm{k}}$ is proportional to the disorder strength
quantified here by $u^{2}_{0}$, making it first order in
$\hbar/(\epsilon_{F}\tau)$. It represents a disorder-dependent correction to
the semiclassical band velocity, which is independent of the applied electric
field. Physically, $\bm{\beta}^{m}_{\bm{k}}$ represents the average value of
the random changes in the carrier velocity that occur every time the carrier
is scattered between bands. Since $\bm{\beta}^{m}_{\bm{k}}$ has units of
velocity and depends on the disorder potential we will refer to it as the
extrinsic inter-band velocity. Given that $\bm{\beta}^{m}_{\bm{k}}$ is
formally of first order in $\hbar/(\epsilon_{F}\tau)$, we are only interested
in its product with the leading term in the distribution function,
$n_{E{\bm{k}}}^{(-1)}$, so that its overall contribution to the current is
formally zeroth order in disorder. Moreover, with $n_{E{\bm{k}}}^{(-1)}$
representing a Fermi surface contribution, the net effect of
$\bm{\beta}^{m}_{\bm{k}}$ can be thought of as a random inter-band walk on the
Fermi surface. Interestingly, ${\bm{\beta}}_{n}$ has the same mathematical
form as the Born approximation scattering term $J_{0}$, except the band off-
diagonal elements of $\bm{\mathcal{R}}$ appear instead of $n_{\bm{k}}$. The
presence of only the band off-diagonal matrix elements of $\bm{\mathcal{R}}$
ensures $\bm{\beta}^{m}_{\bm{k}}$ is gauge covariant. In the examples we study
below we find that $\bm{\beta}^{m}_{\bm{k}}$ is nonzero in systems in which
time reversal symmetry is broken by e.g. a magnetization. It is similar to the
side jump appearing in Ref. [Sinitsyn et al., 2006a], although we stress that
our approach makes no reference to any coordinate shifts, and the formal
position operator is identical to the physical position operator.

Since all contributions to the current density are now expressed in terms of
the distribution function (the equilibrium as well as the leading and sub-
leading terms in an electric field), we are able to write the semiclassical
equations of motion as

$\displaystyle\dot{\bm{r}}_{m}$
$\displaystyle=\frac{1}{\hbar}\frac{\partial\epsilon^{m}_{\bm{k}}}{\partial\bm{k}}-\dot{\bm{k}}_{m}\times\bm{\Omega}^{m}_{\bm{k}}+\bm{\beta}^{m}_{\bm{k}}$
(48) $\displaystyle\hbar\dot{\bm{k}}_{m}$ $\displaystyle=-e\bm{E}.$ (49)

The distribution function is found from the Boltzmann equation, with the
caveat that we require both the leading and subleading order terms in the
disorder strength. The procedure is as follows. First the leading-order term
in the distribution function $n_{E{\bm{k}}}^{(-1)}$ is found from

$J_{0}[n_{E}^{(-1)}]^{m}_{\bm{k}}=\frac{e\bm{E}}{\hbar}\cdot\frac{\partial
n_{FD}(\epsilon^{m}_{\bm{k}})}{\partial\bm{k}},$ (50)

while the sub-leading correction $n^{(0)}_{E}$ is given by

$J_{0}[n^{(0)}_{E}]^{m}_{\bm{k}}=-[J_{E}(n_{FD})]^{m}_{\bm{k}}-J_{sk}[n_{E}^{(-1)}]^{m}_{\bm{k}},$
(51)

where the left hand side is the quantity to be found, and the right hand side
plays the role of a driving term. Finally, we are able to write the full
expectation value of the current in terms of semiclassical quantities

$\langle\bm{j}\rangle=(-e)\,\sum_{m{\bm{k}}}\bigg{\\{}\frac{1}{\hbar}\frac{\partial\epsilon^{m}_{\bm{k}}}{\partial\bm{k}}[n^{(-1)m}_{E{\bm{k}}}+n^{(sk)m}_{E{\bm{k}}}+n^{(sj)m}_{E{\bm{k}}}]+\bigg{(}\frac{e\bm{E}}{\hbar}\times\bm{\Omega}^{m}_{\bm{k}}\bigg{)}n_{FD}(\epsilon^{m}_{\bm{k}})+\bm{\beta}^{m}_{\bm{k}}n^{(-1)m}_{E{\bm{k}}}\bigg{\\}}$
(52)

Let us consider the expectation value of the spin operator in the presence of
an electric field

${\rm
Tr}\\{\bm{s}f\\}=\sum_{m,\bm{k}}\bm{s}^{mm}_{\bm{k}}n^{m}_{E{\bm{k}}}+\sum_{n,m,\bm{k}}\bm{s}^{nm}_{\bm{k}}S^{mn}_{E\bm{k}},$
(53)

where $\bm{s}^{nm}_{\bm{k}}$ represent the matrix elements of the spin
operator. Writing explicitly the off diagonal terms of the density matrix in
the average of the spin operator we can separate the intrinsic and extrinsic
contributions as

$\displaystyle\langle{\bm{s}}\rangle^{\rm int}$
$\displaystyle=\sum_{n,m,\bm{k}}\frac{\hbar{\bm{s}}_{nm}}{i(\epsilon^{m}_{\bm{k}}-\epsilon^{n}_{\bm{k}})}\,D^{mn}_{E\bm{k}}$
(54) $\displaystyle\langle{\bm{s}}\rangle^{\rm ext}$
$\displaystyle=\sum_{n,m,\bm{k}}\frac{\hbar\bm{s}^{nm}_{\bm{k}}}{i(\epsilon^{m}_{\bm{k}}-\epsilon^{n}_{\bm{k}})}\,D^{{}^{\prime}mn}_{E\bm{k}}.$
(55)

If we define the quantity

${\bm{\mathcal{T}}}^{nm}_{\bm{k}}=\frac{\hbar}{i(\epsilon^{n}_{\bm{k}}-\epsilon^{m}_{\bm{k}})}{\bm{s}}^{nm}_{\bm{k}},$
(56)

the intrinsic contribution can be rewritten as

$\displaystyle\langle{\bm{s}}\rangle^{\rm int}$
$\displaystyle=\frac{ie}{\hbar}\sum_{m,\bm{k}}[\bm{\mathcal{T}}_{\bm{k}},\bm{E}\cdot\bm{\mathcal{R}^{\prime}}_{\bm{k}}]^{m}n_{FD}(\varepsilon^{m}_{\bm{k}}),$

where we only have to take off-diagonal components inside the commutator. This
enables us to define an intrinsic spin expectation expectation value as

$\displaystyle{\bm{s}}^{\rm
int}_{\bm{k}}=\frac{ie}{\hbar}[\bm{\mathcal{T}}_{\bm{k}},\bm{E}\cdot\bm{\mathcal{R}^{\prime}}_{\bm{k}}]^{m}.$
(57)

The extrinsic part can be written as

$\displaystyle\langle\bm{s}\rangle^{\rm ext}$
$\displaystyle=\sum_{n,m,\bm{k}}\bm{\mathcal{T}}^{nm}_{\bm{k}}[J_{0}(n^{(-1)}_{E})]^{mn}_{\bm{k}}$
(58)

This is mathematically analogous to the expression for the extrinsic interband
velocity Eq. (44), and using a similar manipulation we can re-express it as

$\displaystyle\langle\bm{s}\rangle^{\rm ext}_{\rm bm}$
$\displaystyle=\sum_{m,\bm{k}}n^{(-1)m}_{E{\bm{k}}}{\bm{\gamma}}^{mm}_{\bm{k}}$
(59)

where $n^{(-1)}_{E{\bm{k}}}$ is the leading order distribution function and we
introduce the new extrinsic spin $\bm{\gamma}^{m}_{\bm{k}}$ given by

$\displaystyle{\bm{\gamma}}^{m}_{\bm{k}}$
$\displaystyle=\frac{1}{2\pi\hbar}\int^{\infty}_{-\infty}d\epsilon\langle[U,G^{A}_{0}(\epsilon)[U,\bm{\mathcal{T}}]G^{R}_{0}(\epsilon)]\rangle^{mm}_{\bm{k}}.$
(60)

In exact analogy with ${\bm{\beta}}^{m}_{\bm{k}}$, since the Born
approximation scattering term vanishes when the distribution function is
replaced by the Fermi-Dirac distribution, ${\bm{\gamma}}^{m}_{\bm{k}}$ also
vanishes in equilibrium in the presence of scalar scattering. This quantity
represents spin rotations during scattering events. Again, the quantity
$\bm{\mathcal{T}}$ only has inter-band matrix elements. Explicitly
$\bm{\gamma}^{m}_{\bm{k}}$ is evaluated as

$\displaystyle{\bm{\gamma}}^{mm}_{\bm{k}}$
$\displaystyle=\frac{\pi}{\hbar}\sum_{n,m^{\prime},\bm{k}^{\prime}}\Big{\\{}\left[\bm{\mathcal{T}}^{m,n}_{\bm{k}}\langle
U^{nm^{\prime}}_{\bm{k},\bm{k}^{\prime}}U^{m^{\prime}m}_{\bm{k}^{\prime},\bm{k}}\rangle+\langle
U^{mm^{\prime}}_{\bm{k},\bm{k}^{\prime}}U^{m^{\prime}n}_{\bm{k}^{\prime},\bm{k}}\rangle\bm{\mathcal{T}}^{n,m}_{\bm{k}}\right]\delta(\epsilon^{m}_{\bm{k}}-\epsilon^{m^{\prime}}_{\bm{k}^{\prime}})$
(61) $\displaystyle-\left[\langle
U^{mn}_{\bm{k},\bm{k}^{\prime}}U^{m^{\prime}m}_{\bm{k}^{\prime},\bm{k}}\rangle\bm{\mathcal{T}}^{n,m^{\prime}}_{\bm{k}^{\prime}}+\langle
U^{mm^{\prime}}_{\bm{k},\bm{k}^{\prime}}U^{nm}_{\bm{k}^{\prime},\bm{k}}\rangle\bm{\mathcal{T}}^{m^{\prime},n}_{\bm{k}^{\prime}}\right]\delta(\epsilon^{n}_{\bm{k}^{\prime}}-\epsilon^{m}_{\bm{k}})\Big{\\}}.$

In analogy to the above, we can write a modified intra-band spin matrix
element as

$\displaystyle{\bm{s}}^{mm}\rightarrow\bm{s}^{mm}_{\bm{k}}+\frac{ie}{\hbar}[\bm{\mathcal{T}}_{\bm{k}},\bm{E}\cdot\bm{\mathcal{R}^{\prime}}_{\bm{k}}]^{m}+\bm{\gamma}^{m}_{\bm{k}}.$
(62)

The first factor is the bare matrix element of the spin operator in band $m$,
while the remaining two represent non-equilibrium corrections, the first being
intrinsic and the second extrinsic. We note that $\bm{\gamma}_{\bm{k}}$ and
$\bm{\beta}_{\bm{k}}$ are mathematically very similar and both contribute at
the Fermi energy. The intrinsic contribution to the spin density from the
Fermi sea appears in the second term in Eq. (62). The full semiclassical
expression for the spin density is given by

$\langle\bm{s}\rangle=\,\sum_{m{\bm{k}}}\bigg{\\{}\bm{s}^{m}_{\bm{k}}[n_{FD}(\varepsilon^{m}_{\bm{k}})+n^{(-1)m}_{E{\bm{k}}}+n^{(sk)m}_{E{\bm{k}}}+n^{(sj)m}_{E{\bm{k}}}]+\frac{ie}{\hbar}[\bm{\mathcal{T}}_{\bm{k}},\bm{E}\cdot\bm{\mathcal{R}^{\prime}}_{\bm{k}}]^{m}n_{FD}(\epsilon^{m}_{\bm{k}})+\bm{\gamma}^{m}_{\bm{k}}n^{(-1)m}_{E{\bm{k}}}\bigg{\\}}$
(63)

##  V Linear response family tree

In closing the methodological discussion we remark briefly on the
relationships between the various linear response theories. The most common
strategy for solving Eq. (11) in an electric field is via Kubo linear response
theory, which is discussed in detail in many textbooksBruus and Flensberg
(2004), hence we only dwell upon its fundamental aspects. Briefly, the
Hamiltonian is decomposed as $H=(H_{0}+U)+V$, and the non-equilibrium part of
the density matrix is likewise singled out as $\rho=\rho_{0}+\rho_{E}$. Then
in linear response one can write

$\frac{\partial\rho_{E}}{\partial
t}+\frac{i}{\hbar}\,[H_{0}+U,\rho_{E}]=-\frac{i}{\hbar}\,[V,\rho_{0}].$ (64)

This equation is solved immediately to yield

$\rho_{E}=\frac{1}{2\pi
i}\int^{\infty}_{0}d\epsilon[G^{R}(\epsilon)VG^{A}(\epsilon),\rho_{0}],$ (65)

where $G^{R}$, $G^{A}$ are the retarded and advanced Green’s functions,
respectively, for the disordered system:

$\displaystyle G^{R}(\epsilon)$
$\displaystyle=-\frac{i}{\hbar}\int^{\infty}_{0}dte^{-i(H_{0}+U)t/\hbar}e^{i\epsilon
t/\hbar}e^{-\eta t}.$ (66)

We refrain from writing out the energy dependence in full. Note that the
Green’s functions have not been averaged over disorder configurations at this
stage. To obtain the customary Kubo formula one must trace over the velocity
operator and average over impurity configurations

$\langle\bm{j}\rangle=-\frac{e}{2\pi i}\int^{\infty}_{0}d\epsilon\,{\rm
tr}\,\langle{\bm{v}}[G^{R}VG^{A},\rho_{0}]\rangle.$ (67)

The procedure is standard, so we do not cover it here in detail, the purpose
of this description is illustrative. The important point to notice is that the
Kubo formula is the integral approach to solving the Liouville equation,
whereas the kinetic equation we follow in this work represents the
differential approach, or integro-differential in view of the complex
scattering term. The Keldysh theory follows a similar path, and although it
takes as its starting point a series of Green’s functions, its ultimate origin
lies in the quantum Liouville equation. The Keldysh theory is formally non-
local in time, although in the vast majority of practical applications the
non-locality is removed and the Keldysh Green’s function, which is analogous
to the density matrix employed in this work, depends only on the difference in
time variables. Thus, for the purposes of the present comparison, the quantum
Boltzmann equation derived in the Keldysh theory is indistinguishable from the
quantum kinetic equation derived from the density matrix. The Kubo formula,
Keldysh theory and quantum kinetic equation may be regarded as holistic
approaches, in which both the carrier dynamics and the carrier distribution
are accounted for in the density matrix (or Keldysh Green’s function), and the
net result is the expectation value of a physical observable. In contrast, the
semiclassical theory involves a separation between the carrier dynamics and
carrier distribution, which can help to build an intuitive picture of the
underlying physics. The relationship between the different approaches is
summarized in the family tree of Fig. 1. The most important observation in
this context is the common origin of all linear response theories, which
reinforces the expectation that they should all lead to the same results.

Figure 1: Family tree of linear response theories.

##  VI Applications

We now turn to applications of the theory, which are intended to illustrate
the way the extrinsic interband velocity, extrinsic spin terms, and additional
scattering terms in the Boltzmann equation appear in the explicit evaluations
of physical observables for model systems. In particular, we emphasize the
relationship between these various contributions and the analogous quantities
appearing in diagrammatic theories, which enables us to reconcile the
semiclassical and diagrammatic results. Our focus will be on topological
insulatorsCulcer (2012), where we discuss the anomalous Hall effect as well as
spin-orbit torques, and compare the semiclassical results with previous work.

###  VI.1 Anomalous Hall effect in topological insulators

In this section we calculate the anomalous Hall conductivity in topological
insulators. The anomalous Hall conductivity is basically expressed in terms of
four contribution: intrinsic contribution that takes into account the whole
Fermi sea of the system, extrinsic contribution due to the extrinsic velocity
$\bm{\beta}^{m}_{\bm{k}}$ at the Fermi surface, side jump like contribution at
the Fermi energy due to an electric field correction to the collision integral
and a skew scattering contribution.

The Hamiltonian that describes low energy excitations in the surface of 3D
topological insulators reads:

$H=\hbar v_{F}(k_{x}\sigma_{y}-k_{y}\sigma_{x})+M\sigma_{z},$ (68)

where $v_{F}$ is the effective Fermi velocity and $\sigma_{i}$ are Pauli
matrices, while $M$ is the magnetization. The eigenvalues are
$\epsilon^{\pm}_{\bm{k}}=\pm\sqrt{\hbar^{2}v^{2}_{F}k^{2}+M^{2}}$, where $\pm$
labels conduction/valence band and the eigenstates are

$|u^{\pm}_{\bm{k}}\rangle=\frac{1}{\sqrt{2}}\begin{pmatrix}e^{-i\theta_{\bm{k}}}\sqrt{1\pm\xi_{k}}\\\
\pm i\sqrt{1\mp\xi_{k}}\end{pmatrix},$ (69)

with a parameter $\xi_{k}=M/\lambda_{\bm{k}}$ with
$\lambda_{\bm{k}}=\sqrt{\hbar^{2}v^{2}_{F}k^{2}+M^{2}}$ and
$\theta_{\bm{k}}=\arctan(k_{y}/k_{x})$.

From the eigenstates we determine the Berry connection vector
$\bm{\mathcal{R}}^{mm^{\prime}}_{\bm{k}}=i\langle
u^{m}_{\bm{k}}|\nabla_{\bm{k}}u^{m^{\prime}}_{\bm{k}}\rangle$. We decompose
this vector into a diagonal and an off-diagonal contribution, namely,
$\bm{\mathcal{R}}_{\bm{k}}=\bm{\mathcal{R}}^{d}_{\bm{k}}+\bm{\mathcal{R}}^{od}_{\bm{k}}$,
with

$\displaystyle\bm{\mathcal{R}}^{d}_{\bm{k}}$
$\displaystyle=\frac{\hat{\bm{\theta}}}{2k}(\sigma_{0}+\xi_{\bm{k}}\sigma_{z})$
(70) $\displaystyle\bm{\mathcal{R}}^{od}_{\bm{k}}$ $\displaystyle=\frac{\hbar
v_{F}}{2\lambda_{\bm{k}}}(\sigma_{x}\hat{\bm{\theta}}-\xi_{\bm{k}}\sigma_{y}\hat{\bm{k}}).$
(71)

We have defined a unit vector parallel to momentum
$\hat{\bm{k}}=(\cos\theta_{\bm{k}},\sin\theta_{\bm{k}})$ and a unit vector
perpendicular to momentum
$\hat{\bm{\theta}}=(-\sin\theta_{\bm{k}},\cos\theta_{\bm{k}})$. Also, these
vectors are related as $\hat{\bm{\theta}}=\hat{\bm{z}}\times\hat{\bm{k}}$.

Let us start by calculating the intrinsic contribution to the Hall
conductivity. It reads:

$\sigma^{\rm
int}_{yx}=\frac{e^{2}}{\hbar}\sum_{m,\bm{k}}\Omega^{m}_{\bm{k},z}n_{FD}(\epsilon^{m}_{\bm{k}}).$
(72)

Explicit evaluation for a topological insulator gives the Berry curvature

$\Omega^{\pm}_{k,z}=\mp\frac{\xi_{k}(1-\xi^{2}_{k})}{2k^{2}}.$ (73)

After a straightforward integration we find the intrinsic contribution to the
Hall conductivity to be

$\displaystyle\sigma^{\rm int}_{yx}$
$\displaystyle=\frac{e^{2}}{4\pi\hbar}\frac{M}{\epsilon_{F}}.$ (74)

The extrinsic inter-band velocity makes the following contribution to the
anomalous Hall conductivity

$\displaystyle\sigma^{\rm ext}_{yx}$
$\displaystyle=-\frac{e}{E_{x}}\sum_{m,\bm{k}}n_{E{\bm{k}}}^{(-1)m}\beta^{m}_{\bm{k},y}.$
(75)

The leading order correction to the distribution function is
$n_{E{\bm{k}}}^{(-1)m}=-e\tau_{tr}\bm{E}\cdot\bm{v}^{m}_{\bm{k}}\delta(\epsilon^{m}_{\bm{k}}-\epsilon_{F})$
with the transport time given by the expression

$\displaystyle\frac{1}{\tau_{tr}}=\frac{1}{2\tau}\left(1+3\xi^{2}_{\bm{k}}\right),$
(76)

with the scattering time defined as $1/\tau=\pi
u^{2}_{0}\rho(\epsilon_{\bm{k}})/\hbar$ and the density of states
$\rho(\epsilon_{\bm{k}})=\lambda_{k}/2\pi\hbar^{2}v^{2}_{F}$. The diagonal
velocity is
$v_{x}=v_{F}\left(1-\xi^{2}_{\bm{k}}\right)^{1/2}\cos\theta_{\bm{k}}$ and the
extrinsic inter-band velocity in the conduction band takes the form

$\displaystyle\bm{\beta}$ $\displaystyle=\sigma_{0}\frac{1}{\tau}\frac{\hbar
v_{F}}{\lambda_{k}}\xi_{k}(1-\xi^{2}_{k})^{1/2}\hat{\bm{\theta}}.$ (77)

The extrinsic inter-band velocity is a transverse velocity since it is
proportional to the unit vector $\hat{\bm{\theta}}$. After explicit
integration we arrive at the expression

$\displaystyle\sigma^{\rm ext}_{yx}$
$\displaystyle=\frac{e^{2}}{2\pi\hbar}\frac{M}{\epsilon_{F}}\frac{\left(1-\xi^{2}_{F}\right)}{\left(1+3\xi^{2}_{F}\right)}.$
(78)

Notice that $n_{E{\bm{k}}}^{(-1)m}$ is inversely proportional to the impurity
density while the extrinsic velocity $\bm{\beta}$ is proportional to the
impurity density. As a result, the overall effect in Eq. (75) is independent
of disorder. The extrinsic inter-band velocity $\bm{\beta}^{m}_{\bm{k}}$
comprises the effect of disorder on carrier dynamics, namely, it can be
interpreted as an effective velocity of the electron after many collisions, in
contrast to the group velocity, which is a velocity between collisions. In
this sense our extrinsic velocity can be associated to the side jump velocity
encountered in previous semiclassical results Sinitsyn et al. (2006b) but with
the difference that $\bm{\beta}^{m}_{\bm{k}}$ is entirely due to band mixing
mediated by the off-diagonal components of the Berry connection vector and
that it is constructed from a collision integral without introducing any
quantity related to coordinate shift.

As we discussed earlier, there are also two contributions to the anomalous
Hall conductivity related to two different diagonal subleading density matrix
functions. Let us first calculate the anomalous Hall conductivity related to
the term $n^{(\rm sj)}_{E{\bm{k}}}$ in the distribution function. It reads

$\displaystyle\sigma^{\rm sj}_{yx}$
$\displaystyle=-\frac{e}{E_{x}}\sum_{\bm{k}}v^{++}_{\bm{k},y}n^{(\rm
sj)++}_{E{\bm{k}}}.$ (79)

The diagonal velocity reads

$\displaystyle v^{++}_{\bm{k},y}$
$\displaystyle=v_{F}(1-\xi^{2}_{\bm{k}})^{1/2}\sin\theta_{\bm{k}},$ (80)

while the correction to the distribution function is

$\displaystyle n_{E{\bm{k}}}^{(\rm sj)++}$ $\displaystyle=2\frac{\hbar
v_{F}}{\lambda_{k}}\xi_{k}\delta(\epsilon_{F}-\epsilon^{+}_{\bm{k}})\frac{(1-\xi^{2}_{k})^{1/2}}{(1+3\xi^{2}_{\bm{k}})}e\bm{E}\cdot\hat{\bm{\theta}}.$
(81)

Replacing all elements we find for the conductivity

$\displaystyle\sigma^{\rm sj}_{yx}$
$\displaystyle=\frac{e^{2}}{2\pi\hbar}\frac{M}{\epsilon_{F}}\frac{(1-\xi^{2}_{F})}{(1+3\xi^{2}_{F})}.$
(82)

This term doubles the contribution in Eq.(78) due to the extrinsic velocity
although in this case the effect of disorder is completely captured by
$n^{(\rm sj)}_{E{\bm{k}}}$.

In previous semiclassical studies an anomalous distribution function was
introduced as a result of a coordinate shift Sinitsyn et al. (2006b). In
contrast, we have derived $n^{(\rm sj)}_{E{\bm{k}}}$ from Eq. (37) without any
need for introducing a coordinate shift.

The contribution to the anomalous Hall conductivity related to the skew
scattering correction to the distribution function, $n^{(\rm
sk)}_{E{\bm{k}}}$, reads

$\displaystyle\sigma^{\rm sk}_{yx}$
$\displaystyle=-\frac{e}{E_{x}}\sum_{\bm{k}}\left[v^{++}_{\bm{k},y}n^{(\rm
sk)++}_{E{\bm{k}}}+v^{--}_{\bm{k},y}n^{(\rm sk)--}_{E{\bm{k}}}\right].$ (83)

The diagonal velocities are $v^{\pm\pm}_{\bm{k},y}=\pm
v_{F}(1-\xi^{2}_{\bm{k}})^{1/2}\sin\theta_{\bm{k}}$, and

$\displaystyle n_{E{\bm{k}}}^{(\rm sk)}$ $\displaystyle=\frac{3}{2}\frac{\hbar
v_{F}}{\lambda_{k}}\xi_{k}\frac{\left(1-\xi^{2}_{F}\right)^{3/2}}{\left(1+3\xi^{2}_{F}\right)^{2}}\delta(\epsilon_{F}-\epsilon^{+}_{\bm{k}})e\bm{E}\cdot\hat{\bm{\theta}}\sigma_{z}.$
(84)

After integration we obtain

$\displaystyle\sigma^{\rm sk}_{yx}$
$\displaystyle=\frac{e^{2}}{2\pi\hbar}\frac{M}{\epsilon_{F}}\frac{3(1-\xi^{2}_{F})^{2}}{2(1+3\xi^{2}_{F})^{2}}.$
(85)

Adding all the contributions to the Hall conductivity we get the final
expression

$\displaystyle\sigma_{yx}$
$\displaystyle=\frac{4e^{2}}{2\pi\hbar}\frac{M}{\epsilon_{F}}\left[\frac{1+\xi^{2}_{F}}{(1+3\xi^{2}_{F})^{2}}\right],$
(86)

in exact agreement with previous results using the non-crossing approximation
and diagrammatic perturbation theory Sinitsyn et al. (2006b); Sinitsyn et al.
(2007); Sinitsyn (2007); Ado et al. (2015); Nagaosa et al. (2010).

###  VI.2 Spin density and spin-orbit torques in topological insulators

In this section we determine the spin density and spin-orbit torques in
topological insulators with an out-of-plane magnetization described by the
effective Hamiltonian Eq. (68). As for the conductivity, the spin density has
five contributions: a dominant contribution from the Fermi surface leading to
the Edelstein effectEdelstein (1990), an intrinsic contribution from the Fermi
sea, a contribution due to the extrinsic spin expectation value
$\bm{\gamma}^{m}_{\bm{k}}$ at the Fermi surface, a side-jump like contribution
at the Fermi energy due to the electric field correction to the collision
integral and a skew scattering contribution.

The leading order contribution to the spin density is

$\langle{\bm{s}}\rangle^{\rm
Edel}=\sum_{m,\bm{k}}{\bm{s}}^{mm}_{\bm{k}}n_{E{\bm{k}}}^{(-1)m}.$ (87)

Using
$n_{E{\bm{k}}}^{(-1)m}=-e\tau_{tr}{\bm{E}}\cdot\bm{v}^{m}_{\bm{k}}\delta(\epsilon^{m}_{\bm{k}}-\epsilon_{F})$
and

$\displaystyle s^{\pm\pm}_{\bm{k},x}$
$\displaystyle=\mp(1-\xi^{2}_{\bm{k}})^{1/2}\sin\theta_{\bm{k}}$ (88)
$\displaystyle s^{\pm\pm}_{\bm{k},y}$
$\displaystyle=\pm(1-\xi^{2}_{\bm{k}})^{1/2}\cos\theta_{\bm{k}},$ (89)

the Edelstein effect contribution to the spin density is

$\langle{\bm{s}}\rangle^{\rm
Edel}=e\tau{\bm{E}}\times\hat{\bm{z}}v_{F}\rho(\epsilon_{F})\frac{1-\xi^{2}_{F}}{1+3\xi^{2}_{F}}.$
(90)

The fraction of the spin density weighted by the intrinsic driving term reads:

$\displaystyle\langle{\bm{s}}\rangle^{\rm int}$
$\displaystyle=\frac{ie\bm{E}}{\hbar}\cdot\sum_{nm\bm{k}}(\bm{\mathcal{T}}^{mn}_{\bm{k}}\bm{\mathcal{R}}^{nm}_{\bm{k}}-\bm{\mathcal{R}}^{mn}_{\bm{k}}\bm{\mathcal{T}}^{nm}_{\bm{k}})n_{FD}(\varepsilon^{m}_{\bm{k}}).$
(91)

The dot product is between the electric field and the Berry connection vector.
With the Berry connection

$\displaystyle\mathcal{R}^{+-}_{\bm{k},x}$ $\displaystyle=-\frac{\hbar
v_{F}}{2\lambda_{\bm{k}}}(\sin\theta_{\bm{k}}-i\xi_{\bm{k}}\cos\theta_{\bm{k}})$
(92) $\displaystyle\mathcal{R}^{+-}_{\bm{k},y}$ $\displaystyle=\frac{\hbar
v_{F}}{2\lambda}(\cos\theta_{\bm{k}}+i\xi_{\bm{k}}\sin\theta_{\bm{k}})$ (93)

and the off-diagonal spin expectation values

$\displaystyle s^{-+}_{\bm{k},x}$
$\displaystyle=(\xi_{\bm{k}}\sin\theta+i\cos\theta)$ (94) $\displaystyle
s^{-+}_{\bm{k},y}$ $\displaystyle=-(\xi_{\bm{k}}\cos\theta-i\sin\theta)$ (95)

we get the intrinsic correction to the spin density

$\displaystyle\langle{\bm{s}}\rangle^{\rm int}$ $\displaystyle=\frac{eM\hbar
v_{F}\rho(\epsilon_{F})}{2\epsilon^{2}_{F}}\,{\bm{E}}.$ (96)

The extrinsic correction is defined as

$\displaystyle\langle{\bm{s}}\rangle^{\rm ext}$
$\displaystyle=\sum_{m,\bm{k}}n_{E{\bm{k}}}^{(-1)m}{\bm{\gamma}}^{m}_{\bm{k}}$
(97) $\displaystyle\bm{\gamma}^{++}_{\bm{k}}$
$\displaystyle=-\sigma_{0}\frac{1}{\tau}\frac{\hbar}{\lambda_{\bm{k}}}\xi_{k}(1-\xi^{2}_{k})^{1/2}\hat{\bm{k}},$
(98)

yielding

$\displaystyle\langle{\bm{s}}\rangle^{\rm ext}=\frac{e\hbar
v_{F}\rho(\epsilon_{F})M}{\epsilon^{2}_{F}}\frac{1-\xi^{2}_{F}}{1+3\xi^{2}_{F}}\,{\bm{E}}.$
(99)

This extrinsic contribution to the spin density is the counterpart of the
extrinsic inter-band velocity contribution in the anomalous Hall effect. The
extrinsic spin $\bm{\gamma}^{m}_{\bm{k}}$ as defined in Eq. (60) is an
interband coherence effect mediated by an effective off-diagonal spin operator
defined in Eq. (56).

The side-jump contribution is given by

$\displaystyle\langle s_{x}\rangle^{\rm sj}$
$\displaystyle=\sum_{\bm{k}}s^{++}_{\bm{k},x}n^{(\rm sj)++}_{E{\bm{k}}}.$
(100)

Using the spin expectations values above and the density function given by
Eq.(81) we find that the side-jump contribution to the spin density reads

$\displaystyle\langle{\bm{s}}\rangle^{\rm sj}$ $\displaystyle=\frac{e\hbar
v_{F}\rho(\epsilon_{F})M}{\epsilon^{2}_{F}}\frac{(1-\xi^{2}_{F})}{(1+3\xi^{2}_{F})}\,{\bm{E}}.$
(101)

This term doubles the contribution of Eq. (99). This is the counterpart of the
anomalous distribution function introduced in the semiclassical theory
Sinitsyn et al. (2006b) and also calculated in Eq. (82) for the anomalous Hall
conductivity.

The skew scattering contribution takes the form

$\displaystyle\langle{\bm{s}}\rangle^{\rm skew}$
$\displaystyle=\sum_{\bm{k}}\left[{\bm{s}}^{++}_{\bm{k},x}n^{(\rm
sk)++}_{E{\bm{k}}}+{\bm{s}}^{--}_{\bm{k}}n^{(\rm sk)--}_{E{\bm{k}}}\right].$
(102)

Using the diagonal spin expectations values and the distribution function
found in Eq.(84) we find for the skew scattering contribution to the spin
density

$\displaystyle\langle{\bm{s}}\rangle^{\rm sk}$
$\displaystyle=\frac{e\hbar\rho(\epsilon_{F})v_{F}\xi_{F}}{\lambda_{F}}\frac{3(1-\xi^{2}_{F})^{2}}{2(1+3\xi^{2}_{F})^{2}}\,{\bm{E}}.$
(103)

This is the counterpart of the skew scattering term in the anomalous Hall
conductivity calculated in Eq. (85).

Finally, the total spin density of the system then reads

$\displaystyle\langle\bm{s}\rangle$ $\displaystyle=-e\tau
v_{F}\rho(\epsilon_{F})\frac{(1-\xi^{2}_{F})}{\left(1+3\xi^{2}_{F}\right)}\hat{\bm{z}}\times\bm{E}$
(104) $\displaystyle+4e\hbar
v_{F}\rho(\epsilon_{F})\frac{M}{\epsilon^{2}_{F}}\frac{(1+\xi^{2}_{F})}{(1+3\xi^{2}_{F})^{2}}\,{\bm{E}}.$

The spin-orbit torque is defined as
$\bm{\tau}=(2M/\hbar)\bm{m}\times\langle\bm{s}\rangle$, where $\bm{m}$ is a
unit vector in the direction of magnetisation that we took here as
$\bm{m}=\hat{\bm{z}}$. Restoring the factor of $\hbar/2$ in the spin matrix
elements and substituting the density of states
$\rho(\epsilon_{F})=\epsilon_{F}/2\pi\hbar^{2}v^{2}_{F}$, the spin-orbit
torque is finally given by

$\displaystyle\bm{\tau}$
$\displaystyle=-\frac{e\tau\epsilon_{F}M\rho(\epsilon_{F})}{2\pi\hbar^{2}v_{F}}\frac{(1-\xi^{2}_{F})}{\left(1+3\xi^{2}_{F}\right)}\bm{m}\times(\hat{\bm{z}}\times\bm{E})$
(105) $\displaystyle+\frac{2e}{\pi\hbar
v_{F}}\frac{M^{2}}{\epsilon_{F}}\frac{(1+\xi^{2}_{F})}{(1+3\xi^{2}_{F})^{2}}{\bm{m}}\times{\bm{E}},$

in agreement with previous results. Ndiaye et al. (2017); Sakai and Kohno
(2014) A peculiarity of topological insulators is that the velocity operator
is directly related to the spin operator as
$\hat{\bm{v}}=-v_{F}\hat{\bm{z}}\times\bm{\sigma}$. Then the current density
is $\langle\bm{j}\rangle=-e\langle\hat{\bm{v}}\rangle$ what implies that the
anomalous Hall conductivity and the spin density are related by
$\sigma_{yx}=J_{y}/E_{x}=ev_{F}\langle s_{x}\rangle$. This fact can indeed be
verified from Eq. (86) and Eq. (104).

##  VII Conclusions and Outlook

We have demonstrated that the semiclassical dynamics of electrons in
disordered solids can be determined using linear response theory by taking the
density matrix and quantum Liouville equation as the starting point. This
results in a disorder-dependent correction to the semiclassical equation of
motion for the carrier position, which we have termed the extrinsic inter-band
velocity $\bm{\beta}^{m}_{\bm{k}}$, and which accounts for the effect of
disorder on the carrier velocity after many collisions. In analogy to the
extrinsic inter-band velocity, an extrinsic correction to the spin expectation
value $\bm{\gamma}^{m}_{\bm{k}}$ is also present, which accounts for spin
rotations during scattering events. This is accompanied by an intrinsic,
electric field dependent non-equilibrium correction to the spin expectation
value, which is analogous to the anomalous velocity present in the
semiclassical equations of motion. At the same time, the Boltzmann equation
must be solved up to sub-leading order in the disorder strength. The collision
integral in the Boltzmann equation includes an electric field dependent
correction which is analogous to the side-jump scattering term, as well as an
additional correction analogous to the skew scattering term in systems with
intrinsic spin-orbit interactions. We have applied this theory to describe the
anomalous Hall effect and spin-orbit torques in topological insulators,
obtaining exact agreement with quantum mechanical results using the Kubo
formula.

The general prescription we have formulated in this work can be
straightforwardly generalised to include extrinsic spin-orbit and magnetic
impurity scattering, as well as magnetic fields, which, however, require more
effort since the description relies on the Wigner function. Our prescription
paves the way towards a systematic semiclassical picture encompassing a host
of effects that conventionally lie beyond the purview of semiclassical theory:
disorder effects beyond the Born approximation, such as weak localization,
electron-electron interactions in the mean-field approximation, as well as
Kondo physics in magnetic systems.

Acknowledgments. DC is supported by the Australian Research Council Future
Fellowship FT190100062. We would like to thank Leonid Glazman, Aydin Keser and
Dmitry Efimkin for a series of educational discussions.

## References

  * Sundaram and Niu (1999) G. Sundaram and Q. Niu, Phys. Rev. B 59, 14915 (1999), URL https://link.aps.org/doi/10.1103/PhysRevB.59.14915.
  * Ashcroft and Mermin (1976) N. W. Ashcroft and N. D. Mermin, _Solid State Physics_ (Holt-Saunders, 1976).
  * Jungwirth et al. (2002) T. Jungwirth, Q. Niu, and A. H. MacDonald, Phys. Rev. Lett. 88, 207208 (2002), URL https://link.aps.org/doi/10.1103/PhysRevLett.88.207208.
  * Matsumoto and Murakami (2011) R. Matsumoto and S. Murakami, Phys. Rev. Lett. 106, 197202 (2011), URL https://link.aps.org/doi/10.1103/PhysRevLett.106.197202.
  * Cheng et al. (2016) R. Cheng, S. Okamoto, and D. Xiao, Phys. Rev. Lett. 117, 217202 (2016), URL https://link.aps.org/doi/10.1103/PhysRevLett.117.217202.
  * Gao and Xiao (2019) Y. Gao and D. Xiao, Phys. Rev. Lett. 122, 227402 (2019), URL https://link.aps.org/doi/10.1103/PhysRevLett.122.227402.
  * Zhang et al. (2019) X. Zhang, Y. Zhang, S. Okamoto, and D. Xiao, Phys. Rev. Lett. 123, 167202 (2019), URL https://link.aps.org/doi/10.1103/PhysRevLett.123.167202.
  * Shindou and Imura (2005) R. Shindou and K.-I. Imura, Nuclear Physics B 720, 399 (2005), ISSN 0550-3213, URL https://www.sciencedirect.com/science/article/pii/S0550321305004360.
  * Dumitrescu et al. (2012) E. Dumitrescu, C. Zhang, D. C. Marinescu, and S. Tewari, Phys. Rev. B 85, 245301 (2012), URL https://link.aps.org/doi/10.1103/PhysRevB.85.245301.
  * Son and Spivak (2013) D. T. Son and B. Z. Spivak, Phys. Rev. B 88, 104412 (2013), URL https://link.aps.org/doi/10.1103/PhysRevB.88.104412.
  * Song et al. (2015) J. C. W. Song, G. Refael, and P. A. Lee, Phys. Rev. B 92, 180204 (2015), URL https://link.aps.org/doi/10.1103/PhysRevB.92.180204.
  * Yang et al. (2015) S. A. Yang, H. Pan, and F. Zhang, Phys. Rev. Lett. 115, 156603 (2015), URL https://link.aps.org/doi/10.1103/PhysRevLett.115.156603.
  * Lensky et al. (2015) Y. D. Lensky, J. C. W. Song, P. Samutpraphoot, and L. S. Levitov, Phys. Rev. Lett. 114, 256601 (2015), URL https://link.aps.org/doi/10.1103/PhysRevLett.114.256601.
  * Chang and Yang (2015) M.-C. Chang and M.-F. Yang, Phys. Rev. B 91, 115203 (2015), URL https://link.aps.org/doi/10.1103/PhysRevB.91.115203.
  * Gorbar et al. (2018) E. V. Gorbar, V. A. Miransky, I. A. Shovkovy, and P. O. Sukhachov, Phys. Rev. B 98, 045203 (2018), URL https://link.aps.org/doi/10.1103/PhysRevB.98.045203.
  * Misaki et al. (2018) K. Misaki, S. Miyashita, and N. Nagaosa, Phys. Rev. B 97, 075122 (2018), URL https://link.aps.org/doi/10.1103/PhysRevB.97.075122.
  * Araki (2018) Y. Araki, Scientific Reports 8, 15236 (2018), URL https://doi.org/10.1038/s41598-018-33655-w.
  * Alexandradinata and Glazman (2018) A. Alexandradinata and L. Glazman, Phys. Rev. B 97, 144422 (2018), URL https://link.aps.org/doi/10.1103/PhysRevB.97.144422.
  * Das and Agarwal (2019) K. Das and A. Agarwal, Phys. Rev. B 99, 085405 (2019), URL https://link.aps.org/doi/10.1103/PhysRevB.99.085405.
  * Nandy et al. (2019) S. Nandy, A. Taraphder, and S. Tewari, Phys. Rev. B 100, 115139 (2019), URL https://link.aps.org/doi/10.1103/PhysRevB.100.115139.
  * Hou et al. (2019) Z. Hou, Y.-F. Zhou, N.-X. Yang, and Q.-F. Sun, Communications Physics 2, 86 (2019), URL https://doi.org/10.1038/s42005-019-0186-9.
  * Zhu et al. (2021) P. Zhu, T. L. Hughes, and A. Alexandradinata, Phys. Rev. B 103, 014417 (2021), URL https://link.aps.org/doi/10.1103/PhysRevB.103.014417.
  * Yokoyama (2021) T. Yokoyama, Sci Rep 11, 12065 (2021), ISSN 2045-2322 (Electronic); 2045-2322 (Linking).
  * Wang et al. (2021) Z. Wang, L. Dong, C. Xiao, and Q. Niu, Phys. Rev. Lett. 126, 187001 (2021), URL https://link.aps.org/doi/10.1103/PhysRevLett.126.187001.
  * Das and Agarwal (2021) K. Das and A. Agarwal, Phys. Rev. B 103, 125432 (2021), URL https://link.aps.org/doi/10.1103/PhysRevB.103.125432.
  * Chaudhary et al. (2021) S. Chaudhary, C. Knapp, and G. Refael, Phys. Rev. B 103, 165119 (2021), URL https://link.aps.org/doi/10.1103/PhysRevB.103.165119.
  * Bhowal and Vignale (2021) S. Bhowal and G. Vignale, Phys. Rev. B 103, 195309 (2021), URL https://link.aps.org/doi/10.1103/PhysRevB.103.195309.
  * Moore and Orenstein (2010) J. E. Moore and J. Orenstein, Phys. Rev. Lett. 105, 026805 (2010), URL https://link.aps.org/doi/10.1103/PhysRevLett.105.026805.
  * Sodemann and Fu (2015) I. Sodemann and L. Fu, Phys. Rev. Lett. 115, 216806 (2015), URL https://link.aps.org/doi/10.1103/PhysRevLett.115.216806.
  * Morimoto et al. (2016) T. Morimoto, S. Zhong, J. Orenstein, and J. E. Moore, Phys. Rev. B 94, 245121 (2016), URL https://link.aps.org/doi/10.1103/PhysRevB.94.245121.
  * Zeng et al. (2019) C. Zeng, S. Nandy, A. Taraphder, and S. Tewari, Phys. Rev. B 100, 245102 (2019), URL https://link.aps.org/doi/10.1103/PhysRevB.100.245102.
  * Golub et al. (2020) L. E. Golub, E. L. Ivchenko, and B. Spivak, Phys. Rev. B 102, 085202 (2020), URL https://link.aps.org/doi/10.1103/PhysRevB.102.085202.
  * Wang et al. (2006) X. Wang, J. R. Yates, I. Souza, and D. Vanderbilt, Phys. Rev. B 74, 195118 (2006), URL https://link.aps.org/doi/10.1103/PhysRevB.74.195118.
  * Wang et al. (2007) X. Wang, D. Vanderbilt, J. R. Yates, and I. Souza, Phys. Rev. B 76, 195109 (2007), URL https://link.aps.org/doi/10.1103/PhysRevB.76.195109.
  * Gradhand et al. (2012) M. Gradhand, D. V. Fedorov, F. Pientka, P. Zahn, I. Mertig, and B. L. Györffy, Journal of Physics: Condensed Matter 24, 213202 (2012), URL https://doi.org/10.1088/0953-8984/24/21/213202.
  * He et al. (2012) Y. He, J. Moore, and C. M. Varma, Phys. Rev. B 85, 155106 (2012), URL https://link.aps.org/doi/10.1103/PhysRevB.85.155106.
  * Chen et al. (2013) Y. Chen, D. L. Bergman, and A. A. Burkov, Phys. Rev. B 88, 125110 (2013), URL https://link.aps.org/doi/10.1103/PhysRevB.88.125110.
  * Bianco et al. (2014) R. Bianco, R. Resta, and I. Souza, Phys. Rev. B 90, 125153 (2014), URL https://link.aps.org/doi/10.1103/PhysRevB.90.125153.
  * Chen et al. (2014) H. Chen, Q. Niu, and A. H. MacDonald, Phys. Rev. Lett. 112, 017205 (2014), URL https://link.aps.org/doi/10.1103/PhysRevLett.112.017205.
  * Olsen and Souza (2015) T. Olsen and I. Souza, Phys. Rev. B 92, 125146 (2015), URL https://link.aps.org/doi/10.1103/PhysRevB.92.125146.
  * Feng et al. (2016) W. Feng, C.-C. Liu, G.-B. Liu, J.-J. Zhou, and Y. Yao, Computational Materials Science 112, 428 (2016), ISSN 0927-0256, computational Materials Science in China, URL https://www.sciencedirect.com/science/article/pii/S0927025615005893.
  * Dai et al. (2017) X. Dai, Z. Z. Du, and H.-Z. Lu, Phys. Rev. Lett. 119, 166601 (2017), URL https://link.aps.org/doi/10.1103/PhysRevLett.119.166601.
  * Martiny et al. (2019) J. H. J. Martiny, K. Kaasbjerg, and A.-P. Jauho, Phys. Rev. B 100, 155414 (2019), URL https://link.aps.org/doi/10.1103/PhysRevB.100.155414.
  * Wuttke et al. (2019) C. Wuttke, F. Caglieris, S. Sykora, F. Scaravaggi, A. U. B. Wolter, K. Manna, V. Süss, C. Shekhar, C. Felser, B. Büchner, et al., Phys. Rev. B 100, 085111 (2019), URL https://link.aps.org/doi/10.1103/PhysRevB.100.085111.
  * Du et al. (2020) S. Du, P. Tang, J. Li, Z. Lin, Y. Xu, W. Duan, and A. Rubio, Phys. Rev. Research 2, 022025 (2020), URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.022025.
  * He et al. (2020) W.-Y. He, D. Goldhaber-Gordon, and K. T. Law, Nature Communications 11, 1650 (2020), URL https://doi.org/10.1038/s41467-020-15473-9.
  * He and Law (2021) W.-Y. He and K. T. Law, Phys. Rev. Research 3, L032012 (2021), URL https://link.aps.org/doi/10.1103/PhysRevResearch.3.L032012.
  * Chang and Niu (2008) M.-C. Chang and Q. Niu, Journal of Physics: Condensed Matter 20, 193202 (2008), URL https://doi.org/10.1088/0953-8984/20/19/193202.
  * Culcer et al. (2017) D. Culcer, A. Sekine, and A. H. MacDonald, Phys. Rev. B 96, 035106 (2017), URL https://link.aps.org/doi/10.1103/PhysRevB.96.035106.
  * Stedman et al. (2019) T. Stedman, C. Timm, and L. M. Woods, New Journal of Physics 21, 103007 (2019), URL https://doi.org/10.1088/1367-2630/ab3f5f.
  * Stedman and Woods (2020) T. Stedman and L. M. Woods, Phys. Rev. Research 2, 033086 (2020), URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.033086.
  * Xiao et al. (2010) D. Xiao, M.-C. Chang, and Q. Niu, Rev. Mod. Phys. 82, 1959 (2010), URL https://link.aps.org/doi/10.1103/RevModPhys.82.1959.
  * Inoue et al. (2004) J. Inoue, G. E. W. Bauer, and L. W. Molenkamp, Phys. Rev. B 70, 041303 (2004).
  * Inoue et al. (2006) J.-i. Inoue, T. Kato, Y. Ishikawa, H. Itoh, G. E. W. Bauer, and L. W. Molenkamp, Phys. Rev. Lett. 97, 046604 (2006).
  * Karplus and Luttinger (1954) R. Karplus and J. Luttinger, Phys. Rev. 95, 1154 (1954).
  * Luttinger (1958) J. M. Luttinger, Phys. Rev. 112, 739 (1958).
  * Smit (1958) J. Smit, Physica 24, 39 (1958).
  * Berger (1970) L. Berger, Phys. Rev. B 2, 4559 (1970).
  * Nozières and Lewiner (1973) P. Nozières and C. Lewiner, J. Phys. (Paris) 34, 901 (1973).
  * Kovalev et al. (2010) A. A. Kovalev, J. Sinova, and Y. Tserkovnyak, Phys. Rev. Lett. 105, 036601 (2010), URL https://link.aps.org/doi/10.1103/PhysRevLett.105.036601.
  * Yang et al. (2011) S. A. Yang, H. Pan, Y. Yao, and Q. Niu, Phys. Rev. B 83, 125122 (2011), URL https://link.aps.org/doi/10.1103/PhysRevB.83.125122.
  * Ado et al. (2015) I. A. Ado, I. A. Dmitriev, P. M. Ostrovsky, and M. Titov, EPL (Europhysics Letters) 111, 37004 (2015), URL https://doi.org/10.1209/0295-5075/111/37004.
  * Ado et al. (2016) I. A. Ado, I. A. Dmitriev, P. M. Ostrovsky, and M. Titov, Phys. Rev. Lett. 117, 046601 (2016), URL https://link.aps.org/doi/10.1103/PhysRevLett.117.046601.
  * Ado et al. (2017) I. A. Ado, I. A. Dmitriev, P. M. Ostrovsky, and M. Titov, Phys. Rev. B 96, 235148 (2017), URL https://link.aps.org/doi/10.1103/PhysRevB.96.235148.
  * Rauch et al. (2018) T. c. v. Rauch, T. Olsen, D. Vanderbilt, and I. Souza, Phys. Rev. B 98, 115108 (2018), URL https://link.aps.org/doi/10.1103/PhysRevB.98.115108.
  * Keser et al. (2019) A. C. Keser, R. Raimondi, and D. Culcer, Phys. Rev. Lett. 123, 126603 (2019), URL https://link.aps.org/doi/10.1103/PhysRevLett.123.126603.
  * Shen et al. (2020) J. Shen, Q. Yao, Q. Zeng, H. Sun, X. Xi, G. Wu, W. Wang, B. Shen, Q. Liu, and E. Liu, Phys. Rev. Lett. 125, 086602 (2020), URL https://link.aps.org/doi/10.1103/PhysRevLett.125.086602.
  * Xiao and Niu (2017) C. Xiao and Q. Niu, Phys. Rev. B 96, 045428 (2017), URL https://link.aps.org/doi/10.1103/PhysRevB.96.045428.
  * Sinitsyn et al. (2005) N. A. Sinitsyn, Q. Niu, J. Sinova, and K. Nomura, Phys. Rev. B 72, 045346 (2005), URL https://link.aps.org/doi/10.1103/PhysRevB.72.045346.
  * Sinitsyn et al. (2007) N. A. Sinitsyn, A. H. MacDonald, T. Jungwirth, V. K. Dugaev, and J. Sinova, Phys. Rev. B 75, 045315 (2007), URL https://link.aps.org/doi/10.1103/PhysRevB.75.045315.
  * Sinitsyn (2007) N. A. Sinitsyn, Journal of Physics: Condensed Matter 20, 023201 (2007), URL https://doi.org/10.1088/0953-8984/20/02/023201.
  * Vasko and Raichev (2005) F. T. Vasko and O. E. Raichev, _Quantum Kinetic Theory and Applications: Electrons, Photons, Phonons_ (Springer-Verlag New York, 2005), 1st ed.
  * Sinitsyn et al. (2006a) N. A. Sinitsyn, Q. Niu, and A. H. MacDonald, Phys. Rev. B 73, 075318 (2006a), URL https://link.aps.org/doi/10.1103/PhysRevB.73.075318.
  * Nagaosa et al. (2010) N. Nagaosa, J. Sinova, S. Onoda, A. H. MacDonald, and N. P. Ong, Rev. Mod. Phys. 82, 1539 (2010), URL https://link.aps.org/doi/10.1103/RevModPhys.82.1539.
  * Silberstein et al. (2020) N. Silberstein, J. Behrends, M. Goldstein, and R. Ilan, Phys. Rev. B 102, 245147 (2020), URL https://link.aps.org/doi/10.1103/PhysRevB.102.245147.
  * Marrazzo and Resta (2017) A. Marrazzo and R. Resta, Phys. Rev. B 95, 121114 (2017), URL https://link.aps.org/doi/10.1103/PhysRevB.95.121114.
  * Lapa and Hughes (2019) M. F. Lapa and T. L. Hughes, Phys. Rev. B 99, 121111 (2019), URL https://link.aps.org/doi/10.1103/PhysRevB.99.121111.
  * Tu et al. (2020) M. W.-Y. Tu, C. Li, and W. Yao, Phys. Rev. B 102, 045423 (2020), URL https://link.aps.org/doi/10.1103/PhysRevB.102.045423.
  * Kozii et al. (2021) V. Kozii, A. Avdoshkin, S. Zhong, and J. E. Moore, Phys. Rev. Lett. 126, 156602 (2021), URL https://link.aps.org/doi/10.1103/PhysRevLett.126.156602.
  * Ndiaye et al. (2017) P. B. Ndiaye, C. A. Akosa, M. H. Fischer, A. Vaezi, E.-A. Kim, and A. Manchon, Phys. Rev. B 96, 014408 (2017), URL https://link.aps.org/doi/10.1103/PhysRevB.96.014408.
  * Culcer et al. (2010) D. Culcer, E. M. Hankiewicz, G. Vignale, and R. Winkler, Phys. Rev. B 81, 125332 (2010), URL https://link.aps.org/doi/10.1103/PhysRevB.81.125332.
  * Sekine et al. (2017) A. Sekine, D. Culcer, and A. H. MacDonald, Phys. Rev. B 96, 235134 (2017), URL https://link.aps.org/doi/10.1103/PhysRevB.96.235134.
  * Allen (1978) P. B. Allen, Phys. Rev. B 17, 3725 (1978), URL https://link.aps.org/doi/10.1103/PhysRevB.17.3725.
  * Sinitsyn et al. (2006b) N. A. Sinitsyn, J. E. Hill, H. Min, J. Sinova, and A. H. MacDonald, Phys. Rev. Lett. 97, 106804 (2006b), URL https://link.aps.org/doi/10.1103/PhysRevLett.97.106804.
  * Bruus and Flensberg (2004) H. Bruus and K. Flensberg, _Many-body quantum theory in condensed matter physics - an introduction_ (Oxford University Press, United States, 2004).
  * Culcer (2012) D. Culcer, Physica E: Low-dimensional Systems and Nanostructures 44, 860 (2012), ISSN 1386-9477, sI:Topological Insulators, URL https://www.sciencedirect.com/science/article/pii/S1386947711003985.
  * Edelstein (1990) V. Edelstein, Solid State Communications 73, 233 (1990), ISSN 0038-1098, URL https://www.sciencedirect.com/science/article/pii/003810989090963C.
  * Sakai and Kohno (2014) A. Sakai and H. Kohno, Phys. Rev. B 89, 165307 (2014), URL https://link.aps.org/doi/10.1103/PhysRevB.89.165307.