# Supplementary information: Nonlinear Hall effect from long-lived valley-
polarizing relaxons

Jae-Mo Lihm<EMAIL_ADDRESS>Cheol-Hwan Park<EMAIL_ADDRESS>Department of Physics and Astronomy, Seoul National University, Seoul 08826,
Korea,  Center for Correlated Electron Systems, Institute for Basic Science,
Seoul 08826, Korea,  Center for Theoretical Physics, Seoul National
University, Seoul 08826, Korea

Supplementary Information: Nonlinear Hall effect from long-lived, valley-
polarizing relaxons

##  Appendix A Boltzmann transport equation and the relaxon decomposition

In the presence of a uniform, time-dependent electric field
${\bm{\mathbf{E}}}(t)$, the BTE in the linear response regime reads [1]

$\frac{\partial\,\delta\\!f_{i}(t)}{\partial
t}=e{\bm{\mathbf{E}}}(t)\cdot{\bm{\mathbf{v}}}_{i}\,f^{(0)\prime}_{i}-\sum_{j}S_{ij}\,\delta\\!f_{j}(t).$
(1)

We use the shorthand notation for the electron eigenstate indices
$i=(n,{\bm{\mathbf{k}}})$ and $j=(n^{\prime},{\bm{\mathbf{k}}}^{\prime})$. The
scattering matrix $S$ is defined as

$S_{ij}=\tau^{-1}_{i}\delta_{ij}-\tau^{-1}_{i\leftarrow j}(1-\delta_{ij})$ (2)

where $\tau^{-1}_{j\leftarrow i}$ is the rate of the scattering from state $i$
to $j$, and

$\tau^{-1}_{i}=\sum_{j\neq i}\tau^{-1}_{j\leftarrow i}$ (3)

is the quasiparticle lifetime of state $i$. The first term on the right-hand
side of Eq. (2) describes the “scattering-out” process in which the electron
at state $i$ decays to another state, while the second term describes the
“scattering-in” process in which an electron from another state scatter into
state $i$.

The partial decay rate for electron-phonon scattering is [2, 3]

$\displaystyle\tau^{-1}_{i\leftarrow j}=\frac{2\pi}{\hbar N_{k}}$
$\displaystyle\sum_{\nu}\absolutevalue{g_{ji\nu}}^{2}\Big{[}\delta(\varepsilon_{i}-\varepsilon_{j}-\hbar\omega_{\mathbf{q}\nu})(n^{(0)}_{\mathbf{q}\nu}+f^{(0)}_{i})$
$\displaystyle+$
$\displaystyle\delta(\varepsilon_{i}-\varepsilon_{j}+\hbar\omega_{\mathbf{q}\nu})(n^{(0)}_{\mathbf{q}\nu}+1-f^{(0)}_{i})\Big{]},$
(4)

where $\varepsilon_{i}$ is the electron energy,
${\bm{\mathbf{q}}}={\bm{\mathbf{k}}}^{\prime}-{\bm{\mathbf{k}}}$ the momentum
transferred to the phonon, $\omega_{\mathbf{q}\nu}$ the phonon energy,
$n^{(0)}_{\mathbf{q}\nu}$ the phonon occupation number, and

$g_{ji\nu}\equiv
g_{n^{\prime}n\nu}({\bm{\mathbf{k}}},{\bm{\mathbf{q}}})=\matrixelement{n^{\prime}{\bm{\mathbf{k}}}+{\bm{\mathbf{q}}}}{\partial_{\mathbf{q}\nu}V}{n{\bm{\mathbf{k}}}}$
(5)

the electron-phonon matrix element. Here, we have defined
$\ket{n{\bm{\mathbf{k}}}}$ the electron eigenstate with momentum
${\bm{\mathbf{k}}}$ and band $n$, and

$\partial_{\mathbf{q}\nu}V=\sum_{\ell,\kappa,\alpha}\sqrt{\frac{\hbar}{2M_{\kappa}\omega_{\mathbf{q}\nu}}}e^{i{\bm{\mathbf{q}}}\cdot{\bm{\mathbf{R}}}_{\ell}}u_{\kappa\alpha,\nu}({\bm{\mathbf{q}}})\frac{\partial
V}{\partial R_{\ell\kappa\alpha}}$ (6)

the phonon perturbation potential where $R_{\ell\kappa\alpha}$ is the
$\alpha$-th Cartesian component of the position of the atom $\kappa$ in the
$\ell$-th unit cell, $M_{\kappa}$ the mass of the atom, and
$u_{\kappa\alpha,\nu}({\bm{\mathbf{q}}})$ the phonon eigendisplacement. We
assume that the phonon occupation is fixed to its value at thermal
equilibrium, i. e. , we ignore the phonon drag effect [4, 5, 6].

The ordinary linear current for a given carrier distribution $\delta\\!f_{i}$
is

${\bm{\mathbf{J}}}_{\rm L}(t)=-\frac{e}{V_{\rm
cell}\,N_{k}}\sum_{i}{\bm{\mathbf{v}}}_{i}\,\delta\\!f_{i}(t).$ (7)

The expression for the NLH current, which leads to the NLH effect, follows by
replacing ${\bm{\mathbf{v}}}_{i}$ with the anomalous velocity
${\bm{\mathbf{v}}}_{i,{\rm anom.}}$ [7]:

${\bm{\mathbf{J}}}_{\rm NLH}(t)=-\frac{e}{V_{\rm
cell}\,N_{k}}\sum_{i}{\bm{\mathbf{v}}}_{i,{\rm anom.}}\,\delta\\!f_{i}(t)\,,$
(8)

where

${\bm{\mathbf{v}}}_{i,{\rm
anom.}}=\frac{e}{\hbar}\,{\bm{\mathbf{E}}}\times{\bm{\mathbf{\Omega}}}_{i}\,.$
(9)

Here, the Berry curvature ${\bm{\mathbf{\Omega}}}_{i}$ is defined as

${\bm{\mathbf{\Omega}}}_{n\mathbf{k}}=-\imaginary\sum_{m}{\bm{\mathbf{A}}}_{nm{\bm{\mathbf{k}}}}\times{\bm{\mathbf{A}}}_{mn{\bm{\mathbf{k}}}}$
(10)

where

${\bm{\mathbf{A}}}_{mn{\bm{\mathbf{k}}}}=i\innerproduct{u_{m\mathbf{k}}}{\frac{\partial
u_{n\mathbf{k}}}{\partial{\bm{\mathbf{k}}}}}$ (11)

is the Berry connection with $\ket{u_{m\mathbf{k}}}$ the periodic part of the
Bloch eigenstate. Note that the $\mathcal{O}(E^{2})$ term of the occupation
function is even in ${\bm{\mathbf{k}}}$ and thus does not contribute to the
conductivity of a time-reversal symmetric material whose band velocity and
Berry curvature are odd in ${\bm{\mathbf{k}}}$ [8].

Let us consider a monochromatic electric field ($\omega=0$ yields the DC
limit)

${\bm{\mathbf{E}}}(t)=\real{\bm{\mathbf{E}}}_{\omega}\,e^{i\omega t},$ (12)

where the linear response occupation can be written as

$\delta\\!f_{i}(t)=\real\sum_{b}E_{\omega}^{b}\,\delta\\!f^{b}_{i}(\omega)\,e^{i\omega
t}.$ (13)

Then, the BTE [Eq. (1)] becomes

$i\omega\delta\\!f^{b}_{i}(\omega)=e\,v^{b}_{i}\,f^{(0)\prime}_{i}-\sum_{j}S_{ij}\,\delta\\!f^{b}_{j}(\omega).$
(14)

Given the occupation, the linear current [Eq. (7)] reads

$\displaystyle J^{a}_{\rm L}(t)$ $\displaystyle=-\frac{e}{V_{\rm
cell}\,N_{k}}\sum_{i}v^{a}_{i}\,\real\sum_{b}E_{\omega}^{b}\,\delta\\!f^{b}_{i}(\omega)\,e^{i\omega
t}$ $\displaystyle=\real\sum_{b}\sigma^{ab}_{\rm
L}\,E^{b}_{\omega}\,e^{i\omega t},$ (15)

with the linear conductivity

$\sigma^{ab}_{\rm L}(\omega)=-\frac{e}{V_{\rm
cell}\,N_{k}}\sum_{i}v^{a}_{i}\,\delta\\!f^{b}_{i}(\omega).$ (16)

For the NLH current, we modify Eqs. (8) and (10) to use the finite-frequency
generalization of the Berry curvature [9]:

${\bm{\mathbf{\Omega}}}_{n\mathbf{k}}(\omega)=-\imaginary\sum_{m}\frac{(\varepsilon_{m\mathbf{k}}-\varepsilon_{n\mathbf{k}})^{2}}{(\varepsilon_{m\mathbf{k}}-\varepsilon_{n\mathbf{k}})^{2}-\hbar^{2}\omega^{2}}\,{\bm{\mathbf{A}}}_{nm{\bm{\mathbf{k}}}}\times{\bm{\mathbf{A}}}_{mn{\bm{\mathbf{k}}}}.$
(17)

Then, the NLH current is

$\displaystyle\quad J^{a}_{\rm NLH}(t)$
$\displaystyle=-\frac{e^{2}}{\hbar\,V_{\rm
cell}\,N_{k}}\sum_{i}({\bm{\mathbf{E}}}(t)\times{\bm{\mathbf{\Omega}}}_{i}(\omega))^{a}\real\sum_{b}E_{\omega}^{b}\,\delta\\!f^{b}_{i}(\omega)\,e^{i\omega
t}$ $\displaystyle=-\frac{e^{2}}{2\,\hbar\,V_{\rm
cell}\,N_{k}}\sum_{i}\sum_{c,d}\Big{\\{}\epsilon_{acd}\,[E^{c}_{\omega}\,e^{i\omega
t}+(E^{c}_{\omega})^{*}e^{-i\omega t}]\,\Omega^{d}_{i}(\omega)$
$\displaystyle\quad\quad\times\real\sum_{b}E_{\omega}^{b}\,\delta\\!f^{b}_{i}(\omega)\,e^{i\omega
t}\Big{\\}}$ $\displaystyle=\real\sum_{b,c}\sigma^{a;bc}_{\rm
NLH}(\omega)\,\big{[}E_{\omega}^{b}\,E^{c}_{\omega}\,e^{2i\omega
t}+E_{\omega}^{b}\,(E^{c}_{\omega})^{*}\big{]}$ (18)

with the AC NLH conductivity

$\sigma^{a;bc}_{\rm NLH}(\omega)=\frac{e^{2}}{2\,\hbar\,V_{\rm
cell}\,N_{k}}\sum_{d}\epsilon_{adc}\sum_{i}\Omega^{d}_{i}(\omega)\,\delta\\!f^{b}_{i}(\omega).$
(19)

Note that the anomalous current has two terms, the second-harmonic term with
frequency $2\omega$ and the static rectification term. In the BTE, the two
terms share the same conductivity tensor.

Now, we study the solution of the BTE in the relaxon basis. By construction,
$x_{i}=1$ is a left eigenvector of the scattering matrix with eigenvalue zero:

$\sum_{i}S_{ij}=\tau_{j}^{-1}-\sum_{i\neq j}\tau^{-1}_{i\leftarrow j}=0.$ (20)

This property easily follows from Eqs. (2) and (3). For the symmetrized
scattering matrix $\widetilde{S}$ [Eq. (LABEL:eq:stilde_def)], the
corresponding null eigenvector is

$\Theta_{i0}=\frac{1}{\sqrt{k_{\rm
B}\,T\,\kappa}}\sqrt{f^{(0)}_{i}\left(1-f^{(0)}_{i}\right)}=\sqrt{-\frac{1}{\kappa}f^{(0)^{\prime}}_{i}}\,.$
(21)

where

$\displaystyle\kappa=\frac{dn}{d\mu}$ $\displaystyle=\frac{1}{V_{\rm
cell}\,N_{k}}\sum_{i}\frac{df^{(0)}_{i}}{d\mu}$ $\displaystyle=\frac{1}{k_{\rm
B}\,T}\frac{1}{V_{\rm
cell}\,N_{k}}\sum_{i}f^{(0)}_{i}\left(1-f^{(0)}_{i}\right)$ (22)

is the charge compressibility. This eigenvector, which is one of the relaxons,
satisfies the normalization condition

$\frac{1}{V_{\rm cell}\,N_{k}}\sum_{i}\Theta_{i0}^{2}=\frac{1}{k_{\rm
B}\,T\,\kappa}\frac{1}{V_{\rm
cell}\,N_{k}}\sum_{i}f^{(0)}_{i}\left(1-f^{(0)}_{i}\right)=1.$ (23)

This relaxon describes the change in the electron occupation due to the
increase in the chemical potential, $\delta\\!f_{i}\propto df^{(0)}_{i}/d\mu$.
Such a perturbation to the occupation cannot be dissipated because of charge
conservation. Hence, the relaxon has a zero scattering rate.

We now define the effective velocity [10, 11] and the Berry curvature of the
relaxons, which is the sum of the band velocity and the Berry curvature
weighted by the relaxon eigenvector and $\Theta_{i0}$:

$\displaystyle\widetilde{{\bm{\mathbf{v}}}}_{p}$
$\displaystyle\equiv\frac{1}{V_{\rm
cell}\,N_{k}}\sum_{i}{\bm{\mathbf{v}}}_{i}\,\Theta_{i0}\,\Theta_{ip}$
$\displaystyle=\frac{1}{V_{\rm cell}\,N_{k}}\frac{1}{\sqrt{k_{\rm
B}\,T\,\kappa}}\sum_{i}{\bm{\mathbf{v}}}_{i}\,\sqrt{f^{(0)}_{i}\left(1-f^{(0)}_{i}\right)}\,\Theta_{ip}$
$\displaystyle=\frac{\sum_{i}{\bm{\mathbf{v}}}_{i}\,\sqrt{f^{(0)}_{i}\left(1-f^{(0)}_{i}\right)}\,\left(\Theta_{ip}/\sqrt{V_{\rm
cell}\,N_{k}}\right)}{\sqrt{\sum_{i}f^{(0)}_{i}\left(1-f^{(0)}_{i}\right)}}\,,$
(24) $\displaystyle\widetilde{{\bm{\mathbf{\Omega}}}}_{p}$
$\displaystyle\equiv\frac{1}{V_{\rm
cell}\,N_{k}}\sum_{i}{\bm{\mathbf{\Omega}}}_{i}\,\Theta_{i0}\,\Theta_{ip}$
$\displaystyle=\frac{1}{V_{\rm cell}\,N_{k}}\frac{1}{\sqrt{k_{\rm
B}\,T\,\kappa}}\sum_{i}{\bm{\mathbf{\Omega}}}_{i}\sqrt{f^{(0)}_{i}\left(1-f^{(0)}_{i}\right)}\Theta_{ip}$
$\displaystyle=\frac{\sum_{i}{\bm{\mathbf{\Omega}}}_{i}\,\sqrt{f^{(0)}_{i}\left(1-f^{(0)}_{i}\right)}\,\left(\Theta_{ip}/\sqrt{V_{\rm
cell}\,N_{k}}\right)}{\sqrt{\sum_{i}f^{(0)}_{i}\left(1-f^{(0)}_{i}\right)}}\,.$
(25)

The last lines of Eq. (A) and Eq. (A), which we obtained using Eq. (A),
clearly show that the effective velocity and the Berry curvature of the
relaxons are the weighted average of those quantities over the thermally
activated single-particle Bloch states near the Fermi level. Also, from the
normalization convention of relaxons, these effective quantities have the same
units as those of the single-particle eigenstates. Our normalization
convention is consistent with Ref. [11] if the charge compressibility is
replaced by the heat capacity.

Now, we rewrite the conductivities in the relaxon basis. We begin by writing
the electron distribution in the relaxon basis:

$\delta\\!f^{b}_{i}(\omega)=\sum_{p}\sqrt{f^{(0)}_{i}\left(1-f^{(0)}_{i}\right)}\,\Theta_{ip}\,\delta\\!\widetilde{f}^{b}_{p}(\omega).$
(26)

Transforming the BTE [Eq. (14)] to the relaxon basis yields

$i\omega\,\delta\\!\widetilde{f}^{b}_{p}(\omega)=-e\,\sqrt{\frac{\kappa}{k_{\rm
B}\,T}}\widetilde{v}_{p}^{b}-\Gamma_{p}\,\delta\\!\widetilde{f}^{b}_{p}(\omega).$
(27)

By solving this equation, we find the stationary-state relaxon occupation

$\delta\\!\widetilde{f}^{b}_{p}(\omega)=-e\,\sqrt{\frac{\kappa}{k_{\rm
B}\,T}}\,\frac{\widetilde{v}_{p}^{b}}{i\omega+\Gamma_{p}}.$ (28)

Transforming back to the electron eigenstate basis, we find

$\displaystyle\delta\\!f^{b}_{i}(\omega)$
$\displaystyle=-e\,\sqrt{\frac{\kappa}{k_{\rm
B}\,T}}\,\sqrt{f^{(0)}_{i}\left(1-f^{(0)}_{i}\right)}\,\sum_{p}\Theta_{ip}\,\frac{\widetilde{v}_{p}^{b}}{(i\omega+\Gamma_{p})}$
$\displaystyle=-e\,\kappa\,\Theta_{i0}\,\sum_{p}\Theta_{ip}\frac{\widetilde{v}_{p}^{b}}{(i\omega+\Gamma_{p})}\,.$
(29)

By substituting this occupation into Eqs. (16, 19), we find the AC
conductivities in the relaxon basis:

$\begin{gathered}\sigma^{ab}_{\rm
L}(\omega)=e^{2}\,\kappa\,\sum_{p}\widetilde{v}^{a}_{p}\,\widetilde{v}^{b}_{p}\,\frac{\Gamma_{p}-i\omega}{\Gamma_{p}^{2}+\omega^{2}},\\\
\sigma^{a;bc}_{\rm
NLH}(\omega)=-\frac{e^{3}\,\kappa}{2\,\hbar}\sum_{d}\epsilon_{adc}\,\sum_{p}\widetilde{\Omega}^{d}_{p}(\omega)\,\widetilde{v}^{b}_{p}\,\frac{\Gamma_{p}-i\omega}{\Gamma_{p}^{2}+\omega^{2}}.\end{gathered}$
(30)

Here, $\widetilde{\Omega}^{d}(\omega)$ is defined analogously to Eq. (A),
using the frequency-dependent Berry curvature $\Omega^{d}_{i}(\omega)$ [Eq.
(17)] instead of the static one $\Omega^{d}_{i}$ [Eq. (10)].

The crystal structure and the field and current directions considered in this
work are shown in Fig. S1.

Figure S1:  Crystal structure of $\alpha$-GeTe and the field and current
directions considered in this work.

##  Appendix B Discrete and continuous eigenvalue spectra of relaxons

In this section, we discuss the eigenvalue spectra of the scattering matrix.
In the limit of an infinitely dense $k$-point grid, the eigenvalues of the
scattering matrix, i.e., the relaxon eigenvalues, are discrete below
$\tau^{-1}_{\rm max}=(\max_{i}\tau_{i})^{-1}$ and continuous above
$\tau^{-1}_{\rm max}$.

Figure S2 compares the relaxon velocity [Fig. S2(a)], relaxon Berry curvature
[Fig. S2(b)], and relaxon contribution to the conductivity and its cumulative
sum [Fig. S2(c,d)] obtained at two different $k$-point grids. Changing the
$k$-point grid from 80$\times$80$\times$80 to 100$\times$100$\times$100
roughly doubles the $k$-point density ($100^{3}/80^{3}=1.95$). Figures S2(a,b)
show that the relaxon velocity and Berry curvature are approximately
$N_{k}$-indepdent in the discrete regime, while decreasing with $N_{k}$ in the
continuum regime. Figures S2(c,d) show that the individual contributions to
the conductivity of the relaxons in the discrete regime
($\Gamma<\tau_{\mathrm{max}}^{-1}$) remain approximately unchanged, while
those in the continuum regime ($\Gamma>\tau_{\mathrm{max}}^{-1}$) become
approximately halved. The cumulative distribution remains almost invariant as
the number of relaxons in the continuum regime is doubled.

Figure S2:  (a) Relaxon effective velocity and (b) relaxon effective Berry
curvature obtained for two different $k$-point grids. (c, d) Relaxon
decomposition and its cumulative sum for the (a) linear and (b) NLH
conductivity of GeTe. We plot only the individual contributions whose absolute
value is greater than $10^{-3}\,$%.

##  Appendix C Temperature dependence of the NLH effect

Figure S3 shows the temperature dependence of the NLH effect (current
reponsivity is the ratio between the NLH conductivity and the linear
conductivity) as well as its relaxon decomposition. We find that the current
responsivitiy decreases with temperature [Fig. S3(a)]. The characteristics of
the relaxon decomposition, namely the sharp frequency dependence of the AC
current responsivity [Fig. S3(b)], and the large contribution of the discrete
relaxons to the NLH conductivity but not to the linear conductivity [Fig.
S3(c,d)] is present at all temperatures.

Figure S3:  (a) Temperature dependence of the current responsivity of GeTe.
(b) Frequency dependence of the current responsivity at various temperatures.
(c) Cumulative sum of the relaxon-decomposed linear conductivity at various
temperatures. (d) Same as (c), but for the NLH conductivity. All results are
for $n=10^{19}~{}\mathrm{cm}^{-3}$.

##  Appendix D Scattering diagnostics

In this section, we develop the “scattering diagnostics” method which allows
one to investigate how each phonon mode contributes to the conductivity. The
name “scattering diagnostics” is inspired by the “fluctuation diagnostics”
method [12] which decomposes self-energy into the contribution of individual
fluctuations. Here, we decompose the scattering matrix into the contributions
of individual phonon wavevectors:

$S=\sum_{{\bm{\mathbf{q}}}}S^{{\bm{\mathbf{q}}}}.$ (31)

To study the effect of a single phonon wavevector ${\bm{\mathbf{q}}}$ on
transport, we subtract its contribution from the scattering matrix and solve
the modified BTE. The resulting modified BTE is

$0=e\,v^{b}_{i}\,f^{(0)\prime}_{i}-\sum_{j}(S_{ij}-S_{ij}^{{\bm{\mathbf{q}}}})(\delta\\!f^{b}_{j}-\Delta^{\bm{\mathbf{q}}}\delta\\!f^{b}_{j}),$
(32)

where $\Delta^{\bm{\mathbf{q}}}\delta\\!f^{b}_{j}$ is the change in the
stationary-state electron distribution due to the scattering with wavevector
${\bm{\mathbf{q}}}$. Solving for $\Delta^{\bm{\mathbf{q}}}\delta\\!f^{b}_{j}$,
we find

$\sum_{j}(S_{ij}-S_{ij}^{{\bm{\mathbf{q}}}})\,\Delta^{\bm{\mathbf{q}}}\delta\\!f^{b}_{j}=-\sum_{j}S_{ij}^{{\bm{\mathbf{q}}}}\,\delta\\!f^{b}_{j}.$
(33)

Given the solution to this linear equation, we calculate the conductivity
using Eq. (LABEL:eq:sigma), but with
$\Delta^{\bm{\mathbf{q}}}\delta\\!f^{b}_{j}$ instead of $\delta\\!f^{b}_{j}$:

$\begin{gathered}\Delta^{\bm{\mathbf{q}}}\sigma^{ab}_{\rm L}=-\frac{e}{V_{\rm
cell}\,N_{k}}\sum_{i}v^{a}_{i}\,\Delta^{\bm{\mathbf{q}}}\delta\\!f^{b}_{i},\\\
\Delta^{\bm{\mathbf{q}}}\sigma^{a;bc}_{\rm NLH}=\frac{e^{2}}{2\hbar V_{\rm
cell}\,N_{k}}\sum_{d}\epsilon_{adc}\sum_{i}\Omega^{d}_{i}\,\Delta^{\bm{\mathbf{q}}}\delta\\!f^{b}_{i}.\end{gathered}$
(34)

The resulting momentum-resolved conductivity describes how sensitive the
conductivity is to the given scattering mechanism. A negative (positive) value
indicates that the scattering tends to suppress (enhance) the conductivity.

If $S^{{\bm{\mathbf{q}}}}$ is much smaller than the total scattering matrix
$S$, which is satisfied in the usual setup for first-principles calculations
because $N_{k}\gg 1$, one can approximate Eq. (33) by neglecting the
$\mathcal{O}((S^{{\bm{\mathbf{q}}}})^{2})$ term:

$\sum_{j}S_{ij}\,\Delta^{\bm{\mathbf{q}}}\delta\\!f^{b}_{j}\approx-\sum_{j}S_{ij}^{{\bm{\mathbf{q}}}}\,\delta\\!f^{b}_{j}.$
(35)

In this approximation, the sum of $\Delta^{\bm{\mathbf{q}}}\delta\\!f^{b}_{j}$
for all momenta ${\bm{\mathbf{q}}}$ satisfies

$\sum_{j}S_{ij}\sum_{\bm{\mathbf{q}}}\Delta^{\bm{\mathbf{q}}}\delta\\!f^{b}_{j}\approx-\sum_{j}S_{ij}\,\delta\\!f^{b}_{j}$
(36)

and so

$\sum_{\bm{\mathbf{q}}}\Delta^{\bm{\mathbf{q}}}\delta\\!f^{b}_{j}=-\delta\\!f^{b}_{j}\,.$
(37)

Thus, interestingly, we find

$\sum_{\bm{\mathbf{q}}}\Delta^{\bm{\mathbf{q}}}\sigma^{ab}_{\rm
L}=-\sigma^{ab}_{\rm
L},\quad\sum_{\bm{\mathbf{q}}}\Delta^{\bm{\mathbf{q}}}\sigma^{a;bc}_{\rm
NLH}=-\sigma^{a;bc}_{\rm NLH}.$ (38)

But, we emphasize that this analysis is based on the subtraction and re-
addition of the individual scattering channels to the full scattering matrix.
The effect of all scattering channels cannot be disentangled, and
$\sum_{\bm{\mathbf{q}}}\Delta^{\bm{\mathbf{q}}}\sigma$ should not be viewed as
the contribution of phonons with wavevector ${\bm{\mathbf{q}}}$ to the
scattering. Instead, it should be understood as a way to quantify how
sensitive the conductivity is to the corresponding scattering process.

Figure S4:  (a) Scattering diagnostics for the linear and NLH conductivities
of GeTe for phonons with momentum ${\bm{\mathbf{q}}}=(q_{x},0,0)$. (b) Zoom-in
of (a). $N_{q}$ is the size of the $q$-point grid, which we choose to be
identical to $N_{k}$.  Figure S5:  (a) Scattering diagnostics for the linear
and NLH conductivities of GeTe for phonons with momentum $q$ near $\Gamma$.
(b) Zoom-in of (a).

Figure S4 shows the momentum-resolved conductivity along the $q_{x}$ axis.
From Fig. S4(a), we see that the small-momentum phonon modes suppress the
conductivity the most, as in many polar semiconductors [13]. Figure S4(b)
shows the contribution of phonons with larger momentum, which are relevant for
intervalley scattering. For the linear conductivity, the momentum-resolved
contribution is always negative, indicating that the scattering dissipates the
current and enhances the resistivity. In contrast, for the NLH effect,
large-$q$ phonons with momentum slightly smaller than the intervalley distance
enhance the NLH conductivity. This result is consistent with our
interpretation based on valley polarization. Phonons with momentum smaller
than the intervalley distance increase the anisotropy of the scattering (by
increasing $\gamma_{2}$ in Fig. S6 where $\gamma_{2}>\gamma_{1}$), thus
enhancing the valley polarization and the NLH conductivity. Phonons with
momentum larger than the intervalley distance tend to decrease the anisotropy
of the scattering (by increasing $\gamma_{1}$ in Fig. S6) and thus suppress
the NLH conductivity.

Figure S5 shows the momentum-resolved conductivity for $q$ points near
$\Gamma$. For the inter-minivalley scattering, it is more difficult to analyze
the role of each scattering vector, since a single scattering vector
${\bm{\mathbf{q}}}$ contributes to many different scatterings
${\bm{\mathbf{k}}}\to{\bm{\mathbf{k+q}}}$. Nevertheless, we find that the
momentum-resolved contributions to linear conductivity and NLH conductivity
show a large difference near the inter-minivalley distance, implying the
importance of inter-minivalley scatterings on enhancing the NLH conductivity.

##  Appendix E Scattering-induced valley polarizability

In order to demonstrate that the scattering-induced valley polarization is
common in multi-valley semiconductors, we define and compute the valley
polarizability. Although here we consider only the ordinary valleys, not the
minivalleys, the same analysis can be applied to the minivalleys if they are
present. For a given valley $V$ centered around the $k$ point
${\bm{\mathbf{k}}}_{V}$, we define the valley polarization as

$\delta n_{V}=\frac{1}{V_{\rm cell}\,N_{k}}\sum_{{\bm{\mathbf{k}}}\in
V}\sum_{n}\delta\\!f_{n\mathbf{k}}.$ (39)

We define the valley polarizability ${\bm{\mathbf{\alpha}}}_{V}$ as the valley
polarization per current ${\bm{\mathbf{J}}}$:

$\delta n_{V}=\sum_{a}\alpha_{V}^{a}J^{a}\,.$ (40)

For the valleys to be well-defined, different valleys must not overlap with
each other and the nonequilibrium electron occupation $\delta\\!f$ outside the
valleys must be negligible.

In nonmagnetic semiconductors, one requirement for nonzero valley
polarizability is that the valley is not at a time-reversal invariant
momentum:

${\bm{\mathbf{\alpha}}}_{V}={\bm{\mathbf{0}}}\text{ if
}{\bm{\mathbf{k}}}_{V}=-{\bm{\mathbf{k}}}_{V}+{\bm{\mathbf{G}}}\,,$ (41)

where ${\bf G}$ is a reciprocal lattice vector. This condition is required
because the scattering-induced valley polarization transfers electrons from a
valley to its time-reversal pair. Formally, if ${\bm{\mathbf{k}}}_{V}$ is a
time-reversal invariant momentum, the corresponding valley polarization
$\delta n_{V}$ should remain unchanged under time reversal, which reverses the
sign of ${\bm{\mathbf{J}}}$. Thus, Eq. (40) gives $\delta n_{V}=-\delta
n_{V}=0$.

Also, if a valley $V$ is mapped to a valley $V^{\prime}$ (which may or may not
be indentical to $V$) by a spatial symmetry operation $s$ that maps a vector
${\bm{\mathbf{k}}}$ to $s({\bm{\mathbf{k}}})$, the valley polarizibility
should have the same symmetry:

${\bm{\mathbf{\alpha}}}_{V^{\prime}}=s({\bm{\mathbf{\alpha}}}_{V})\text{ if
}{\bm{\mathbf{k}}}_{V^{\prime}}=s({\bm{\mathbf{k}}}_{V})+{\bm{\mathbf{G}}}\,.$
(42)

The reason is that the valley polarization at $V^{\prime}$ under a current
$s({\bm{\mathbf{J}}})$ should equal the valley polarization at $V$ under a
current ${\bm{\mathbf{J}}}$:

${\bm{\mathbf{\alpha}}}_{V^{\prime}}\cdot
s({\bm{\mathbf{J}}})={\bm{\mathbf{\alpha}}}_{V}\cdot{\bm{\mathbf{J}}}\,.$ (43)

This condition forbids scattering-induced valley polarization in pristine
graphene due to the C3 symmetry [14], because no nonzero in-plane vector can
satisfy
$\mathrm{C}_{3}\,{\bm{\mathbf{\alpha}}}_{K}={\bm{\mathbf{\alpha}}}_{K}$.

Having an inversion symmetry is compatible with nonzero valley polarization if
the above conditions are satisfied. Consider an inversion-symmetric material
in which inversion maps a valley $V$ to a different valley $V^{\prime}$
satisfying
${\bm{\mathbf{k}}}_{V^{\prime}}=-{\bm{\mathbf{k}}}_{V}+{\bm{\mathbf{G}}}$.
Then, the constraint the inversion symmetry imposes on the valley polarization
is ${\bm{\mathbf{\alpha}}}_{V^{\prime}}=-{\bm{\mathbf{\alpha}}}_{V}$, and this
value can be nonzero. For an electric current ${\bm{\mathbf{J}}}$, the valley
polarizations read $\delta n_{V}=-\delta
n_{V^{\prime}}={\bm{\mathbf{\alpha}}}_{V}\cdot{\bm{\mathbf{J}}}$. If the
current direction is inverted, the valley polarizations are also inverted:
$\delta n_{V}=-\delta
n_{V^{\prime}}=-{\bm{\mathbf{\alpha}}}_{V}\cdot{\bm{\mathbf{J}}}$.

In the CRTA, the valley polarizability is always negligible because the number
of carriers within each valley is almost conserved:

$\displaystyle\delta n_{V}^{\rm CRTA}$
$\displaystyle\propto\tau\,{\bm{\mathbf{E}}}\cdot\int_{V}d{\bm{\mathbf{k}}}{\bm{\mathbf{\nabla}}}f(\varepsilon_{n\mathbf{k}})$
$\displaystyle=\tau\,{\bm{\mathbf{E}}}\cdot\int_{\partial
V}{\bm{\mathbf{dS}}}f(\varepsilon_{n\mathbf{k}})$ $\displaystyle\approx 0,$
(44)

where $\partial V$ is the surface integral for the boundary of the valley. The
surface integral vanishes because $f(\varepsilon_{n\mathbf{k}})$ is, to a good
approximation, a constant at the boundary of the valley: 1 for the hole-doped
case and 0 for the electron-doped case.

Figure S6:  A four-state model for valley polarization.

To study how a finite valley polarization can emerge, let us consider the
four-state model shown in Fig. S6. This model provides a minimal description
of the scattering-induced valley polarization. It is inspired by the model of
Ref. [14], which has two circular Fermi pockets and scattering matrix elements
that vary sinusoidally with the electron wavevector direction and was used to
describe the scattering-induced valley polarization switching of twisted
bilayer graphene.

The BTE of the four-state model reads

$\displaystyle\begin{pmatrix}\Gamma+\gamma_{2}+\gamma_{3}&-\Gamma&-\gamma_{3}&-\gamma_{2}\\\
-\Gamma&\Gamma+\gamma_{1}+\gamma_{3}&-\gamma_{1}&-\gamma_{3}\\\
-\gamma_{3}&-\gamma_{1}&\Gamma+\gamma_{1}+\gamma_{3}&-\Gamma\\\
-\gamma_{2}&-\gamma_{3}&-\Gamma&\Gamma+\gamma_{2}+\gamma_{3}\end{pmatrix}$
$\displaystyle\times$ $\displaystyle\begin{pmatrix}\delta\\!f_{\rm L2}\\\
\delta\\!f_{\rm L1}\\\ \delta\\!f_{\rm R1}\\\ \delta\\!f_{\rm
R2}\end{pmatrix}=\begin{pmatrix}c\\\ -c\\\ c\\\ -c\end{pmatrix}.$ (45)

Here, $c$ is a constant proportional to the band velocity and the external
field. Setting $c>0$ denotes a situation where the electrons move to the right
and the current flows to the left. This $4\times 4$ scattering matrix has rank
3 due to charge conservation: the sum of each column is 0. To find the
physical stationary-state solution, we need to impose charge conservation

$\delta\\!f_{\rm L1}+\delta\\!f_{\rm L2}+\delta\\!f_{\rm R1}+\delta\\!f_{\rm
R2}=0.$ (46)

Let us assume that the intravalley scattering rate $\Gamma$ is much higher
than the intervalley scattering rates $\gamma_{1}$, $\gamma_{2}$, and
$\gamma_{3}$. Then, the stationary-state solution of the BTE is

$\displaystyle\delta\\!f_{\rm L1}=-\delta\\!f_{\rm
R1}\approx-\frac{\gamma_{2}+\gamma_{3}}{\Gamma(\gamma_{1}+\gamma_{2}+2\gamma_{3})}c$
$\displaystyle\delta\\!f_{\rm L2}=-\delta\\!f_{\rm
R2}\approx\frac{\gamma_{1}+\gamma_{3}}{\Gamma(\gamma_{1}+\gamma_{2}+2\gamma_{3})}c\,.$
(47)

The corresponding valley polarization $\delta n=\delta\\!f_{\rm
L1}+\delta\\!f_{\rm L2}-\delta\\!f_{\rm R1}-\delta\\!f_{\rm R2}$ is

$\displaystyle\delta
n\approx\frac{2(\gamma_{1}-\gamma_{2})}{\Gamma(\gamma_{1}+\gamma_{2}+2\gamma_{3})}c\,.$
(48)

Thus, if the intervalley scattering rate is anisotropic, i. e. ,
$\gamma_{1}\neq\gamma_{2}$, the electric field can generate a nonzero valley
polarization.

Figure S7:  (a-e) Valley polarizability, (f-j) contribution of each relaxon to
the valley polarizability, (k-o) band structure, and (p-t) electron occupation
for five semiconductors at 300 K. For the cubic semiconductors, the valley
polarizabilities are shown for carrier densities low enough such that the
chemical potential is below the X point energy and the CBM valleys are
separated. The dashed vertical lines in (k-t) indicate the position of the
band extrema. The $k$-point paths used in (k-t) are the one-dimensional cut
shown in Fig. LABEL:fig:relaxon_dist(a) for GeTe (k, p), and the
$\mathrm{\Gamma}$-X-$\mathrm{\Gamma}$ path (a straight line in momentum space,
parallel to the electric field) with the X point at $k=0$ for other systems
(l-o, q-t).

Finally, we compute the valley polarizability from first principles for the
hole-doped GeTe and four electron-doped cubic semiconductors, diamond,
silicon, GaP, and AlSb, which have conduction-band valleys around $k$ points
that are not time-reversal invariant momenta. For the cubic semiconductors,
the valley polarizability vector is parallel to the valley momentum due to
cubic symmetry:
${\bm{\mathbf{\alpha}}}_{V}=\alpha{\bm{\mathbf{k}}}_{V}/\absolutevalue{{\bm{\mathbf{k}}}_{V}}$.
For GeTe, the valley polarizability is isotropic only within the $xy$ plane:

${\bm{\mathbf{\alpha}}}_{V}=\alpha_{\parallel}\,{\bm{\mathbf{k}}}_{V,\parallel}\,/\,\absolutevalue{{\bm{\mathbf{k}}}_{V,\parallel}}+(0,0,\alpha_{\perp})\,k_{V,z}\,/\,\absolutevalue{k_{V,z}}\,,$
(49)

where ${\bm{\mathbf{k}}}_{V}=({\bm{\mathbf{k}}}_{V,\parallel},k_{V,z})$. For
GeTe, we only compute the in-plane component $\alpha_{\parallel}$.

Figures S7(a-e) show the valley polarizability of the five materials. We find
that the valley polarizability is of the same order of magnitude in all cases,
demonstrating that valley polarization generally occurs in mutivalley
semiconductors, regardless of the presence of polar phonons or spin-orbit
coupling. From the relaxon decomposition in Figs. S7(f-j), we see that the
long-lived relaxons make a significant and often dominant contribution to the
valley polarizability. Figures S7(k-o) show the relevant band structures, and
Figs. S7(p-t) show the corresponding steady-state electron occupation obtained
by solving the BTE.

This valley polarization occurs in the linear response regime and describes a
transfer of electrons between a valley and its time-reversal pair. This
phenomenon is different from the polarization between valleys on different
crystal axes under high electric fields studied in Refs. [15, 16]. We note
that while Ref. [16] states that the six valleys of silicon are equally
occupied in the linear-response regime, our analysis shows that valley
polarization between the time-reversal pairs occurs in silicon already at the
linear order. Also, the valley polarization proposed in this work occurs for
electric fields along any direction. For example, if the electric field is
applied along the [111] direction, the three valleys along the positive [100],
[010], [001] directions will be less occupied than their inversion partners.

##  Appendix F Computational details

Figure S8:  Workflow of the calculation.

Here, we discuss the details of the calculation performed in this work. Figure
S8 summarizes the workflow of the calculation.

We used the Quantum ESPRESSO package [17] for the density functional theory
and density functional perturbation theory calculations. We used a
20$\times$20$\times$20 unshifted $k$-point grid, a kinetic energy cutoff of 80
Ry, fully relativistic norm-conserving pseudopotentials [18] from the
PseudoDojo library (v0.4) [19], and the Perdew-Burke-Ernzerhof (PBE) exchange-
correlation functional [20]. The lattice parameters and atomic coordinates
were optimized until the stresses and forces were less than $1.0\times
10^{-6}$ Ry/Bohr3 and $1.0\times 10^{-4}$ Ry/Bohr, respectively. The optimized
hexagonal lattice parameters, $a$=4.233 Å, and $c$=10.921 Å, and the Te shift
$\tau$=0.34 Å are slightly overestimated than the experimental values,
$a$=4.16 Å, $c$=10.67 Å, and $\tau$=0.31 Å [21].

We used the Wannier90 [22] and EPW [23] codes to construct the Wannier
functions and compute the real-space matrix elements. The Brillouin zone was
sampled with 12$\times$12$\times$12 and 6$\times$6$\times$6 grids for
electrons and phonons, respectively. We constructed 16 Wannier functions using
the atomic $s$ and $p$ orbitals as initial guesses. To preserve the crystal
symmetries, we performed only the disentanglement step and not the maximal
localization step. We set the upper bound of the frozen (inner) and
disentanglement (outer) windows at 3.3 eV and 7.7 eV above the valence band
maximum, respectively.

We performed the Wannier interpolation, solved the BTE, and computed the
transport coefficients using an in-house developed code ElectronPhonon.jl,
which is written in the Julia programming language [24]. Details of the
implementation will be discussed elsewhere. For the solution of the BTE, we
included only the states with energy in the range $[E_{\rm
VBM}-0.4~{}\mathrm{eV},\,E_{\rm VBM}]$, with $E_{\rm VBM}$ the valence-band
maximum energy. We sampled the Brillouin zone using unshifted
80$\times$80$\times$80 $k$\- and $q$-point grids, and sampled only the
irreducible Brillouin zone for the $k$ points utilizing the crystal and time-
reversal symmetries. To calculate the relaxons, we unfolded the scattering
matrix from the irreducible Brillouin zone to the full Brillouin zone. We
solved the linear Boltzmann transport equations using the GMRES (generalized
minimum residual) method as implemented in the IterativeSolvers.jl package
[25]. A fixed numerical Gaussian smearing of 15 meV was applied to the energy-
conserving delta functions in Eq. (A). We ignored scattering processes that
violate the total energy conservation by more than 2.5 times the smearing
parameter. Due to the higher computational cost of the relaxon calculation, we
used a narrower window $[E_{\rm VBM}-0.2~{}\mathrm{eV},\,E_{\rm VBM}]$, which
is accurate enough for carrier densities below 1020 cm-3. We used
100$\times$100$\times$100 $k$\- and $q$-point grids to obtain the results
shown in Figs. LABEL:fig:relaxon_decomp and LABEL:fig:relaxon_dist. One-
dimensional plots of the relaxon eigenvector, velocity, and Berry curvature in
Figs. LABEL:fig:relaxon_dist were obtained by cubic spline interpolation.

To calculate the band velocity, we directly Wannier interpolated the velocity
matrix elements, as done in Refs. [26, 27]. To calculate the Berry curvature,
we used Eq. (A1) of Ref. [28]. We regularized the denominator including
intermediate states to avoid divergence due to a vanishing denominator:

$\displaystyle{\bm{\mathbf{\Omega}}}_{n\mathbf{k}}(\omega)=-\sum_{m}$
$\displaystyle\real\frac{(\varepsilon_{m\mathbf{k}}-\varepsilon_{n\mathbf{k}})^{2}}{(\varepsilon_{m\mathbf{k}}-\varepsilon_{n\mathbf{k}})^{2}-(\hbar\omega+i\eta)^{2}}$
$\displaystyle\times\imaginary({\bm{\mathbf{A}}}_{nm{\bm{\mathbf{k}}}}\times{\bm{\mathbf{A}}}_{mn{\bm{\mathbf{k}}}}).$
(50)

The regularization parameter was set to $\eta=15~{}\mathrm{meV}$. The
frequency factor and the regularization parameter were only included in the
second and third terms of Eq. (A1) of Ref. [28]. The frequency dependence of
the first term can be safely neglected since it describes the coupling between
a low-energy state near the Fermi level and a high-energy state outside the
Wannierization window, so that
$\absolutevalue{\varepsilon_{m\mathbf{k}}-\varepsilon_{n\mathbf{k}}}\gg\hbar\omega,\eta$.

Doping was treated using the rigid-band approximation. The screening of the
electron-phonon coupling from the free carriers was included [29]. The
momentum- and frequency-dependent dielectric function
$\epsilon(\omega,{\bm{\mathbf{q}}})$ was calculated using the random phase
approximation [1] and the material-specific electronic structure:

$\chi_{0}(\omega,{\bm{\mathbf{q}}})=\frac{1}{V_{\rm
cell}\,N_{k}}\sum_{{\bm{\mathbf{k}}}}\sum_{n,m=1}^{N_{\rm
W}}\frac{(f_{n\mathbf{k}}-f_{m\mathbf{k}+\mathbf{q}})\absolutevalue{M_{mn}({\bm{\mathbf{k}}},{\bm{\mathbf{q}}})}^{2}}{\varepsilon_{n\mathbf{k}}-\varepsilon_{m\mathbf{k}+\mathbf{q}}+\hbar\omega+i\eta}$
(51)

and

$\epsilon(\omega,{\bm{\mathbf{q}}})=1-\frac{4\pi
e^{2}}{{\bm{\mathbf{q}}}\cdot{\bm{\mathbf{\epsilon^{\infty}}}}\cdot{\bm{\mathbf{q}}}}\chi_{0}(\omega,{\bm{\mathbf{q}}})\,,$
(52)

where

$M_{mn}({\bm{\mathbf{k}}},{\bm{\mathbf{q}}})=[U^{\dagger}({\bm{\mathbf{k+q}}})\,U({\bm{\mathbf{k}}})]_{mn}$
(53)

is the overlap of the eigenvectors at ${\bm{\mathbf{k}}}$ and
${\bm{\mathbf{k+q}}}$ in the Wannier basis, and
${\bm{\mathbf{\epsilon^{\infty}}}}$ the infinite-frequency susceptibility
tensor of the undoped system. This method is similar to Eq. (8) of Ref. [30],
with the differences that in our work (i) the frequency dependence of the
dielectric function is taken into account, and (ii) the matrix element is
approximated as
$\innerproduct{u_{m{\bm{\mathbf{k+q}}}}}{u_{n{\bm{\mathbf{k}}}}}\approx
M_{mn}({\bm{\mathbf{k}}},{\bm{\mathbf{q}}})$ instead of
$\innerproduct{u_{m{\bm{\mathbf{k+q}}}}}{u_{n{\bm{\mathbf{k}}}}}\approx\delta_{mn}$.
We used $\eta$=15 meV which was chosen based on the calculated electron
quasiparticle lifetimes. We screened only the long-range Fröhlich part, not
the Wannier-interpolated short-range part, because the free-carrier screening
dominantly acts on the long-range (small-$q$) component of the potential, and
the dielectric function $\epsilon(\omega,{\bm{\mathbf{q}}})$ calculated using
Eqs. (51-53) is intended to work only in that limit. Concretely, we used [30]

$g^{\rm scr.}_{mn\nu}({\bm{\mathbf{k}}},{\bm{\mathbf{q}}})=g^{\rm short-
range}_{mn\nu}({\bm{\mathbf{k}}},{\bm{\mathbf{q}}})+g^{\rm
Fr\ddot{o}hlich}_{mn\nu}({\bm{\mathbf{k}}},{\bm{\mathbf{q}}})/\epsilon(\omega_{\mathbf{q}\nu},{\bm{\mathbf{q}}}).$
(54)

The quadrupole contribution to the phonon dynamical matrix [31] and the
electron-phonon coupling [32, 33] was not included.

For the calculation of cubic semiconductors, diamond, silicon, GaP, and AlSb,
we used the calculation parameters used in Ref. [26] unless otherwise
specified. For the cubic semiconductors, we neglected spin-orbit coupling and
used scalar-relativistic pseudopotentials. The Wannier functions were obtained
with a 12$\times$12$\times$12 grid and a 6$\times$6$\times$6 grid for
electrons and phonons, respectively. To calculate the valley polarization with
the BTE, we sampled the $k$ points with an $N_{k}\times N_{k}\times N_{k}$
grid with $N_{k}=150$ and used 15 meV smearing for the energy-conserving delta
function. For the relaxon calculation we used $N_{k}=120$ for silicon and
AlSb, and $N_{k}=100$ for diamond and GaP after a convergence study.

The Brillouin zone plots in Fig. LABEL:fig:relaxon_dist(a, b) were generated
using Brillouin.jl [34]. Cubic spline interpolation was performed using
Interpolations.jl [35].

We note that Ref. [36] also provides a first-principles calculation of the
mobility of GeTe considering electron-phonon scattering using the momentum
relaxation time approximation. The linear mobility of GeTe reported there is
about 50% smaller than that of our calculation. A possible explanation is that
the calculation of the band velocity neglecting the contribution of the
nonlocal pseudopotential in Ref. [36] led to an underestimation of the band
velocity and mobility. It has been shown in Ref. [26] that such an
approximation can lead to a sizable relative error in the mobility, ranging
from $-$71% to $+$37%, depending on the material. This explanation is
supported by the fact that the scattering rates of Ref. [36] are in good
agreement with our results, while the mean-free paths, which are the band
velocities divided by the scattering rates, are 2 to 3 times smaller than our
result.

## References

  * Mahan [2000] G. D. Mahan,  _Many-particle physics_ , 3rd ed. (Kluwer Academic / Plenum Publishers, New York, 2000).
  * Ziman [2001] J. M. Ziman,  _Electrons and phonons: the theory of transport phenomena in solids_ (Oxford university press, Oxford, 2001).
  * Poncé  _et al._ [2020] S. Poncé, W. Li, S. Reichardt, and F. Giustino, First-principles calculations of charge carrier mobility and conductivity in bulk semiconductors and two-dimensional materials, Rep. Prog. Phys. 83, 036501 (2020).
  * Zhou  _et al._ [2015] J. Zhou, B. Liao, B. Qiu, S. Huberman, K. Esfarjani, M. S. Dresselhaus, and G. Chen, Ab initio optimization of phonon drag effect for lower-temperature thermoelectric energy conversion, PNAS 112, 14777 (2015).
  * Fiorentini and Bonini [2016] M. Fiorentini and N. Bonini, Thermoelectric coefficients of $n$-doped silicon from first principles via the solution of the Boltzmann transport equation, Phys. Rev. B 94, 085204 (2016).
  * Protik and Kozinsky [2020] N. H. Protik and B. Kozinsky, Electron-phonon drag enhancement of transport properties from a fully coupled ab initio Boltzmann formalism, Phys. Rev. B 102, 245202 (2020).
  * Xiao  _et al._ [2010] D. Xiao, M.-C. Chang, and Q. Niu, Berry phase effects on electronic properties, Rev. Mod. Phys. 82, 1959 (2010).
  * Sodemann and Fu [2015] I. Sodemann and L. Fu, Quantum nonlinear Hall effect induced by Berry curvature dipole in time-reversal invariant materials, Phys. Rev. Lett. 115, 216806 (2015).
  * Tsirkin  _et al._ [2018] S. S. Tsirkin, P. A. Puente, and I. Souza, Gyrotropic effects in trigonal tellurium studied from first principles, Phys. Rev. B 97, 035158 (2018).
  * Hardy [1970] R. J. Hardy, Phonon Boltzmann equation and second sound in solids, Phys. Rev. B 2, 1193 (1970).
  * Cepellotti and Marzari [2016] A. Cepellotti and N. Marzari, Thermal transport in crystals as a kinetic theory of relaxons, Phys. Rev. X 6, 041013 (2016).
  * Gunnarsson  _et al._ [2015] O. Gunnarsson, T. Schäfer, J. P. F. LeBlanc, E. Gull, J. Merino, G. Sangiovanni, G. Rohringer, and A. Toschi, Fluctuation diagnostics of the electron self-energy: Origin of the pseudogap physics, Phys. Rev. Lett. 114, 236402 (2015).
  * Zhou and Bernardi [2016] J.-J. Zhou and M. Bernardi,  _Ab Initio_ electron mobility and polar phonon scattering in GaAs, Phys. Rev. B 94, 201201 (2016).
  * Ying  _et al._ [2021] X. Ying, M. Ye, and L. Balents, Current switching of valley polarization in twisted bilayer graphene, Phys. Rev. B 103, 115436 (2021).
  * Isberg  _et al._ [2013] J. Isberg, M. Gabrysch, J. Hammersberg, S. Majdi, K. K. Kovi, and D. J. Twitchen, Generation, transport and detection of valley-polarized electrons in diamond, Nat. Mater. 12, 760 (2013).
  * Maliyov  _et al._ [2021] I. Maliyov, J. Park, and M. Bernardi, Ab initio electron dynamics in high electric fields: Accurate prediction of velocity-field curves, Phys. Rev. B 104, L100303 (2021).
  * Giannozzi  _et al._ [2017] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. Buongiorno Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. Dal Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. Otero-de-la-Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, Advanced capabilities for materials modelling with Quantum ESPRESSO, J. Condens. Matter Phys. 29, 465901 (2017).
  * Hamann [2013] D. R. Hamann, Optimized norm-conserving Vanderbilt pseudopotentials, Phys. Rev. B 88, 085117 (2013).
  * van Setten  _et al._ [2018] M. van Setten, M. Giantomassi, E. Bousquet, M. Verstraete, D. Hamann, X. Gonze, and G.-M. Rignanese, The PseudoDojo: Training and grading a 85 element optimized norm-conserving pseudopotential table, Comput. Phys. Commun. 226, 39 (2018).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Krempaský  _et al._ [2016] J. Krempaský, H. Volfová, S. Muff, N. Pilet, G. Landolt, M. Radović, M. Shi, D. Kriegner, V. Holý, J. Braun, H. Ebert, F. Bisti, V. A. Rogalev, V. N. Strocov, G. Springholz, J. Minár, and J. H. Dil, Disentangling bulk and surface Rashba effects in ferroelectric $\alpha$ -GeTe, Phys. Rev. B 94, 205111 (2016).
  * Pizzi  _et al._ [2020] G. Pizzi, V. Vitale, R. Arita, S. Blügel, F. Freimuth, G. Géranton, M. Gibertini, D. Gresch, C. Johnson, T. Koretsune, J. Ibañez-Azpiroz, H. Lee, J.-M. Lihm, D. Marchand, A. Marrazzo, Y. Mokrousov, J. I. Mustafa, Y. Nohara, Y. Nomura, L. Paulatto, S. Poncé, T. Ponweiser, J. Qiao, F. Thöle, S. S. Tsirkin, M. Wierzbowska, N. Marzari, D. Vanderbilt, I. Souza, A. A. Mostofi, and J. R. Yates, Wannier90 as a community code: New features and applications, J. Condens. Matter Phys. 32, 165902 (2020).
  * Poncé  _et al._ [2016] S. Poncé, E. Margine, C. Verdi, and F. Giustino, EPW: Electron–phonon coupling, transport and superconducting properties using maximally localized Wannier functions, Comput. Phys. Commun. 209, 116 (2016).
  * Bezanson  _et al._ [2017] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, Julia: A fresh approach to numerical computing, SIAM Review 59, 65 (2017).
  * Jul [2021] IterativeSolvers.jl (v0.9.2), https://github.com/JuliaLinearAlgebra/IterativeSolvers.jl (2021).
  * Poncé  _et al._ [2021] S. Poncé, F. Macheda, E. R. Margine, N. Marzari, N. Bonini, and F. Giustino, First-principles predictions of Hall and drift mobilities in semiconductors, Phys. Rev. Res. 3, 043022 (2021).
  * Lihm and Park [2022] J.-M. Lihm and C.-H. Park, Comprehensive theory of second-order spin photocurrents, Phys. Rev. B 105, 045201 (2022).
  * Wang  _et al._ [2006] X. Wang, J. R. Yates, I. Souza, and D. Vanderbilt, Ab initio calculation of the anomalous Hall conductivity by wannier interpolation, Phys. Rev. B 74, 195118 (2006).
  * Verdi  _et al._ [2017] C. Verdi, F. Caruso, and F. Giustino, Origin of the crossover from polarons to Fermi liquids in transition metal oxides, Nat. Commun. 8, 15769 (2017).
  * Macheda  _et al._ [2022] F. Macheda, P. Barone, and F. Mauri, Electron-phonon interaction and longitudinal-transverse phonon splitting in doped semiconductors, Phys. Rev. Lett. 129, 185902 (2022).
  * Royo and Stengel [2019] M. Royo and M. Stengel, First-Principles Theory of Spatial Dispersion: Dynamical Quadrupoles and Flexoelectricity, Phys. Rev. X 9, 021050 (2019).
  * Brunin  _et al._ [2020] G. Brunin, H. P. C. Miranda, M. Giantomassi, M. Royo, M. Stengel, M. J. Verstraete, X. Gonze, G.-M. Rignanese, and G. Hautier, Electron-phonon beyond Fröhlich: Dynamical quadrupoles in polar and covalent solids, Phys. Rev. Lett. 125, 136601 (2020).
  * Jhalani  _et al._ [2020] V. A. Jhalani, J.-J. Zhou, J. Park, C. E. Dreyer, and M. Bernardi, Piezoelectric electron-phonon interaction from ab initio dynamical quadrupoles: Impact on charge transport in wurtzite GaN, Phys. Rev. Lett. 125, 136602 (2020).
  * Jul [2022a] Brillouin.jl (v0.5.4), https://github.com/thchr/Brillouin.jl (2022a).
  * Jul [2022b] Interpolations.jl (v0.14.0), https://github.com/JuliaMath/Interpolations.jl (2022b).
  * Askarpour and Maassen [2019] V. Askarpour and J. Maassen, Unusual thermoelectric transport anisotropy in quasi-two-dimensional rhombohedral GeTe, Phys. Rev. B 100, 075201 (2019).