# Theory of nonlinear optical response

Amin Maleki Sheikhabadi  Department of Physics, Shahid Beheshti University,
Evin, Tehran 1983969411, Iran  Zahra Bagheri  School of Particles and
Accelerators, Institute for Research in Fundamental Sciences (IPM), P.O. Box
19395-5531, Tehran, Iran  Ali Sadeghi  Department of Physics, Shahid Beheshti
University, Evin, Tehran 1983969411, Iran

###### Abstract

We present a general formalism for investigating the second-order optical
response of solids to an electric field in weakly disordered crystals with
arbitrarily complicated band structures based on density-matrix equations of
motion, on a Born approximation treatment of disorder, and on an expansion in
scattering rate to leading non-trivial order. One of the principal aims of our
work is to enable extensive transport theory applications that accounts fully
for the interplay between electric-field-induced interband and intraband
coherence, and Bloch-state scattering. The quasiparticle bands are treated in
a completely general manner that allows for arbitrary forms of the intrinsic
spin-orbit coupling (SOC) and could be extended to the extrinsic SOC.
According to the previous results, in the presence of the disorder potential,
the interband response in conductors in addition to an intrinsic contribution
due to the entire Fermi sea that captures, among other effects, the Berry
curvature contribution to wave-packet dynamics includes an anomalous
contribution caused by scattering that is sensitive to the presence of the
Fermi surface. To demonstrate the rich physics captured by our theory, the
relaxation time matrix for different strength order is considered and at the
same time we explicitly solve for some electric-field response properties of
simple disordered Rashba model that are known to be dominated by interband
coherence contributions. The expressions we present are amenable for numerical
calculations, and we demonstrate this by performing a full band-structure
calculation of the interband contribution, even in metals.

††preprint: APS/123-QED

##  I Introduction

The study of steady state response of a solids to external electric fields is
subject of intensive theoretical and experimental investigations Refs. [1, 2,
3], but there are no reasons why this effect would not consider in a more
complex form like in metals [4, 5, 6, 7]. Originally the effect was first
observed in the very intensive laser light in 60s allowed irradiated systems
to be taken out of equilibrium that currently investigates to many electrical
applications [8]. The optical response of solids is fundamental theory of the
electric polarizations based on polarization suseptibilities [9],
characterizing not only linear response but also nonlinear of the materials to
the electromagnetic field [1]. Roughly speaking, although the susceptibilities
determine all the optical response of solids, but recently we have tarted to
probe this information generate by all the optical response of metals. An
example of these applications is to design of novel photovoltaic cells [10,
11, 7]. The reason for not examining this effect in metallic systems is that
the steady-state response of a metal or a doped semiconductor to external
magnetic field is complicated in the band structure [12, 13] and also the
essential roles are played by host crystal due to disorder or thermal
excitation of lattice degree-of-freedom [14]. Currently, fascinating ongoing
experiments start to observe this effect in the surface state of topological
insulators and Weyl semimetals [11, 4, 15, 6].

In this present work, we start with supposing that the Bloch state scattering
mechanisms that plays dominate role in the limit of the repopulation of states
close to the Fermi level in the direction of the electric field. The require
consideration to use Boltzmann equation is to assume that the dominate Bloch
state scattering mechanisms is weak [16] because of the finite size of the
system including many interfaces with strong spin-orbit splitting, including
metals with semiconductor giant SOC or insulators [17, 18, 19]. To do so, by
using quantum kinetic equation, we present a general formulation to describe
the Bloch state density matrix. The diagonal and off-diagonal terms in density
matrix are responsible for charge and spin density to spatially uniform and
constant field which appears at the leading sub-dominate order in a weak
scattering expansion. To the second order response, the diagonal term in
density matrix which is responsible for longitudinal response dominates and
captured by the relaxation time matrix. However the off-diagonal terms in
density matrix is arising from interband coherence terms that contains two
major parts, band structure and off-diagonal terms [12, 13], coming from the
entire Fermi surface including the Berry phase contributions and scattering
anomalous deriving term on the Fermi surface. Hence we analytically derive
both parts by using the Boltzmann equation and practically show that which
part under what conditions can be dominated. Generally we present a general
quantum kinetic theory of linear and nonlinear response to an electric field
applied to any model systems that requires an approach working for any band
structure with arbitrarily complicated band structures including the interband
and interaband coherence response and the Bloch-state repopulation responses
on an equal footing. Although this theory is a very powerful approach, it is
very useful to also show how the same results can be obtained independently by
using the different methods such as the SU(2) gauge theory formulation [19],
diagrammatic approach of the Kubo linear response theory [20, 21].

The outline of the rest of the paper follows. In the section [II], we consider
the simplest approach, known as “independent particle approximation” to
neglect the electron-electron interaction except insofar as it is included in
the calculation of the band structure of the solid. Clearly, quasiparticle
corrections and exciton effects associated with the creation of electron-hole
pairs, as well as local field corrections, should be considered in this
approach. In section [III], we develop a perturbation expansion in the case of
the Born approximation treatment of disorder, and on an expansion in
scattering rate to leading nontrivial order. The Bloch state density matrix of
a crystal to a spatially constant electric field assumes that the Wannier
representation Bloch Hamiltonian is present as a function of wave vector and
treats as short range impurity potential. Analysis shows that the verity of
different times in the approximations for disorder in crystal leads to the
density matrix response. Then our calculation focuses to the case when only
off-diagonal response of the Bloch state density matrix in metallic systems
present and in the simple case of disordered Rashba model our theory is
investigated [IV]. Finally, we state our conclusions in Section V.

##  II THE DIPOLE HAMILTONIAN

Full many-particle Hamiltonian in the presence of disorder potential is

$\displaystyle\hat{H}=\int{\boldmath{\bf\psi}}^{\dagger}({\bf
r},t)\left(\frac{{k}^{2}}{2m}+U_{b}({\bf r})+V({\bf r})\right){\bf\psi}({\bf
r},t),$ (1)

describing the Bloch electrons in disordered system. In Eq. 4, $U_{b}({\bf
r})$ is the periodic potential $U_{b}({\bf r}+a)$ with $a$ being any lattice
vector, and ${k}=-i\hbar\nabla$ the momentum operator. In the simplest case of
Gaussian white noise disorder potential, $V({\bf r})$ is assumed to have zero
average and Gaussian distribution given by $\langle V({\bf r})V({\bf
r}^{\prime})\rangle=(\hbar/(2\pi N_{0}\tau_{0}))\delta({\bf r}-{\bf
r}^{\prime})\rangle$, where $\langle\cdots\rangle$ denotes the angular average
over the momentum direction, $N_{0}=m/2\hbar^{2}\pi$, $n_{i}$ and $v_{0}$ are
the single-particle density of states per spin in the absence of SOC, the
impurity concentration and the scattering amplitude, respectively. $\tau_{0}$
is the elastic scattering time at the level of the Fermi Golden Rule. From now
on, we work with units such that $\hbar=c=1$, $e=-|e|$ for simplicity’s sake.

The homogeneous electric field is connected to the Hamiltonian with the
minimal substitution, ${{\bf k}}\rightarrow e{\bf A}(t)$, where ${\bf A}(t)$
is a vector potential that describes the (macroscopic Maxwell) electric field
as ${\bf E}=\dot{\bf A}(t)$. Hence we define a new field operator in the
Heisenberg picture as

$\displaystyle\psi({\bf r},t)\rightarrow\tilde{\psi}({\bf r},t)\equiv\psi({\bf
r},t)e^{-ie{\bf A}(t)\cdot{\bf r}}$ (2)

and it is easy to verify that in the Heisenberg picture it formally satisfies
the equations of motion

$\displaystyle i\frac{d\tilde{\psi}({\bf r},t)}{dt}=[\tilde{\psi}({\bf
r},t),H_{eff}(t)]$ (3)

then after the gauge transformation, the dipole Hamiltonian is written as

$\displaystyle{H}_{eff}(t)={H}_{0}+{H}_{\bf E}(t)+{H}_{imp}({\bf r}),$ (4)

where

$\displaystyle H_{0}=\int\tilde{\psi}^{\dagger}({\bf
r},t)\big{(}\frac{{k}^{2}}{2m}+U_{b}({\bf r})\big{)}\tilde{\psi}({\bf r},t)$
(5)

and the electric dependent part is given by

$\displaystyle{H}_{\bf E}=-e\int d{\bf r}{\bf\psi}^{\dagger}{\bf r}\cdot{\bf
E}{\bf\psi}\equiv-{\bf P}\cdot{\bf E},$ (6)

with ${\bf P}$ being the polarization operator. From Bloch’s theorem, the
eigenfunction of ${H}_{0}$ can be chosen of the form

$\displaystyle\psi_{n}({\bf r},t)=u_{n}({\bf k},{\bf r})e^{i{\bf k}\cdot{\bf
r}}$ (7)

with the band index $n$ and crystal momentum $\bf k$, where ${u}_{n}({\bf
k},{\bf r})={u}_{n}({\bf k},{\bf r}+{\bf a})$ is the period part of the Bloch
function. Hence in the Bloch states, the new field operation can be expanded
as

$\displaystyle\tilde{\psi}({\bf r},t)=\sum_{n{{\bf k}}}{\bf\psi}_{n}({\bf
r},t)a_{n}({\bf k})$ (8)

with $a_{n}^{\dagger}(\bf k)$ and $a_{n}(\bf k)$ are the creation and
annihilation operators, that obeys the anticommutation rules
$\\{a_{n}^{\dagger}({\bf k}),a_{n}({\bf k}^{\prime})\\}=\delta_{nm}\delta({\bf
k}-{\bf k}^{\prime})$. In this basis we have

$\displaystyle{H}_{0}=\sum_{n{\bf k}}\omega_{n}({\bf k})a_{n}^{\dagger}({\bf
k})a_{n}({\bf k})$ (9)

with $\omega_{n}$ the energy of band $n$. The model of disorder potential
arising from a serious randomly distributed impurities reads

$\displaystyle V({\bf r})=\sum_{\bf R}v({\bf r}-{\bf R})$ (10)

where $v$ is the potential of a single impurity placed at the origin of the
coordinate system, and $\bf R$ labels the coordinates. In the Bloch
representation we can easily show

$\displaystyle V^{mm^{\prime}}_{{\bf k}{\bf k}^{\prime}}=v^{mm^{\prime}}_{{\bf
k}{\bf k}^{\prime}}\sum_{\bf R}e^{-i({\bf k}-{\bf k}^{\prime})\cdot{\bf R}},$
(11)

where we have used $v^{mm^{\prime}}_{{\bf k}{\bf k}^{\prime}}\equiv\langle
m{\bf k}|V|m^{\prime}{\bf k}^{\prime}\rangle$ for the disorder potential
matrix elements in a band eigenstates representation. We refer to the
representation provided by the $|m{\bf k}\rangle$ basis as the eigenstate
representation and write $|m{\bf k}\rangle=e^{-i{\bf k}\cdot{\bf r}}|u_{m{\bf
k}}\rangle$. In the limit of Born approximation, it is straightforward to show
that $\langle V^{mm^{\prime}}_{{\bf k}{\bf
k}^{\prime}}V^{m^{\prime\prime}m^{\prime\prime\prime}}_{{\bf k}^{\prime}{\bf
k}}\rangle=n_{i}v^{mm^{\prime}}_{{\bf k}{\bf
k}^{\prime}}v^{m^{\prime\prime}m^{\prime\prime\prime}}_{{\bf k}^{\prime}{\bf
k}}$ [12]. Hence, in the Bloch representation we can write

$\displaystyle H_{imp}=\sum_{nm}\int d{\bf k}d{\bf k}^{\prime}v_{{\bf k}{\bf
k}^{\prime}}^{mn}a_{n}^{\dagger}({\bf k})a_{m}({\bf k}^{\prime}).$ (12)

To proceed with Eq. 6, we need to identify the electric polarization ${\bf
P}=e{\bf r}$ in the basis of Bloch fields as

$\displaystyle{\bf{P}}=e\sum_{nm{\bf k}{\bf k}^{\prime}}\langle n{{\bf
k}}|{\bf r}|m{\bf k}^{\prime}\rangle a_{n}^{\dagger}({\bf k})a_{m}({\bf
k}^{\prime}),$ (13)

where to proceed with above equation, one should derive an expression for the
matrix element

$\displaystyle\langle n{{\bf k}}|{\bf r}|m{\bf k}^{\prime}\rangle\equiv\int
d{\bf r}{\bf\psi}_{n}^{\dagger}({\bf k},{\bf r}){\bf r}{\bf\psi}_{m}({\bf
k}^{\prime},{\bf r}).$ (14)

For Eq.14, Yet Bount [22] first identified these matrix elements as

$\displaystyle\langle n{{\bf k}}|{\bf r}|m{\bf k}^{\prime}\rangle$
$\displaystyle=$ $\displaystyle\delta_{nm}\left[\delta({\bf k}-{\bf
k}^{\prime})\xi_{nn}({\bf k})+i\partial_{{\bf k}}\delta({\bf k}-{\bf
k}^{\prime})\right]$ (15) $\displaystyle+$
$\displaystyle(1-\delta_{nm})\delta({{\bf k}-{\bf k}^{\prime}})\xi_{nm}$

containing the correlation $\xi_{mn}$ as

$\displaystyle i\frac{\partial u_{n}({\bf k},{\bf r})}{\partial{\bf
k}}=\sum_{m}u_{m}({\bf k},{\bf r})\xi_{mn}$ (16)

and for nondegenerate points in the Brillouin zone

$\displaystyle\xi_{mn}=\frac{{\bf v}_{mn}({\bf k})}{i\omega_{mn}({\bf k})}.$
(17)

The velocity matrix elements ${\bf v}_{mn}({\bf k})$ is defined by

$\displaystyle\frac{1}{m}\int\psi_{n}^{\dagger}({\bf k},{\bf
x})(-i\nabla)\psi_{m}({\bf k}^{\prime},{\bf x})d{\bf x}={\bf v}_{mn}({\bf
k})\delta({\bf k}-{\bf k}^{\prime})$ (18)

and $m$ is the free electron mass. In order to make it convenient, we can use

$\displaystyle r_{nm}({\bf k})$ $\displaystyle\equiv$
$\displaystyle\xi_{nm}({\bf k})\qquad\qquad\textnormal{if}\quad n\neq m$ (19)
$\displaystyle\equiv$ $\displaystyle
0\quad\,\,\,\qquad\quad\qquad\textnormal{if}\quad n=m.$

According to Eq. (15), the polarization operator presented in Eq. 6 can be
separated into “interband” polarization $\hat{\bf P}_{i}$ for
$(1-\delta_{nm})$ component and the “intraband” polarization $\hat{\bf P}_{e}$
polarization for $\delta_{nm}$ component

$\displaystyle\hat{\bf P}$ $\displaystyle=$ $\displaystyle\hat{\bf
P}_{e}+\hat{\bf P}_{i}$ (20) $\displaystyle=$ $\displaystyle e\sum_{nm{\bf
k}}r_{nm}({\bf k})a_{n}^{\dagger}({\bf k})a_{n}({\bf k})$ $\displaystyle+$
$\displaystyle ie\sum_{n{\bf k}}a_{n}^{\dagger}({\bf k})a_{n,b}({\bf k})$

In Eq. 20, the interaband polarization dependents on the covariant derivative

$\displaystyle a_{n,b}({\bf k})\equiv\frac{\partial a_{n}({\bf
k})}{\partial{{k}^{b}}}-i(\xi_{mm}^{b}-\xi^{b}_{nn})a_{n}({\bf k}),$ (21)

where here and throughout, the Latin superscripts stand for Cartesian
components $b=x,y,z$. The main result of this paper is to generate quantum
kinetic equation that accounts for disorder and electric field. Throughout
this paper, we take disorder into account within the standard model of white-
noise disorder potential, implicitly assuming therefore that disorder is weak.
Although the extrinsic SOC can also be considered in the near future work.

##  III Current density operator

To establish some notations, in this part, we first consider the case when the
external electric and magnetic fields are absent and then, in the next
subsection, we will evaluate the Bloch equation in the more general when the
electric fields are present. We start with a single-particle density matrix
${\rho}=e{\bf\psi}^{\dagger}\psi$ described by the quantum Liouville equation

$\displaystyle\frac{\partial{\rho}}{\partial t}=\frac{1}{i}[H,{\rho}],$ (22)

where $H$ is the total Hamiltonian Eq. 4 of the complicated system containing
the band Hamiltonian, disorder potential and a perturbation $H_{E}$ due to the
external electric fields and all other interactions. To grasp proceed with the
scattering process, it is convenient to divide the density matrix into two
parts

$\displaystyle\rho=\langle\rho\rangle+{g},$ (23)

where $\langle\rho\rangle$ denotes the average over impurity configuration and
${g}$ is the reminder where its average over impurity configuration vanishes
$\langle{g}\rangle$. To keep the discussion as simple as possible, we confine
first to the case when the electric field is absence. The extension results to
$H_{E}\neq 0$ is straightforward. In Ref. [12, 13], a general quantum kinetic
theory of linear response to an external electric field applied to solids with
arbitrarily complicated band structures including the interband and intraband
coherence response and the formalism was based on density-matrix equations of
motion limited to the Born approximation treatment of disorder. To present a
consistent analysis we first recall the key steps of their derivations. The
starting point is the quantum Liouville Eq. 22 expands to the expressions for
$\langle\rho\rangle$ and $g$ as follows

$\displaystyle\frac{d\langle\rho\rangle}{dt}+{i}[H_{0},\langle\rho\rangle]+{i}\langle[v,g]\rangle=0,$
(24)
$\displaystyle\frac{d{g}}{dt}+{i}[H_{0},g]+{i}[v,\langle\rho\rangle]=-{i}[v,\hat{g}]+{i}\langle[v,{g}]\rangle.$
(25)

In Eq. 25 we can ignore both terms on the right hand side since they are
beyond the Born approximation, we can ignore them. Then from Eqs. 25 and Eq.
24, we can obtain [14]

$\displaystyle
g(t)={-i}\int_{0}^{\infty}dt^{\prime}[e^{-iH_{0}t^{\prime}}Ve^{iH_{0}t^{\prime}},\langle\rho(t)\rangle],$
(26)

where we have used the disorder-free expression for the time evolution of the
density matrix as follows

$\displaystyle\langle\rho(t-t^{\prime})\rangle=e^{-iH_{0}t^{\prime}}\langle\rho(t)\rangle
e^{iH_{0}t^{\prime}}.$ (27)

Hence by inserting Eq. 26 into Eq. 24 we can obtain the full kinetic equation
as

$\displaystyle\frac{d\langle\hat{\rho}\rangle}{dt}+\frac{i}{\hbar}[H_{0},\langle\hat{\rho}\rangle]+J(\langle\hat{\rho}\rangle)=0,$
(28)

where the electron-impurity collision integral is given by

$J(\langle\hat{\rho}\rangle)=\int_{0}^{\infty}dt^{\prime}\langle[V,[e^{-iH_{0}t^{\prime}}ve^{iH_{0}t^{\prime}},\langle\rho(t)\rangle]]\rangle.$
(29)

In this equation we notice that the dependence of $\langle\rho(t)\rangle$ to
the disorder potential is no longer explicitly and at this moment we can take
into account the translational symmetry recovered after impurity averaging by
working in momentum space where $J_{{\bf k}}(\langle\rho\rangle)$ is diagonal.
In the Bloch representation the collision integral at the momentum ${\bf k}$
is given by

$\displaystyle J_{{\bf k}}(\langle\rho\rangle)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}d{t^{\prime}}\sum$
$\displaystyle\bigl{\\{}\langle U^{mm^{\prime}}_{{\bf k}{\bf
k}^{\prime}}U^{m^{\prime}m^{\prime\prime}}_{{\bf k}^{\prime}{\bf k}}\rangle
e^{-i(\omega^{m^{\prime}}_{{\bf k}^{\prime}}-\omega^{m^{\prime\prime}}_{{\bf
k}})t^{\prime}}\langle\rho\rangle^{m^{\prime\prime}m^{\prime\prime\prime}}_{{\bf
k}}$ $\displaystyle-$ $\displaystyle\langle U^{mm^{\prime}}_{{\bf k}{\bf
k}^{\prime}}U^{m^{\prime\prime}m^{\prime\prime\prime}}_{{\bf k}^{\prime}{\bf
k}}\rangle\langle\rho\rangle^{m^{\prime}m^{\prime\prime}}_{{\bf
k}^{\prime}}e^{-i(\omega^{m^{\prime\prime}}_{{\bf
k}^{\prime}}-\omega^{m^{\prime\prime\prime}}_{{\bf k}})t^{\prime}}$
$\displaystyle-$ $\displaystyle\langle U^{mm^{\prime}}_{{\bf k}{\bf
k}^{\prime}}U^{m^{\prime\prime}m^{\prime\prime\prime}}_{{\bf k}^{\prime}{\bf
k}}\rangle e^{-i(\omega^{m}_{{\bf k}}-\omega^{m^{\prime}}_{{\bf
k}^{\prime}})t^{\prime}}\langle\rho\rangle^{m^{\prime}m^{\prime\prime}}_{{\bf
k}^{\prime}}$ $\displaystyle+$ $\displaystyle\langle
U^{m^{\prime}m^{\prime\prime}}_{{\bf k}{\bf
k}^{\prime}}U^{m^{\prime\prime}m^{\prime\prime\prime}}_{{\bf k}^{\prime}{\bf
k}}\rangle\langle\rho\rangle^{mm^{\prime}}_{{\bf
k}}e^{-i(\omega^{m^{\prime}}_{{\bf k}}-\omega^{m^{\prime\prime}}_{{\bf
k}^{\prime}})t^{\prime}}\bigr{\\}}.$

The integral over $t^{\prime}$ can be obtained by regularizing the time
integral by inserting a convergence factor $e^{\eta t^{\prime}}$ as follows

$\displaystyle\int_{0}^{\infty}$ $\displaystyle dt^{\prime}$ $\displaystyle
e^{-i(\omega^{m}_{{\bf k}}-\omega^{m^{\prime}}_{{\bf
k}^{\prime}}-i\eta)t^{\prime}}$ $\displaystyle\equiv$
$\displaystyle\mathcal{P}[\frac{1}{i(\omega^{m}_{{\bf
k}}-\omega^{m^{\prime}}_{{\bf k}^{\prime}})}]+{\pi}\delta(\omega^{m}_{{\bf
k}}-\omega^{m^{\prime}}_{{\bf k}^{\prime}})$

where $\delta$ part give us the conservation law, whereas the principle part
presents the disorder-induced level repulsion. Now we are able to back to Eq.
28 and solve the kinetic equation by decomposing the disorder averaged part of
the density matrix into two parts, $\langle\rho\rangle=n+S$, where $n$ is
diagonal in band index and $S$ is off-diagonal in band index. By using Eq.
(III), the diagonal part for the energy conservation term reads

$\displaystyle[J_{d}(\rho)]^{mm}_{{\bf k}}$ $\displaystyle=2\pi
n_{i}\sum_{m^{\prime}{\bf k}^{\prime}}u^{mm^{\prime}}_{{\bf k}{\bf
k}^{\prime}}u^{m^{\prime}m}_{{\bf k}^{\prime}{\bf k}}(\rho^{mm}_{{\bf
k}}-\rho^{m^{\prime}m^{\prime}}_{{\bf k}^{\prime}})\delta(\omega^{m}_{{\bf
k}}-\omega^{m^{\prime}}_{{\bf k}^{\prime}})$
$\displaystyle[J_{d}(\rho)]^{mm^{\prime\prime}}_{{\bf k}}$ $\displaystyle={\pi
n_{i}}\sum_{m^{\prime}{\bf k}^{\prime}}u^{mm^{\prime}}_{{\bf k}{\bf
k}^{\prime}}u^{m^{\prime}m^{\prime\prime}}_{{\bf k}^{\prime}{\bf k}}$ (33)
$\displaystyle\quad\quad\left[(\rho^{mm}_{{\bf
k}}-\rho^{m^{\prime}m^{\prime}}_{{\bf k}^{\prime}})\delta(\omega^{m}_{{\bf
k}}-\omega^{m^{\prime}}_{{\bf k}^{\prime}})\right.$
$\displaystyle\quad\qquad\quad\left.+(\rho^{m^{\prime\prime}m^{\prime\prime}}_{{\bf
k}}-\rho^{m^{\prime}m^{\prime}}_{{\bf
k}^{\prime}})\delta(\omega^{m^{\prime\prime}}_{{\bf
k}}-\omega^{m^{\prime}}_{{\bf k}^{\prime}})\right]$

where Eq. III presents the Fermi’s golden role when $m=m^{\prime\prime}$ and
Eq. III is raised when for $m\neq m^{\prime\prime}$. For principle part of the
collision integral term we get

$\displaystyle[J_{pp}(\rho)]^{mm^{\prime\prime}}_{{\bf k}}$ $\displaystyle=$
$\displaystyle{n_{i}}\sum u^{mm^{\prime}}_{{\bf k}{\bf
k}^{\prime}}u^{m^{\prime}m^{\prime\prime}}_{{\bf k}^{\prime}{\bf k}}$
$\displaystyle\left[(\rho^{m^{\prime\prime}m^{\prime\prime}}_{{\bf
k}}-\rho^{m^{\prime}m^{\prime}}_{{\bf
k}^{\prime}})\mathcal{P}(\frac{1}{\omega^{m^{\prime}}_{{\bf
k}^{\prime}}-\omega^{m^{\prime\prime}}_{{\bf k}}})\right.$ $\displaystyle+$
$\displaystyle\left.(\rho^{mm}_{{\bf k}}-\rho^{m^{\prime}m^{\prime}}_{{\bf
k}^{\prime}})\frac{1}{\omega^{m}_{{\bf k}}-\omega^{m^{\prime}}_{{\bf
k}^{\prime}}}\right]$

Notice that the diagonal elements $m=m^{\prime\prime}$ vanishes for the
principle part. In our analysis the principle terms are not considered and
only we take into account delta function of Eq. III since the real part has
been absorbed in just a shift of the Fermi level and the $\delta$-function
part becomes important near points in momentum space where different bands
touch.

###  III.1 Applying the electric field into the Kinetic equation

In this subsection, we add the external electric field to the system. As we
mention before, when $E\neq 0$, the term ${H}_{\bf E}$ adds to the total
Hamiltonian. We assume that the electric field is small like disorder term,
and with this assumption, we are able to use the perturbation theory.
According to Eq. 22, the quantum Liouville equation has the form

$\displaystyle\frac{d\rho}{dt}+{i}[H_{0},\rho]+J(\rho)=-{i}[H_{\bf E},\rho]$
(35)

Now we can solve the quantum kinetic perturbatively by applying a double
expansion on the impurity concentration $n_{i}$. Note that in this present
work we confine the results in the Born approximation and we neglect higher
order in impurities like the skew scattering effect and side-jump effect. To
do so, we derive the density matrix first to linear order to the electric
field and then we extend our results to second order. In linear response the
density matrix $\langle\rho\rangle$ can decompose into the equilibrium
$\langle\rho_{0}\rangle$, and the correlation to the equilibrium density
matrix to first order in the electric field $\langle\rho_{E}\rangle$. Hence
$\langle\rho\rangle=\langle\rho_{0}\rangle+\langle\rho_{E}\rangle$. By
linearizing the quantum kinetic equation, we obtain

$\displaystyle\frac{d\langle\rho_{E}\rangle}{dt}+{i}[H_{0},\langle\rho_{E}\rangle]+J(\langle\rho_{E}\rangle)=-{i}[H_{E},\langle\rho_{0}\rangle]-{i}[H_{E},\langle\rho_{E}\rangle]$
(36)

where
$\langle\rho_{E}\rangle=\lambda^{1}\langle\rho^{1}\rangle+\lambda^{2}\langle\rho^{2}\rangle+\cdots$
with expanding in the external electric field. Notice that in the response
theory, we can ignore the latest term on right hand side of Eq. 36 because it
contributes one order beyond of the electric field. Now we can start with
driving term on the right hand side of above equation by using the Bloch
eigenstate representation

$\displaystyle{-i}[H_{E},\langle\rho_{0}\rangle]$ $\displaystyle=$
$\displaystyle{-e}\sum_{b}E^{b}(t)\langle\rho_{0}\rangle_{mn,b}$
$\displaystyle+$
$\displaystyle{ie}\sum_{bp}E^{b}(t)\left(r^{b}_{mp}\langle\rho_{0}\rangle_{pn}-\langle\rho_{0}\rangle_{mp}r^{b}_{pn}\right).$

In the Bloch eigenstate representation, the ground state representation for
the density matrix was given by

$\displaystyle\langle\rho_{0}\rangle_{mn}=f_{0}(\omega^{m}_{{\bf
k}})\delta_{mn}\equiv n^{mm}_{0{\bf k}},$ (38)

where $f_{0}(\omega_{m}^{{\bf k}})$ presents the Fermi-Dirac distribution
function evaluated at energy $\omega_{m}^{{\bf k}}$. Therefore by using Eq.
(38) we can immediately derive

$\displaystyle{-i}[H_{E},\langle\rho_{0}\rangle_{nm}]$ $\displaystyle=$
$\displaystyle{-e{\bf E}}\cdot\left[\delta_{mn}\frac{\partial
f_{0}(\omega^{n}_{{\bf k}})}{\partial{{\bf k}}}\right.$
$\displaystyle\qquad-i\left.\sum_{p}{\bf r}_{nm}\left(f_{0}(\omega^{n}_{{\bf
k}})-f_{0}(\omega^{m}_{{\bf k}})\right)\right]$

The later term in square brackets denotes off-diagonal electric response when
$m\neq m^{\prime}$, and is responsible for the interband coherent
contributions to the electric response of solids. By using Eq. III.1 and
limiting to the Born approximation for the collision integral, we can derive
the density matrix for the diagonal terms and then for off-diagonal part is
straight forward. Hence by using Eq. 36 for the diagonal part of the density
matrix $m=m^{\prime}$, we have

$\displaystyle
i\frac{d\langle\rho^{1}\rangle_{mm}}{dt}-[H_{0},\langle\rho^{1}\rangle]_{mm}-\langle[v,\rho^{1}]\rangle_{mm}=[H_{E},\langle\rho^{0}\rangle]_{mm}$
(40)

where the second term on the left hand side of above equation is purely off-
diagonal in the band index and it can be neglected for diagonal part of the
density matrix. The later term can replace by Eq. III.1 where just the first
term on that equation is diagonal in matrix space. The collision term on the
later term of Eq. 40, can be seen as a matrix operator acting on the density
matrix. For the diagonal response we can define it according to

$\displaystyle{i}\langle[v,\rho^{1}]\rangle_{mm}\equiv
J_{d}(\rho^{1})^{mm}_{{\bf k}}=n_{i}\sum_{m^{\prime}{\bf
k}^{\prime}}\mathcal{K}^{mm^{\prime}}_{{\bf k}{\bf
k}^{\prime}}\rho^{1m^{\prime}m^{\prime}}_{{\bf k}^{\prime}}$ (41)

where $\mathcal{K}$ acts on the density matrix as a vector with components
labeled by band and wave vector. In comparison with Eq. III, we have

$\displaystyle\mathcal{K}^{mm^{\prime}}_{{\bf k}{\bf k}^{\prime}}$
$\displaystyle=$ $\displaystyle{2\pi n_{i}}\left[\delta_{{\bf k}{\bf
k}^{\prime}}\delta_{mm^{\prime}}\sum_{m^{\prime\prime},{\bf
k}^{\prime\prime}}u^{mm^{\prime\prime}}_{{\bf k}{\bf
k}^{\prime\prime}}u^{m^{\prime\prime}m}_{{\bf k}^{\prime\prime}{\bf
k}}\delta(\omega^{m}_{{\bf k}}-\omega^{m^{\prime}}_{{\bf k}^{\prime}})\right.$
(42) $\displaystyle\qquad\qquad\quad\quad-\left.u^{mm^{\prime}}_{{\bf k}{\bf
k}^{\prime}}u^{m^{\prime}m}_{{\bf k}^{\prime}{\bf k}}\delta(\omega^{m}_{{\bf
k}}-\omega^{m^{\prime}}_{{\bf k}^{\prime}})\right]$

where for the simplicity of the calculations from now we consider
$\langle\hat{\rho}^{1}\rangle=\hat{\rho}^{1}$. Ansatz for the diagonal matrix
in longer time limit,

$\rho_{{\bf k},{\bf k}}^{1mm}=\sum_{b\beta}\widetilde{\rho}^{1,b\beta}_{m{\bf
k},m{\bf k}}E^{b}_{\beta}e^{-i\omega_{\beta}t},$ (43)

where monochromatic electric field is defined by

$E^{b}=\sum_{\beta}E_{\beta}^{b}e^{-i\omega_{\beta}t}.$ (44)

By taking the derivation of Eq. 43, and inserting Eq. 41 and Eq. III.1 inside
Eq. 40, we can easily show that

$(-i\omega_{\beta})\hat{\rho}^{1,b\beta}_{m{\bf k},m{\bf
k}}+n_{i}\sum_{m^{\prime}{\bf k}^{\prime}}\mathcal{K}_{{\bf k}{\bf
k}^{\prime}}^{mm^{\prime}}\hat{\rho}^{1,b\beta}_{m^{\prime}{\bf
k}^{\prime},m^{\prime}{\bf k}^{\prime}}={e}\frac{\partial{f({\omega^{m}_{{\bf
k}}})}}{\partial{k}^{b}}$ (45)

By inserting Eq. 43 to Eq. 40, we read

$\hat{\rho}_{{\bf k},{\bf k}}^{1mm}={e}\sum_{b\beta}\sum_{m^{\prime}{\bf
k}^{\prime}}(\Gamma^{mm^{\prime}}_{{\bf k}{\bf
k}^{\prime}})^{-1}\frac{\partial{f({\omega_{{\bf
k}}^{m}})}}{\partial{k^{\prime}}^{b}}E^{b}_{\beta}e^{-i\omega_{\beta}t}$ (46)

where

$\displaystyle\Gamma^{mm^{\prime}}_{{\bf k}{\bf
k}^{\prime}}=-i\omega_{\beta}\delta({{\bf k}-{\bf
k}^{\prime}})\delta_{mm^{\prime}}+n_{i}\mathcal{K}^{mm^{\prime}}_{{\bf k}{\bf
k}^{\prime}},$ (47) $\displaystyle\frac{\partial{f({\omega_{{\bf
k}}^{m}})}}{\partial{k}^{b}}=\hat{k}_{b}v^{m}_{k}\delta(f({\omega_{{\bf
k}}^{m}})-\omega_{F}),$ (48) $\displaystyle\hat{{\bf v}}_{{\bf
k}}=\hat{k}\partial_{{\bf
k}}-\frac{\hat{\psi}}{k}\partial_{\psi},\qquad\text{with \ \
}\hat{\psi}=(-\sin\psi,\cos\psi),$ (49)

and $\omega_{F}$ is the Fermi level energy. It is clear that in the absence of
the disorder potential, the diagonal term in density matrix is zero and the
there is only a contribution to $\hat{\rho}_{{\bf k},{\bf k}}^{1mm^{\prime}}$
from the interband coupling when $m\neq m^{\prime}$. This is in agreement with
the results reported in Ref. [23, 24]. The collision term in Eq. 41 can be
seen as a matrix operator acting on the density matrix. By solving for
$\rho_{{\bf k}}^{(-1),mm}$ to leading order in $n_{i}$ according to Eq. 40, we
have

$\displaystyle\rho_{{\bf k}}^{(-1),mm}=n_{i}^{-1}e{\bf E}\cdot\sum_{m{\bf
k}^{\prime}}({\mathcal{K}^{mm^{\prime}}_{{\bf k}{\bf
k}^{\prime}}})^{-1}\cdot{\bf v}^{m^{\prime}}_{{\bf k}^{\prime}}\frac{\partial
f_{0}(\omega^{m^{\prime}}_{{\bf
k}^{\prime}})}{\partial\omega^{m^{\prime}}_{{\bf k}^{\prime}}}$ (50)

where we have used ${\bf v}^{m}_{{\bf k}}=\partial\omega_{{\bf
k}}^{m}/\partial{{\bf k}}$. In the simplest case of metals and insulators we
have

$\displaystyle\frac{\partial f_{0}(\omega_{{\bf k}}^{m})}{\partial\omega_{{\bf
k}}^{m}}$ $\displaystyle=$ $\displaystyle-\delta(\omega_{{\bf
k}}^{m}-\omega_{F})\quad\textnormal{for metals}$ (51) $\displaystyle=$
$\displaystyle 0\ \ \qquad\qquad\qquad\textnormal{for insulators}$ (52)

It is clear that the band-diagonal part of the density matrix Eq. 50 is a
contribution to balance the intraband driving term and scattering on the Fermi
surface and in the limit of very weak disorder scattering diverges. Now we can
derive the off-diagonal terms in the density matrix. Since
$[H_{0},\langle\rho_{0}\rangle]$ is purely off-diagonal, we can conclude that
the diagonal part on the density matrix is never divergence even in the
absence of disorder potential. Hence the leading response is independent of
disorder strength $n_{i}$ as

$\displaystyle\frac{\partial S^{1,mm^{\prime}}}{\partial{{\bf
k}}}+{i}[H_{0{\bf k}},S^{1}_{{\bf k}}]^{mm^{\prime}}=D^{mm^{\prime}}_{{\bf
k}{\bf k}^{\prime}}+D^{{}^{\prime}mm^{\prime}}_{{\bf k}{\bf k}^{\prime}}$ (53)

where the first term on the right hand side is intrinsic off-diagonal part
shown in later term of Eq. (III.1)

$\displaystyle D_{Ek}^{mm^{\prime}}={ie{\bf E}}\cdot\sum_{{\bf k}}{\bf
r}_{mm^{\prime}}\left(f_{0}(\omega_{{\bf k}}^{m})-f_{0}(\omega_{{\bf
k}}^{m^{\prime}})\right)$ (54)

and $D^{{}^{\prime}mm^{\prime}}_{{\bf k}{\bf k}^{\prime}}$ is responsible for
disorder scattering $D^{{}^{\prime}mm^{\prime}}_{{\bf k}{\bf
k}^{\prime}}=-J_{d}[\rho_{\bf k}^{1}]_{mm^{\prime}}$ which is non zero at
Fermi surface. By using the matrix element $J_{d}[\rho_{\bf
k}^{1}]_{mm^{\prime}}$ from Eq. III, $D^{{}^{\prime}mm^{\prime}}_{{\bf k}{\bf
k}^{\prime}}$ takes the form

$\displaystyle D^{{}^{\prime}mm^{\prime}}_{{\bf k}{\bf k}^{\prime}}$
$\displaystyle=$ $\displaystyle{\pi n_{i}}\sum_{m^{\prime}{\bf
k}^{\prime}}u^{mm^{\prime}}_{{\bf k}{\bf
k}^{\prime}}u^{m^{\prime}m^{\prime\prime}}_{{\bf k}^{\prime}{\bf k}}$
$\displaystyle\left[(\rho^{-1,mm}_{{\bf
k}}-\rho^{-1,m^{\prime}m^{\prime}}_{{\bf k}^{\prime}})\delta(\omega^{m}_{{\bf
k}}-\omega^{m^{\prime}}_{{\bf k}^{\prime}})\right.$ $\displaystyle+$
$\displaystyle\left.(\rho^{-1,m^{\prime\prime}m^{\prime\prime}}_{{\bf
k}}-\rho^{-1,m^{\prime}m^{\prime}}_{{\bf
k}^{\prime}})\delta(\omega^{m^{\prime\prime}}_{{\bf
k}}-\omega^{m^{\prime}}_{{\bf k}^{\prime}})\right].$

By solving Eq. 53, we have

$\displaystyle S^{1,mm^{\prime}}_{{\bf
k}}=\int_{0}^{\infty}dt^{\prime}e^{-iH_{0}({\bf
k})t^{\prime}}(D^{mm^{\prime}}_{{\bf k}{\bf
k}^{\prime}}+D^{{}^{\prime}mm^{\prime}}_{{\bf k}{\bf
k}^{\prime}})e^{iH_{0}({\bf k})t^{\prime}}$

After integrating over time, the off-diagonal density matrix in linear order
of the electric field yields

$\displaystyle S^{1,mm^{\prime}}_{{\bf k}}$ $\displaystyle=$
$\displaystyle-i\mathcal{P}\left(\frac{D^{mm^{\prime}}_{{\bf k}{\bf
k}^{\prime}}+D^{{}^{\prime}mm^{\prime}}_{{\bf k}{\bf
k}^{\prime}}}{\omega_{{\bf k}}^{m}-\omega_{{\bf
k}^{\prime}}^{m^{\prime}}}\right)$ $\displaystyle+$
$\displaystyle\pi\left(D^{mm^{\prime}}_{{\bf k}{\bf
k}^{\prime}}+D^{{}^{\prime}mm^{\prime}}_{{\bf k}{\bf
k}^{\prime}}\right)\delta(\omega_{{\bf k}}^{m}-\omega_{{\bf
k}^{\prime}}^{m^{\prime}})$

As a result, we can write the linear response of the density matrix in the
present of disorder potential accordingly

$\rho^{1,mm^{\prime}}_{{\bf k}{\bf k}^{\prime}}=\rho^{1,mm^{\prime}}_{{\bf
k}}\delta_{mm^{\prime}}+\rho^{1,mm^{\prime}}_{{\bf k}{\bf k}^{\prime}},$ (58)

where the first term describes the diagonal elements of the density matrix to
the first order in the electric field and finite frequency in long term limit
derived by Eq. 46, and the second term is off-diagonal density matrix to
finite frequency when $m\neq m^{\prime}$ is

$\displaystyle\rho^{1,mm^{\prime}}_{{\bf k}{\bf
k}}=-i\sum_{b\beta}\frac{\widetilde{D}^{b\beta,mm^{\prime}}_{{\bf k}{\bf
k}}+\widetilde{I}^{b\beta,mm^{\prime}}_{{\bf k}{\bf k}}}{(\omega^{m}_{{\bf
k}}-\omega^{m^{\prime}}_{{\bf
k}})-\omega_{\beta}}E^{b}_{\beta}e^{-i\omega_{\beta}t}$ (59)

with

$\displaystyle\widetilde{D}^{b\beta,mm^{\prime}}_{{\bf k}{\bf k}}$
$\displaystyle=$ $\displaystyle{ie}r^{b}_{mm^{\prime}}({\bf
k})\big{(}f(\omega^{m}_{{\bf k}})-f(\omega^{m^{\prime}}_{{\bf k}})\big{)}$
(60) $\displaystyle\widetilde{I}^{b\beta,mm^{\prime}}_{{\bf k}{\bf
k}^{\prime}}$ $\displaystyle=$ $\displaystyle{\pi
n_{i}}\sum_{m^{\prime\prime},{\bf k}^{\prime\prime}}u^{mm^{\prime}}_{{\bf
k}{\bf k}^{\prime}}u^{m^{\prime}m^{\prime\prime}}_{{\bf k}^{\prime}{\bf k}}$
$\displaystyle\big{[}(\widetilde{\rho}_{{\bf
k}}^{1b\beta,mm}-\widetilde{\rho}_{{\bf
k}^{\prime}}^{1b\beta,m^{\prime}m^{\prime}})\delta(\omega^{m}_{{\bf
k}}-\omega^{m^{\prime}}_{{\bf k}^{\prime}})$ $\displaystyle+$
$\displaystyle(\widetilde{\rho}_{{\bf
k}}^{1b\beta,m^{\prime\prime}m^{\prime\prime}}-\widetilde{\rho}_{{\bf
k}^{\prime}}^{1b\beta,m^{\prime}m^{\prime}})\delta(\omega^{m^{\prime\prime}}_{{\bf
k}}-\omega^{m^{\prime}}_{{\bf k}^{\prime}})\big{]}$

The off-diagonal matrix presenting in Eq. 59 consists of two main parts: First
part shown in Eq. 60 is independent of disorder character and is an intrinsic
band-structure property involving the entire Fermi sea [24, 23]. Second one
shown in Eq. III.1 is an extrinsic contribution containing both the band off-
diagonal and the band-diagonal parts of the density matrix. Obviously it is
finite in the weak disorder limit according to Eq. 50 and originates from the
collision kernel acting on the leading Fermi surface response.

###  III.2 Second-order Electric Response

In this subsection, we extend our results to the second order electric
response for the density matrix. To do so, we start from the equation of the
motion for ${\rho}^{2}$ accordingly

$\displaystyle\frac{d\langle\rho^{2}\rangle}{dt}+{i}[H_{0},\langle\rho^{2}\rangle]$
$\displaystyle+$ $\displaystyle J(\langle\rho^{2}\rangle)=$
$\displaystyle-{i}[H_{E},\langle\rho^{1}\rangle]-{i}[H_{E},\langle\rho^{2}\rangle],$

where the last term on the right hand side is beyond second order and is
negligible. To derive the diagonal term, also we can neglect the second term
on Eq. III.2 since it is purely off-diagonal in matrix space. As we derived in
the linear response, we must derive all the terms presented in Eq. III.2.
Hence the electric dependency is derived by

$\displaystyle{-i}\langle mk|[H_{E},\rho^{1}]|m^{\prime}k^{\prime}\rangle$
$\displaystyle=$ $\displaystyle{-e{\bf
E}}\cdot\left[\delta_{mm^{\prime}}\frac{\partial\langle\rho^{1}\rangle^{mm^{\prime}}_{{\bf
k}}}{\partial{{\bf k}}}\right.$ (63)
$\displaystyle\quad-\left.i[\mathcal{R}_{{\bf k}},\langle\rho^{1}\rangle_{{\bf
k}}]^{mm^{\prime}}\right],$

where $\mathcal{R}^{mm^{\prime}}_{\bf k}=\langle u^{m}_{{\bf
k}}|i\frac{\partial u^{m^{\prime}}_{{\bf k}}}{\partial{{\bf k}}}\rangle$ and
the first term is derived by inserting the diagonal term of density matrix
shown in Eq. 58 and take the derivation of the momentum to the second order.
Hence with doing some calculations the diagonal term to second order response
of the density matrix yields

$\displaystyle\rho_{m{\bf k}m{\bf k}}^{(2),b\beta c\gamma}$ $\displaystyle=$
$\displaystyle{e}\sum_{m^{\prime}{\bf k}^{\prime}}({\Gamma^{mm^{\prime}}_{{\bf
k}{\bf k}^{\prime}}})^{-1}\big{[}\partial_{k_{c}}\sum_{m^{\prime\prime}{\bf
k}^{\prime\prime}}\big{(}({\Gamma^{m^{\prime}m^{\prime\prime}}_{{\bf
k}^{\prime}{\bf k}^{\prime\prime}}})^{-1}\frac{\partial
f(\omega^{m^{\prime\prime}}_{{\bf k}^{\prime\prime}})}{\partial_{{\bf
k}^{\prime\prime}}}\big{)}$ (64) $\displaystyle-$
$\displaystyle{i}{e}\sum_{m^{\prime}{\bf
k}^{\prime}}{(\Gamma^{mm^{\prime}}_{{\bf k}{\bf
k}^{\prime}})}^{-1}\sum_{n}\big{(}r^{c}_{mn}B^{b}_{nm}-B^{b}_{mn}r^{c}_{nm}$
$\displaystyle\quad\qquad\qquad\qquad\qquad+r^{c}_{mn}\widetilde{I}^{{}^{\prime}b\beta}_{nm}-\widetilde{I}^{{}^{\prime}b\beta}_{mn}r^{c}_{nm}\big{)}\big{]}$

where

$\displaystyle
B^{b}_{nm}=\frac{ef_{mn}r^{b}_{mn}}{\omega^{mn}-\omega_{\beta}},\quad\widetilde{I}^{{}^{\prime}b\beta}_{nm}=\frac{\widetilde{I}^{b\beta}_{nm}}{\omega^{mn}-\omega_{\beta}},$
(65)

including $f_{nm}=f(\omega^{n}_{{\bf k}})-f(\omega^{m}_{{\bf k}})$ and
$\omega^{mn}=\omega^{m}_{{\bf k}}-\omega^{n}_{{\bf k}}$. For the off-diagonal
terms, we can start from Eq. III.2 with a point that in this part we have to
consider the second term on the left hand side since this term is purely off-
diagonal and need to be considered. Hence we have

$\displaystyle{i}[H_{0},\langle\rho^{2}\rangle]^{mm^{\prime}}$
$\displaystyle=$ $\displaystyle i\omega^{mm^{\prime}}_{{\bf k}{\bf
k}^{\prime}}\rho^{2,mm^{\prime}}_{{\bf k}{\bf k}^{\prime}}$ (66)

It is obviously that the above equation is purely of diagonal. The next term
on Eq. III.2 is written as

$\displaystyle-{i}[H_{E},\langle\rho^{(1)}\rangle]^{mm^{\prime\prime}}={e}\sum_{b\beta}\sum_{c\gamma}\big{[}\partial_{k_{c}}\big{(}B^{b}_{mm^{\prime\prime}}+\widetilde{I}^{{}^{\prime}b\beta}_{mm^{\prime\prime}}\big{)}$
$\displaystyle\qquad\qquad\quad-i\big{(}r^{c}_{mm^{\prime\prime}}-r^{c}_{m^{\prime\prime}m}\big{)}\sum_{m^{\prime}{\bf
k}^{\prime}}(\Gamma^{mm^{\prime}}_{{\bf k}{\bf
k}^{\prime}})^{-1}\frac{\partial f(\omega^{m^{\prime}}_{{\bf
k}^{\prime}})}{\partial{k^{{}^{\prime}b}}}$
$\displaystyle\qquad\quad\qquad+\sum_{m^{\prime}}\bigg{(}r^{c}_{mm^{\prime}}B^{b}_{m^{\prime}m^{\prime\prime}}-B^{b}_{mm^{\prime}}r^{c}_{m^{\prime}m^{\prime\prime}}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad+i(r^{c}_{mm^{\prime}}\widetilde{I}^{{}^{\prime}b\beta}_{m^{\prime}m^{\prime\prime}}-\widetilde{I}^{{}^{\prime}b\beta}_{mm^{\prime}}r^{c}_{m^{\prime}m^{\prime\prime}})\bigg{)}$
$\displaystyle\quad\qquad\qquad-\big{(}\omega^{c}_{mm}-\omega^{c}_{m^{\prime\prime}m^{\prime\prime}}\big{)}\big{(}iB^{b}_{mm^{\prime\prime}}+\widetilde{I}^{{}^{\prime}b\beta}_{mm^{\prime\prime}}\big{)}\big{]}$
$\displaystyle\qquad\qquad\qquad\qquad\times
E^{b}_{\beta}E^{c}_{\gamma}e^{-i(\omega_{\beta}+\omega_{\gamma})t}$
$\displaystyle\qquad\qquad\quad\quad\qquad\qquad\equiv\widetilde{\eta}^{b\beta
c\gamma}_{mm^{\prime\prime}}E^{b}_{\beta}E^{c}_{\gamma}e^{-i(\omega_{\beta}+\omega_{\gamma})t}$
(67)

In above equation, when $v_{imp}=0$, the latter term vanishes and we have

$\displaystyle-{i}[H_{E},\langle\rho^{(1)}\rangle]^{mm^{\prime\prime}}$
$\displaystyle=$
$\displaystyle{e}\sum\big{[}r^{c}_{mm^{\prime}}B^{b}_{m^{\prime}m^{\prime\prime}}-B^{b}_{mm^{\prime}}r^{c}_{m^{\prime}m^{\prime\prime}}\big{]}$
(68) $\displaystyle\qquad\qquad\times
E^{b}_{\beta}E^{c}_{\gamma}e^{-i(\omega_{\beta}+\omega_{\gamma})t}$

This is in agreement with Refs. [24, 23] when the second order of the electric
response is recorded in clean limit. Now the collision integral is founded as

$\displaystyle J_{imp}(\rho^{2})$ $\displaystyle=$ $\displaystyle{\pi
n_{i}}\sum_{m^{\prime}{\bf k}^{\prime}}u_{{\bf k}{\bf
k}^{\prime}}^{mm^{\prime}}u_{{\bf k}^{\prime}{\bf k}}^{m^{\prime}m}$ (69)
$\displaystyle\times$ $\displaystyle\big{[}(\rho^{2,b\beta c\gamma}_{m{\bf
k}m{\bf k}}-\rho^{2,b\beta c\gamma}_{m^{\prime}{\bf k}^{\prime}m^{\prime}{\bf
k}^{\prime}})\delta(\omega^{m}_{{\bf k}}-\omega^{m^{\prime}}_{{\bf
k}^{\prime}})$ $\displaystyle+$ $\displaystyle(\rho^{2,b\beta
c\gamma}_{m^{\prime\prime}{\bf k}m^{\prime\prime}{\bf k}}-\rho^{2,b\beta
c\gamma}_{m^{\prime}{\bf k}^{\prime}m^{\prime}{\bf
k}^{\prime}})\delta(\omega^{m^{\prime\prime}}_{{\bf
k}}-\omega^{m^{\prime}}_{{\bf k}^{\prime}})\big{]}$ $\displaystyle\times$
$\displaystyle
E^{b}_{\beta}E^{c}_{\gamma}e^{-i(\omega_{\beta}+\omega_{\gamma})t}$
$\displaystyle\equiv$ $\displaystyle\mathcal{I}^{b\beta c\gamma}_{m{\bf
k}m^{\prime}{\bf
k}^{\prime}}E^{b}_{\beta}E^{c}_{\gamma}e^{-i(\omega_{\beta}+\omega_{\gamma})t}$

Now we are able to back to Eq. III.2. By inserting Eqs. 69 and 67 to Eq.
III.2, we have

$\displaystyle\rho^{2mm^{\prime\prime}}_{{\bf
k}}=-\sum_{b\beta}\sum_{c\gamma}\frac{\widetilde{\eta}^{b\beta
c\gamma}_{mm^{\prime\prime}}+\widetilde{I}^{b\beta
c\gamma}_{mm^{\prime\prime}}}{(\omega_{{\bf k}}^{m}-\omega_{{\bf
k}}^{m^{\prime\prime}})-\omega_{\Sigma}}E^{b}_{\beta}E^{c}_{\gamma}e^{-i\omega_{\Sigma}t}$
(70)

with $\omega_{\Sigma}=\omega_{\beta}+\omega_{\gamma}$. In order to compare our
results with previous results derived by Refs. [24, 23] when $V_{imp}=0$, we
have

$\displaystyle\widetilde{\eta}^{b\beta
c\gamma}_{mm^{\prime\prime}}={e}\big{(}B^{b}_{mm^{\prime\prime},c}+\sum_{m^{\prime}}(r^{c}_{mm^{\prime}}B^{b}_{m^{\prime}m^{\prime\prime}}-B^{b}_{mm^{\prime}}r^{c}_{m^{\prime}m^{\prime\prime}})\big{)}$
(71)

where

$\displaystyle B^{b}_{mm^{\prime\prime},c}=\frac{\partial
B^{b}_{mm^{\prime\prime}}}{\partial{k^{c}}}-i(\omega^{mm}_{c}({\bf
k})-\omega^{m^{\prime\prime}m^{\prime\prime}}_{c}({\bf k})).$ (72)

By inserting Eq. 71 in Eq. 70, when $V_{imp}=0$ we integrate to find

$\displaystyle\rho^{2mm^{\prime\prime}}_{{\bf k}}$ $\displaystyle=$
$\displaystyle\frac{-ie}{(\omega^{m}_{{\bf k}}-\omega^{m^{\prime\prime}}_{{\bf
k}})-\omega_{\Sigma}}\sum_{b\beta}\sum_{c\gamma}\big{(}B^{b}_{mm^{\prime\prime},c}$
$\displaystyle+\sum_{m^{\prime}}(r^{c}_{mm^{\prime}}B^{b}_{m^{\prime}m^{\prime\prime}}-B_{mm^{\prime}}^{b}r^{c}_{m^{\prime}m^{\prime\prime}})\big{)}E^{b}_{\beta}E^{c}_{\gamma}e^{-i\omega_{\Sigma}t}$

This is exactly in agreement with the previous results reported in [24, 23].
Hence, density matrix to the second order response reads as

$\rho^{2mm^{\prime}}_{{\bf k}{\bf k}^{\prime}}=\rho^{2mm^{\prime}}_{{\bf
k}}\delta_{mm^{\prime}}+\rho^{2mm^{\prime}}_{{\bf k}{\bf k}^{\prime}},$ (74)

where the first part is the diagonal part derived in Eq. 64 and the second one
is off-diagonal reported explicitly on Eq. 70.

##  IV Disordered Rashba model for a nonmagnetic semiconductor

The Hamiltonian model for a 2DEG in the presence of Rashba SOC (RSOC) reads

$\displaystyle
H=\frac{k^{2}}{2m}+\alpha(k_{y}\sigma^{x}-k_{x}\sigma^{y})+V_{imp}(\bf r)$
(75)

with $m^{\star}$ the electron effective mass, $\alpha$ the spin-orbit strength
and $\sigma^{i}$ the Pauli matrices $i=x,y,z$. For Hamiltonian the
eigenvectors are

$\displaystyle|mk\rangle\equiv|\pm{\bf k}\rangle=|u_{{\bf
k}}^{\pm}\rangle=\frac{1}{\sqrt{2}}\begin{pmatrix}e^{-i\theta}\\\ \mp
i\end{pmatrix}$ (76)

where $\theta$ is the polar angle of the wave vector. The energy splitting for
each band read as

$\omega^{\pm}_{{\bf k}}=\frac{k^{2}}{2m^{\star}}\pm\alpha k.$ (77)

where $m^{\star}$ is the effective mass of electron. In order to derive the
relaxation time for interband and intraband, we could consider two different
limits, short range impurity potential and long range impurity potential. In
the short range of impurity we assume that the $V^{mm^{\prime}}_{{\bf k}{\bf
k}^{\prime}}$ does not depend on ${\bf q}={\bf k}-{\bf k}^{\prime}$ and one
can replace $V^{mm^{\prime}}_{{\bf k}{\bf k}^{\prime}}$ to $v_{0}$. On other
side, when we assumed the impurity is large in its potential, according to Eq.
76 we can derive

$\displaystyle\langle V^{\pm\mp}_{{\bf k}{\bf k}^{\prime}}V^{\mp\pm}_{{\bf
k}^{\prime}{\bf k}}\rangle$ $\displaystyle=$
$\displaystyle\frac{{v}^{2}_{0}}{2}(1-cos\gamma)$ (78) $\displaystyle\langle
V^{\pm\pm}_{{\bf k}{\bf k}^{\prime}}V^{\mp\mp}_{{\bf k}^{\prime}{\bf
k}}\rangle$ $\displaystyle=$ $\displaystyle\frac{{v}^{2}_{0}}{2}(1+cos\gamma)$
(79) $\displaystyle\langle V^{++}_{{\bf k}{\bf k}^{\prime}}V^{+-}_{{\bf
k}^{\prime}{\bf k}}\rangle$ $\displaystyle=$
$\displaystyle\frac{{v}^{2}_{0}}{2}sin\gamma$ (80) $\displaystyle\langle
V^{+-}_{{\bf k}{\bf k}^{\prime}}V^{--}_{{\bf k}^{\prime}{\bf k}}\rangle$
$\displaystyle=$ $\displaystyle\frac{-{V}^{2}_{0}}{2}sin\gamma$ (81)

In order to derive the matrix of relaxation time, the starting point is
getting back to Eq. 42 for $\mathcal{K}^{mm^{\prime}}_{{\bf k}{\bf
k}^{\prime}}$. In short range impurity potential, as the scattering amplitude
is constant and we can replace $V_{{\bf k}{\bf k}^{\prime}}^{mm^{\prime}}$ to
${v}_{0}$ as

$\displaystyle\frac{1}{\tau_{mm^{\prime}}}=\mathcal{K}^{mm^{\prime}}_{{\bf
k}{\bf k}^{\prime}}={2\pi n_{i}{v}^{2}_{0}}\sum_{{\bf
k}^{\prime}}(1-cos\gamma)\delta(\omega^{m}_{{\bf k}}-\omega^{m^{\prime}}_{{\bf
k}^{\prime}})$ (82)

that with integrating over momentum we read

$\displaystyle\frac{1}{\tau_{mm^{\prime}}}={n_{i}{v}^{2}_{0}}\int
k^{\prime}dk^{\prime}\delta(\omega^{m}_{{\bf k}}-\omega^{m^{\prime}}_{{\bf
k}^{\prime}}).$ (83)

In this case for the relaxation matrix, only term survived when $m=m^{\prime}$
in diagonal terms and for off-diagonal terms $m\neq m^{\prime}$ since $\gamma$
goes to zero, these terms vanishes. By considering the long-range impurity
potential, by solving Eq. 42 to the case when the scattering amplitude
dependents on $\gamma$ we can show that

$\displaystyle\frac{1}{\tau_{mm^{\prime}}}$ $\displaystyle=$
$\displaystyle{2\pi n_{i}{v}^{2}_{0}}\big{(}\frac{1}{4}\sum_{{\bf
k}^{\prime}}\delta(\omega^{m}_{{\bf k}}-\omega^{m^{\prime}}_{{\bf
k}^{\prime}})$ $\displaystyle\qquad\qquad\qquad-\frac{1}{4}\sum_{{\bf
k}^{\prime}}cos(2\gamma)\delta(\omega^{m}_{{\bf k}}-\omega^{m^{\prime}}_{{\bf
k}^{\prime}})\big{)}$

By taking the integration over the momentum, Eq. IV reads as

$\displaystyle\frac{1}{\tau_{\pm\pm}}\equiv\frac{1}{\tau_{\pm}}$
$\displaystyle=$ $\displaystyle{n_{i}{v}^{2}_{0}}\big{(}\frac{1}{4}\int
k^{\prime}dk^{\prime}\delta(\omega^{\pm}_{{\bf k}}-\omega^{\pm}_{{\bf
k}^{\prime}})$ (85) $\displaystyle\qquad\qquad+\frac{3}{4}\int
k^{\prime}dk^{\prime}\delta(\omega^{\pm}_{{\bf k}}-\omega^{\mp}_{{\bf
k}^{\prime}})\big{)}$

In the long impurity regime, our equation shows that when $\gamma$ is taken
into account, the interband relaxation time $m=m^{\prime}$ is three times
smaller than the intraband relaxation time $m\neq m^{\prime}$. In order make
our method as simple as possible, we confine the system to the short-range
impurity potential. From now on, we limit the system to the short-range
impurities and as considered before in this regime, only one relaxation time
for each band is considered. More precisely, in the presence of RSOC, the SO
splitting and the disorder broadening are much smaller than the Fermi energy
$\omega_{F}$

$\displaystyle\omega_{F}\gg\frac{1}{\tau_{0}},\quad\epsilon_{F}\ll 2\alpha
p_{F}$ (86)

where $p_{F}$ is the momentum at the Fermi surface. In the combination of the
disorder broadening and spin-orbit splitting also we can define two different
regime depending on which one dominates. The first one is the diffusive
regime, corresponding to a high impurity concentration, i.e.,
$(\alpha/v_{f})\epsilon_{F}\tau\ll 1$. The second regime, which occurs at
$(\alpha/v_{f})\epsilon_{F}\tau\gg 1$, goes beyond the diffusive limit and
describes the opposite situation of a relatively low concentration of
impurities, where the spin-relaxation time is close to $\tau_{0}$. For short-
range impurities, when the intraband of the relaxation matrix is neglected,
and we can use Eq. 83 to derive the relaxation time in the presence of RSOC.
By inserting Eq. 83 in Eq. 83, and taking the integration over momentum, we
receive [19, 21]

$\displaystyle\frac{1}{\tau_{\pm}}\equiv\mathcal{K}_{\pm{\bf k},\pm{\bf
k}^{\prime}}\approx\frac{1}{\tau_{0}}\left(\frac{m^{\star}}{1\pm\frac{m^{\star}}{k}}\right)$
(87)

In a single matrix in a 2D we can write the relaxation time as

$\displaystyle\hat{\tau}=\tau_{0}\left(\sigma^{0}-\frac{m^{\star}\alpha}{k}\sigma^{z}\right)$
(88)

By adding above equation inside Eq. (50), we have

$\displaystyle\rho_{{\bf k}}^{(-1),\pm\pm}=-\frac{e{\bf E}\cdot\hat{{\bf
k}}}{m^{\star}}\tau_{\pm}k_{\pm}\delta(\omega^{\pm}_{k}-\omega_{F}),$ (89)

where we have used

$\displaystyle\frac{\partial f_{0}(\omega^{m^{\prime}}_{{\bf
k}^{\prime}})}{\partial\omega^{m^{\prime}}_{{\bf
k}^{\prime}}}=-\delta(\omega^{\pm}_{k}-\omega_{F})$ (90)

and to the lowest order in the spin-orbit splitting, for the momentum we have

$\displaystyle k_{\pm}=k(1\pm\frac{\alpha m^{\star}}{k}).$ (91)

Hence we have

$\displaystyle\rho_{{\bf k}}^{(-1),mm}=\frac{em^{\star}\alpha{\bf
E}\cdot\hat{{\bf
k}}\tau_{0}}{k}\sigma^{z}\left(\delta(k-k_{F})-k\frac{\partial}{\partial{k}}\delta(k-k_{F})\right)$

In the eigenstate basic, the spin operator along $y$-direction is given by
$s_{y}=-1/2(cos\theta\sigma^{z}+sin\theta\sigma^{y})$. When the electric field
applied along ${\bf E}=E\hat{x}$, the expectation value of the spin
polarization $s^{y}$ is obtained as

$\displaystyle\langle s^{y}\rangle=\textnormal{Tr}[s^{y}\rho_{{\bf
k}}^{(-1),mm}]=-N_{0}e\alpha\tau_{0}E_{x}$ (93)

where the trace symbol includes the summation over the spin indices and
momentum. This is the standard term for the spin polarization reported by Ref.
[25].

##  V CONCLUSIONS

We have derived a general theory of evaluating the second order response of a
crystal to an electric field accounting for both interband and intraband
contributions. Our theory automatically includes momentum-space Berry phase
effect which is often discussed in the context of transport theory. We have
derived a quantum kinetic equation by applying the Winger transformation to
the quantum Liouvile equation. From this equation, we have obtained that the
intraband response is derived by the band diagonal part of the density matrix
proportional to the strengthen of disorder potential responsible for
properties like longitudinal conductivity. As the system confinement to the
Born approximation, the relaxation time matrix is proportional to the disorder
potential to leading order. Interband response is captured by off-diagonal
part of the density matrix and has sizable contributions to the repose of a
crystals. The interplay between diagonal and off-diagonal parts of the density
matrix is important in determining the “guage invariant” operator Eq. 20 which
is a macroscopic polarization potential to calculate the macroscopic current-
density operator $\textbf{J}(t)=d{\textbf{P}}/dt$. This analysis we presented
here can be played essential role in future, more accurate calculations of
nonlinear optical response.

## References

  * Boyd [1992] R. W. Boyd, Nonlinear optics, academic press, San Diego 155 (1992).
  * Chang  _et al._ [2016] K. Chang, J. Liu, H. Lin, N. Wang, K. Zhao, A. Zhang, F. Jin, Y. Zhong, X. Hu, W. Duan, Q. Zhang, L. Fu, Q.-K. Xue, X. Chen, and S.-H. Ji, Discovery of robust in-plane ferroelectricity in atomic-thick snte, Science 353, 274 (2016), https://www.science.org/doi/pdf/10.1126/science.aad8609 .
  * Fridkin [2001] V. Fridkin, Crystallogr. rep (2001).
  * Osterhoudt  _et al._ [2019] G. B. Osterhoudt, L. K. Diebel, M. J. Gray, X. Yang, J. Stanco, X. Huang, B. Shen, N. Ni, P. J. W. Moll, Y. Ran, and K. S. Burch, Colossal mid-infrared bulk photovoltaic effect in a type-i weyl semimetal, Nature Materials 18, 471 (2019).
  * Neupane  _et al._ [2016] M. Neupane, I. Belopolski, M. M. Hosen, D. S. Sanchez, R. Sankar, M. Szlawska, S.-Y. Xu, K. Dimitri, N. Dhakal, P. Maldonado, P. M. Oppeneer, D. Kaczorowski, F. Chou, M. Z. Hasan, and T. Durakiewicz, Observation of topological nodal fermion semimetal phase in zrsis, Phys. Rev. B 93, 201104 (2016).
  * McIver  _et al._ [2012] J. W. McIver, D. Hsieh, H. Steinberg, P. Jarillo-Herrero, and N. Gedik, Control over topological insulator photocurrents with light polarization, Nature Nanotechnology 7, 96 (2012).
  * Rappe  _et al._ [2017] A. M. Rappe, I. Grinberg, and J. E. Spanier, Getting a charge out of hybrid perovskites, Proceedings of the National Academy of Sciences 114, 7191 (2017), https://www.pnas.org/doi/pdf/10.1073/pnas.1708154114 .
  * Boyd [2003] R. Boyd, Non linear optics, academic press, san diego, Cal., USA , 528 (2003).
  * Rabe  _et al._ [2007a] K. M. Rabe, M. Dawber, C. Lichtensteiger, C. H. Ahn, and J.-M. Triscone, Modern physics of ferroelectrics:essential background, in  _Physics of Ferroelectrics: A Modern Perspective_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 2007) pp. 1–30.
  * Spanier  _et al._ [2016] J. E. Spanier, V. M. Fridkin, A. M. Rappe, A. R. Akbashev, A. Polemi, Y. Qi, Z. Gu, S. M. Young, C. J. Hawley, D. Imbrenda, G. Xiao, A. L. Bennett-Jackson, and C. L. Johnson, Power conversion efficiency exceeding the shockley-queisser limit in a ferroelectric insulator, Nature Photonics 10, 611 (2016).
  * Rakita  _et al._ [2017] Y. Rakita, O. Bar-Elli, E. Meirzadeh, H. Kaslasi, Y. Peleg, G. Hodes, I. Lubomirsky, D. Oron, D. Ehre, and D. Cahen, Tetragonal ch¡sub¿3¡/sub¿nh¡sub¿3¡/sub¿pbi¡sub¿3¡/sub¿ is ferroelectric, Proceedings of the National Academy of Sciences 114, E5504 (2017), https://www.pnas.org/doi/pdf/10.1073/pnas.1702429114 .
  * Culcer  _et al._ [2017] D. Culcer, A. Sekine, and A. H. MacDonald, Interband coherence response to electric fields in crystals: Berry-phase contributions and disorder effects, Physical Review B 96, 035106 (2017).
  * Sekine  _et al._ [2017] A. Sekine, D. Culcer, and A. H. MacDonald, Quantum kinetic theory of the chiral anomaly, Phys. Rev. B 96, 235134 (2017).
  * Vasko and Raichev [2006a] F. T. Vasko and O. E. Raichev,  _Quantum Kinetic Theory and Applications: Electrons, Photons, Phonons_ (Springer Science & Business Media, 2006).
  * Yang  _et al._ [2017] X. Yang, K. Burch, and Y. Ran, Divergent bulk photovoltaic effect in weyl semimetals, arXiv preprint arXiv:1712.09363 (2017).
  * Allen [1978] P. B. Allen, New method for solving boltzmann’s equation for electrons in metals, Phys. Rev. B 17, 3725 (1978).
  * Karube  _et al._ [2016] S. Karube, K. Kondou, and Y. Otani, Experimental observation of spin-to-charge current conversion at non-magnetic metal interfaces, Applied Physics Express 9, 033001 (2016).
  * Chen  _et al._ [2016] L. Chen, M. Decker, M. Kronseder, R. Islinger, M. Gmitra, D. Schuh, D. Bougeard, J. Fabian, D. Weiss, and C. H. Back, Robust spin-orbit torque and spin-galvanic effect at the fe/gaas (001) interface at room temperature, Nature Communications 7, 13802 (2016).
  * Gorini  _et al._ [2017] C. Gorini, A. Maleki Sheikhabadi, K. Shen, I. V. Tokatly, G. Vignale, and R. Raimondi, Theory of current-induced spin polarization in an electron gas, Phys. Rev. B 95, 205424 (2017).
  * Shen  _et al._ [2014] K. Shen, G. Vignale, and R. Raimondi, Microscopic theory of the inverse edelstein effect, Phys. Rev. Lett. 112, 096601 (2014).
  * Maleki Sheikhabadi  _et al._ [2018] A. Maleki Sheikhabadi, I. Miatka, E. Y. Sherman, and R. Raimondi, Theory of the inverse spin galvanic effect in quantum wells, Phys. Rev. B 97, 235412 (2018).
  * Blount [1962] E. Blount, Formalisms of band theory (Academic Press, 1962) pp. 305–373.
  * Sipe and Shkrebtii [2000] J. E. Sipe and A. I. Shkrebtii, Second-order optical response in semiconductors, Phys. Rev. B 61, 5337 (2000).
  * Sipe and Ghahramani [1993] J. E. Sipe and E. Ghahramani, Nonlinear optical response of semiconductors in the independent-particle approximation, Phys. Rev. B 48, 11705 (1993).
  * Edelstein [1990] V. Edelstein, Spin polarization of conduction electrons induced by electric current in two-dimensional asymmetric electron systems, Solid State Communications 73, 233 (1990).
  * Rabe  _et al._ [2007b] K. M. Rabe, K. M. Rabe, C. H. Ahn, and J.-M. Triscone, Physics of ferroelectrics: a modern perspective (2007b).
  * Rabe  _et al._ [2007c] K. M. Rabe, K. M. Rabe, C. H. Ahn, and J.-M. Triscone, Physics of ferroelectrics: a modern perspective (2007c).
  * Vasko and Raichev [2006b] F. T. Vasko and O. E. Raichev,  _Quantum Kinetic Theory and Applications: Electrons, Photons, Phonons_ (Springer Science & Business Media, 2006).
  * Ma  _et al._ [2017] Q. Ma, S.-Y. Xu, C.-K. Chan, C.-L. Zhang, G. Chang, Y. Lin, W. Xie, T. Palacios, H. Lin, S. Jia, P. A. Lee, P. Jarillo-Herrero, and N. Gedik, Direct optical detection of weyl fermion chirality in a topological semimetal, Nature Physics 13, 842 (2017).