#  _Ab-initio_ predictions of spin relaxation, dephasing and diffusion in
solids

Junqing Xu  Department of Physics, Hefei University of Technology, Hefei,
Anhui, China  Department of Chemistry and Biochemistry, University of
California, Santa Cruz, CA 95064, USA  Yuan Ping<EMAIL_ADDRESS>Department
of Materials Science and Engineering, University of Wisconsin-Madison, WI,
53706, USA  Department of Physics, University of California, Santa Cruz, CA,
95064, USA

###### Abstract

Spin relaxation, dephasing and diffusion are at the heart of spin-based
information technology. Accurate theoretical approaches to simulate spin
lifetimes ($\tau_{s}$), determining how fast the spin polarization and phase
information will be lost, are important to the understandings of underlying
mechanism of these spin processes, and invaluable to search for promising
candidates of spintronic materials. Recently, we develop a first-principles
real-time density-matrix (FPDM) approach to simulate spin dynamics for general
solid-state systems. Through the complete first-principles’ descriptions of
light-matter interaction and scattering processes including electron-phonon,
electron-impurity and electron-electron scatterings with self-consistent spin-
orbit coupling, as well as _ab initio_ Landé $g$-factor, our method can
predict $\tau_{s}$ at various conditions as a function of carrier density and
temperature, under electric and magnetic fields. By employing this method, we
successfully reproduce experimental results of disparate materials and
identify the key factors affecting spin relaxation, dephasing, and diffusion
in different materials. Specifically, we predict that germanene has long
$\tau_{s}$ ($\sim$100 ns at 50 K), a giant spin lifetime anisotropy and spin-
valley locking effect under electric fields, making it advantageous for spin-
valleytronic applications. Based on our theoretical derivations and _ab
initio_ simulations, we propose a new useful electronic quantity, named spin-
flip angle $\theta^{\uparrow\downarrow}$, for the understanding of spin
relaxation through intervalley spin-flip scattering processes. Our method can
be further applied to other emerging materials and extended to simulate
exciton spin dynamics and steady-state photocurrents due to photogalvanic
effect.

##  I Introduction

In last two decades, spintronics in unconventional semiconductors and metals,
including two-dimensional materials and their heterostructures[1, 2],
topological and magnetic materials[3, 4], hybrid perovskites[5, 6] etc., have
drawn significant interests owing to its unprecedented potentials in
microelectronics and next generation low-power electronics. Spin, as a pure
quantum mechanical object, is the fundamental information carrier instead of
charge, with much less energy dissipation. Ideally one wants such information
preserved as long as possible for stable manipulation. Therefore,
understanding how spins relax and transport is of central importance in
spintronics.

Spin is an unconserved quantity in solids due to its coupling with other
quantities, such as electron orbital. Therefore, after excess spins being
generated, spin can loose its polarization (relaxation) and phase (dephasing)
due to coupling with environment. One critical parameter describing the
timescale of such processes is spin lifetime $\tau_{s}$ including T1
(relaxation) and T2 (dephasing), which is often required to be sufficiently
long for stable detection and manipulation of spin. Accurate and reliable
theoretical approaches to simulate $\tau_{s}$ are demanded for the detailed
understandings of spin dynamics and transport phenomena, and designing new
spintronics materials and devices.

Previously, methods based on model Hamiltonian with empirical parameters[7, 8,
1] were extensively employed for simulation of spin relaxation, dephasing, and
diffusion in solids. While these methods provide some mechanistic insight,
they do not serve as predictive tools for the design of new materials and may
sometimes lead to qualitatively incorrect predictions due to the use of
simplified electronic structures and interactions.

On the other hand, the existing first-principles methodology for $\tau_{s}$
has been mostly based on Fermi’s Golden rule[9, 10] considering spin-flip
transitions, which is only applicable to systems with inversion symmetry or
high spin polarization, not suitable to lots of materials with broken
inversion symmetry, promising for quantum computing and spintronics
applications [7, 1]. Other first-principles techniques like Time-Dependent
Density Functional Theory (TDDFT)[11] are challenging for crystalline systems
due to high computational cost for describing phonon relaxations that require
large supercells. More importantly, long simulation time over nanoseconds
often required by spin relaxation is a major difficulty for TDDFT, which is
only practical for tens to a few hundred fs. While spin dynamics based on
TDDFT has been recently performed for ultrafast demagnetization of magnetic
systems within tens of fs[12, 13, 14], the intrinsic time scale and supercell
limitations mentioned above remain. A recent development of spin phonon
relaxation based on spin-spin correlation function may be a promising pathway,
where its applicability for various spin relaxation and dephasing pathways
remains to be explored [15].

To overcome these challenges, we recently developed a first-principles real-
time density-matrix (FPDM) approach[16, 17, 18] to simulate spin dynamics and
pump-probe Kerr-rotation for general solid-state systems. The approach is free
from empirical parameters and is thus of great predictive power. Through the
complete first-principles descriptions of light-matter interaction and
scattering processes including electron-phonon, electron-impurity and
electron-electron scatterings with self-consistent spin-orbit coupling (SOC)
and Landé $g$-factor, our method can predict $\tau_{s}$ as a function of
temperature and carrier density, under electric and magnetic fields. The
method was applied to disparate materials, including semiconductors and
metals, with and without inversion symmetry, in good agreement with
experimental results[17, 16, 19, 20, 18, 21].

In this article, we briefly introduce the theory and implementation of our
method, present its applications and show its predictive power using a two-
dimensional Dirac material - Germanene as a showcase. Through detailed
theoretical analysis, we show how _ab-initio_ simulations improve our
understandings of spin relaxation mechanisms and are used to identify the key
quantities and/or factors to spin relaxation, dephasing, and diffusion. At the
end, we discuss how our FPDM approach can be generalized/extended to simulate
spin dynamics of excitons, instead of free carriers, and transport properties
such as photocurrent and spin currents in broken inversion systems
(photogalvanic effect).

##  II Theory

###  II.1 Density-matrix (DM) master equation

####  II.1.1 Quantum master equation

To provide a general formulation of quantum dynamics in solid-state materials,
we start from the Liouville-von Neumann equation in the Schrödinger picture,

$\displaystyle\frac{d\rho\left(t\right)}{dt}$
$\displaystyle=-\frac{i}{\text{$\hbar$}}[H,\rho\left(t\right)],$ (1)
$\displaystyle H$ $\displaystyle=H_{0}+H^{\prime},$ (2)

where $H$, $H_{0}$, and $H^{\prime}$ are total, unperturbed and perturbation
Hamiltonians respectively. In this work,

$\displaystyle H_{0}=$ $\displaystyle
H_{e}+H_{\mathrm{ph}}+H_{\mathrm{photon}},$ (3) $\displaystyle H^{\prime}=$
$\displaystyle
H_{\mathrm{e-light}}+H_{\mathrm{e-ph}}+H_{\mathrm{e-i}}+H_{\mathrm{e-e}},$ (4)

where $H_{e}$, $H_{\mathrm{ph}}$ and $H_{\mathrm{photon}}$ are unperturbed
single-particle electronic, phonon, and photon Hamiltonian respectively.
$H_{\mathrm{e-light}}$ is the light-matter interaction term.
$H_{\mathrm{e-ph}}$, $H_{\mathrm{e-i}}$ and $H_{\mathrm{e-e}}$ describe the
electron-phonon (e-ph), electron-impurity (e-i), and electron-electron (e-e)
interaction, respectively.

In practice, the many-body density matrix in Eq. 1 is reduced to one-particle
density matrix for electrons, where the environmental degree of freedom is
traced out[22], with a proper truncated BBGKY hierarchy [23, 24]. The total
rate of change of the electronic DM $\rho$ is separated into terms related to
different parts of Hamiltonian,

$\displaystyle\frac{d\rho}{dt}=$
$\displaystyle\frac{d\rho}{dt}|_{\mathrm{coh}}+\frac{d\rho}{dt}|_{\mathrm{scatt}},$
(5)

where $\left[H,\rho\right]=H\rho-\rho H$. $\frac{d\rho}{dt}|_{\mathrm{coh}}$
describes the coherent dynamics of electrons including the free-particle
dynamics, the field-induced dynamics, etc.
$\frac{d\rho}{dt}|_{\mathrm{scatt}}$ captures the scattering between electrons
and other particles. Their detailed forms are given in the subsections below.

To obtain Eq. 5 which involves only the dynamics of electrons or the
electronic subsystem, we have assumed the environmental subsystem is
characterized by a huge number of degrees of freedom and is not perturbed by
the electronic subsystem; in this work it means there is no dynamics of
phonons and photons. Formally this is the prerequisite for the Markovian
approximation. This assumption is valid when the system is not far from
initial equilibrium, e.g., the excited carrier density is not very high. The
inclusion of dynamics of the environment, e.g. phonon degrees of freedom in
the dynamics, or formally non-Markovian approximation by taking into account
the memory effect of phonon bath, has been discussed in details in Refs. 25,
22.

####  II.1.2 Coherent terms and external fields

In general, a coherent term corresponding to a single-particle electronic
Hamiltonian $H_{\mathrm{sp}}$ reads

$\displaystyle\frac{d\rho}{dt}|_{\mathrm{sp}}=$
$\displaystyle-\frac{i}{\text{$\hbar$}}\left[H_{\mathrm{sp}},\rho\right].$ (6)

In absence of external fields, the coherent dynamics contains only one term

$\displaystyle\frac{d\rho}{dt}|_{\mathrm{free}}=$
$\displaystyle-\frac{i}{\text{$\hbar$}}\left[H_{e},\rho\right].$ (7)

$H_{e}$ is the unperturbed electronic Hamiltonian, practically computed at the
mean-field level, i.e. from density functional theory (DFT). With the
electronic eigen-basis,

$\displaystyle H_{e,kmn}=$ $\displaystyle\epsilon_{kn}\delta_{mn},$ (8)

where $k$ is k-point index, $m$ and $n$ are band indices, $\epsilon$ is the
single-particle energy, and $\delta_{mn}$ is Kronecker delta function.

Therefore, we have

$\displaystyle\left(\frac{d\rho}{dt}|_{\mathrm{free}}\right)_{kmn}=$
$\displaystyle-\frac{i}{\text{$\hbar$}}\left(\epsilon_{km}-\epsilon_{kn}\right)\rho_{kmn}.$
(9)

We note that the exchange-correlation potential in electronic Hamiltonian
could be time-dependent, which will introduce an additional term in coherent
dynamics[8, 26]. We will not discuss in detail here.

The inclusion of external fields in the DM approach for solid-state materials
is not trivial and was studied by many theorists since 1950s[27, 28, 29, 30,
31, 32, 33]. Here, we consider three spatially homogeneous fields as follows.

##### 2.1 A laser field.

In this work, we approximate that the light-mater interaction
$H_{\mathrm{e-light}}$ consists of two parts - a semiclassical part
$H_{\mathrm{laser}}$ describing electrons moving in a laser field (i.e., an
electromagnetic field caused by the incident laser) and a quantum part
describing the electron-photon interaction in the vacuum without considering
the laser. The semiclassical part corresponds to light absorption and
stimulated emission under a laser field, while the quantum part corresponds to
spontaneous emission. We discuss the coherent term due to the semiclassical
part first but discuss the quantum part in the next subsection.

The Hamiltonian of a laser with frequency $\omega_{\mathrm{laser}}$ is
approximately[16]

$\displaystyle H_{\mathrm{laser},kmn}\left(\omega,t\right)=$
$\displaystyle\frac{e}{m_{e}}{\bf A}_{\mathrm{laser}}\left(t\right)\cdot{\bf
p}_{kmn}\mathrm{exp}\left(i\omega_{\mathrm{laser}}t\right)+H.C.,$ (10)

where H.C. is Hermitian conjugate. ${\bf A}_{\mathrm{laser}}\left(t\right)$ is
the amplitude. ${\bf A}_{\mathrm{laser}}\left(t\right)$ is real (complex) for
linearly (circularly) polarized light. ${\bf p}_{kmn}$ is the momentum
operator matrix element. Note that to include the contribution from the
nonlocal part of the pseudopotentials, momentum operator is computed with
commutator $[{\bf r},H_{e}]=\nabla$.

For a Gaussian pump pulse centered at time $t_{\mathrm{center}}$ with width
$\tau_{\mathrm{pump}}$,

$\displaystyle{\bf A}_{\mathrm{laser}}\left(t\right)$ $\displaystyle={\bf
A}_{\mathrm{laser}}\frac{\mathrm{1}}{\sqrt{\sqrt{\pi}\tau_{\mathrm{pump}}}}\mathrm{exp}\left[-\frac{\left(t-t_{\mathrm{center}}\right)^{2}}{2\tau_{\mathrm{pump}}^{2}}\right].$
(11)

Note that the corresponding pump power is
$I_{\mathrm{laser}}=\omega_{\mathrm{laser}}^{2}|{\bf
A}_{\mathrm{laser}}|^{2}/\left(8\pi\alpha\right)$, where $\alpha$ is fine
structure constant.

The corresponding dynamics is

$\displaystyle\frac{d\rho}{dt}|_{\mathrm{laser}}=$
$\displaystyle-\frac{i}{\text{$\hbar$}}\left[H_{\mathrm{laser}},\rho\right].$
(12)

##### 2.2 A static electric field along a non-periodic direction.

Such electric field along non-periodic direction $E_{\mathrm{np}}$ is directly
included in the DFT calculation and is modeled by a ramp or saw-like
potential. Under such treatment, $H_{e}$ is $E_{\mathrm{np}}$-dependent. As
all terms of the master equation (Eq. 5) are built on eigensystems from
$H_{e}$, they are all $E_{\mathrm{np}}$-dependent. However, if electric field
is applied along the periodic direction, one needs to properly treat the
periodic boundary condition (e.g. using the modern theory of polarization [34,
35]).

##### 2.3 A magnetic field.

We describe the effects of an external magnetic field ${\bf B}^{\mathrm{ext}}$
using Zeeman Hamiltonian

$\displaystyle H_{Z,k}\left({\bf B}^{\mathrm{ext}}\right)=$
$\displaystyle\mu_{B}{\bf B}^{\mathrm{ext}}\cdot\left({\bf L}_{k}+g_{0}{\bf
S}_{k}\right),$ (13)

where $\mu_{B}$ is Bohr magneton; $g_{0}$ is the free-electron $g$-factor;
${\bf S}$ and ${\bf L}$ are the spin and orbital angular momentum
respectively. The simulation of ${\bf L}$ is nontrivial for periodic systems.
With the Blöch basis, the orbital angular momentum reads

$\displaystyle{\bf L}_{k,mn}=$ $\displaystyle i\left\langle\frac{\partial
u_{km}}{\partial\mathrm{{\bf k}}}\right|\times\left(\widehat{H}_{e}\left({\bf
B}^{\mathrm{ext}}=0\right)-\overline{\epsilon}_{kmn}\right)\left|\frac{\partial
u_{kn}}{\partial\mathrm{{\bf k}}}\right\rangle,$ (14)
$\displaystyle\overline{\epsilon}_{kmn}=$
$\displaystyle\frac{\epsilon_{km}+\epsilon_{kn}}{2},$ (15)

where $u$ is the periodic part of the single-particle wavefunction;
$\widehat{H}_{e}\left({\bf B}^{\mathrm{ext}}=0\right)$ is the zero-field
Hamiltonian operator. Eqs. 14 \- 15 can be proven equivalent to ${\bf
L=0.5*\left({\bf r\times{\bf p-{\bf p\times{\bf r}}}}\right)}$ with ${\bf r}$
the position operator. The detailed implementation of Eq. 14 is described in
Ref. 36.

There are two ways to consider magnetic-field effects. The first way is
including $H_{Z,k}\left({\bf B}^{\mathrm{ext}}\right)$ in $H_{e}$
perturbatively (instead of self-consistently), then the new eigensystems can
be obtained by diagonalizing $H_{e}\left({\bf B}^{\mathrm{ext}}\neq 0\right)$.
The second way is including $H_{Z,k}\left({\bf B}^{\mathrm{ext}}\right)$ in
$H^{\prime}$ and the corresponding coherent dynamics is
$-\frac{i}{\text{$\hbar$}}\left[H_{Z,k}\left({\bf
B}^{\mathrm{ext}}\right),\rho\right]$. In practice, the two approaches lead to
nearly the same dynamical quantities such as lifetimes, since
$H_{Z,k}\left({\bf B}^{\mathrm{ext}}\right)$ is rather weak - e.g. Zeeman
splitting under 1 Tesla is of order 0.1 meV for many solid-state systems. In
addition, we note that we do not consider the effect of very strong magnetic
field such as the appearance of Landau level in this work.

####  II.1.3 Scattering terms

The scattering part of quantum master equation can be separated into
contributions from several scattering channels,

$\displaystyle\frac{d\rho}{dt}|_{\mathrm{scatt}}=$
$\displaystyle\sum_{c}\frac{d\rho}{dt}|_{c},$ (16)

where $c$ labels a scattering channel corresponding to an interaction
Hamiltonian $H_{c}$. Under Born-Markov approximation and neglecting
renormalization of single particle energies, in general we have[37]

$\displaystyle\frac{d\rho_{12}}{dt}|_{c}=$
$\displaystyle\frac{1}{2}\sum_{345}\left[\begin{array}[]{c}\left(I-\rho\right)_{13}P_{32,45}^{c}\rho_{45}\\\
-\left(I-\rho\right)_{45}P_{45,13}^{c,*}\rho_{32}\end{array}\right]+H.C.,$
(19)

where $P^{c}$ is the generalized scattering matrix and H.C. denotes Hermitian
conjugate. The subindex, e.g., “1”, is the combined index of k-point and band.
The weights of k points must be considered when summing over k points.

##### 3.1 Electron-phonon.

The electron-phonon (e-ph) scattering matrix is given by[37]

$\displaystyle P_{1234}^{\mathrm{e-ph}}=$
$\displaystyle\sum_{q\lambda\pm}A_{13}^{q\lambda\pm}A_{24}^{q\lambda\pm,*},$
(20) $\displaystyle A_{13}^{q\lambda\pm}=$
$\displaystyle\sqrt{\frac{2\pi}{\hbar}}g_{13}^{q\lambda\pm}\sqrt{\delta_{\sigma}^{G}\left(\omega_{13}\pm\omega_{q\lambda}\right)}\sqrt{n_{q\lambda}^{\pm}},$
(21) $\displaystyle\omega_{13}=$ $\displaystyle\epsilon_{1}-\epsilon_{3},$
(22)

where $q$ and $\lambda$ are phonon wavevector and mode, $g^{q\lambda\pm}$ is
the e-ph matrix element, resulting from the absorption ($-$) or emission ($+$)
of a phonon, $n_{q\lambda}^{\pm}=n_{q\lambda}+0.5\pm 0.5$ in terms of phonon
Bose factors $n_{q\lambda}$, and $\delta_{\sigma}^{G}$ represents an energy
conserving $\delta$-function broadened to a Gaussian of width $\sigma$. The
explicit derivation of Lindblandian dynamics and scattering matrix can be
found in Ref. [37]. This form guarantees the positive definition of density
matrix. $g^{q\lambda\pm}$ is computed fully from first-principles with self-
consistent spin-orbit coupling based on density-functional perturbation theory
and Wannier interpolation methods[38].

##### 3.2 Electron-impurity.

Similar to the e-ph scattering, we write electron-impurity scattering as

$\displaystyle P_{1234}^{\mathrm{e-i}}=$ $\displaystyle
A_{13}^{i}A_{24}^{i,*},$ (23) $\displaystyle A_{13}^{i}=$
$\displaystyle\sqrt{\frac{2\pi}{\hbar}}g_{13}^{i}\sqrt{\delta_{\sigma}^{G}\left(\omega_{13}\right)}\sqrt{n_{i}V_{\mathrm{cell}}},$
(24) $\displaystyle g_{13}^{i}=$ $\displaystyle\left\langle 1\right|\Delta
V^{i}\left|3\right\rangle,$ (25) $\displaystyle\Delta V^{i}=$ $\displaystyle
V^{i}-V^{0},$ (26)

where $n_{i}$ and $V_{\mathrm{cell}}$ are impurity density and unit cell
volume, respectively, $V^{i}$ is the potential of the impurity system and
$V^{0}$ is the potential of the pristine system. Here we assume impurities are
randomly distributed and impurity density is sufficiently low so that the
average distance between neighboring impurities is sufficiently long with
nearly no interactions among impurities.

In this work, $g^{i}$ of neutral and ionized impurities are computed
differently as follows.

For neutral impurities, $V^{i}$ is computed with SOC using a large supercell
including an impurity. To speed up the supercell convergence, we used the
potential alignment method developed in Ref. 39.

For ionized impurities, we approximate $\Delta V^{i}$ as the potential of
point charge and is simply the product of the impurity charge $Z$ and the
screened Coulomb potential[40]. Such approximate $\Delta V^{i}$ is accurate in
the long-range limit,[40] i.e., its Fourier transform $\Delta
V^{i}\left(q\right)$ is accurate when $q\rightarrow 0$, so that it is most
suitable for intravalley relaxation process due to electron-impurity
scatterings. The justification of this approximation is that the long-range
part is the most dominant in the electrostatic interaction of ionzied
impurities[40].

##### 3.3 Electron-electron interaction.

The e-e scattering matrix is given by[37]

$\displaystyle P_{1234}^{\mathrm{e-e}}=$ $\displaystyle
2\sum_{56,78}\left(I-\rho\right)_{65}\mathscr{A}_{15,37}\mathscr{A}_{26,48}^{*}\rho_{78},$
(27) $\displaystyle\mathscr{A}_{1234}=$
$\displaystyle\frac{1}{2}\left(A_{1234}-A_{1243}\right),$ (28) $\displaystyle
A_{1234}=$
$\displaystyle\frac{1}{2}\sqrt{\frac{2\pi}{\hbar}}\left[g_{1234}^{e-e}(\delta_{\sigma,1234}^{G})^{1/2}+g_{2143}^{e-e}(\delta_{\sigma,2143}^{G})^{1/2}\right],$
(29) $\displaystyle g_{1234}^{\mathrm{e-e}}=$ $\displaystyle\left\langle
1\left(r\right)\right|\left\langle
2\left(r^{\prime}\right)\right|V\left(r-r^{\prime}\right)\left|3\left(r\right)\right\rangle\left|4\left(r^{\prime}\right)\right\rangle,$
(30)

where $V\left(r-r^{\prime}\right)$ is the screened Coulomb potential and
$\delta_{\sigma,1234}^{G}=\delta_{\sigma}^{G}\left(\epsilon_{1}+\epsilon_{2}-\epsilon_{3}-\epsilon_{4}\right)$
is a Gaussian-broadened energy conservation function. The screening is
described by Random-Phase-Approximation (RPA) dielectric function. We note
that unlike the e-ph and e-i channels, $P^{\mathrm{e-e}}$ is a function of
$\rho$ and needs to be updated during time evolution of $\rho$. This is a
clear consequence of the two-particle nature of e-e scattering.
$P^{\mathrm{e-e}}$ can be written as the difference between a direct term
($P^{\mathrm{e-e},d}$) and an exchange term ($P^{\mathrm{e-e},x}$),

$\displaystyle P^{\mathrm{e-e}}=$ $\displaystyle
P^{\mathrm{e-e},d}-P^{\mathrm{e-e},x},$ (31) $\displaystyle
P^{\mathrm{e-e},d}=$
$\displaystyle\sum_{56,78}\left(I-\rho\right)_{65}A_{15,37}A_{26,48}^{*}\rho_{78},$
(32) $\displaystyle P^{\mathrm{e-e},x}=$
$\displaystyle\sum_{56,78}\left(I-\rho\right)_{65}A_{15,37}A_{26,84}^{*}\rho_{78}.$
(33)

According to Ref. 22, the direct term is expected to dominate the dynamical
scattering processes among conduction electrons or valence electrons, allowing
us to neglect the exchange term here.

Currently, we use the static Random Phase Approximation (RPA) dielectric
function for the screening without local-field effects. Given the intraband
transition is more important, we take the approximate dielectric function as:

$\displaystyle\varepsilon\left({\bf q}\right)=$
$\displaystyle\varepsilon_{s}\varepsilon^{\mathrm{intra}}\left({\bf
q}\right),$ (34)

where $\epsilon_{s}$ is the static dielectric constant well-defined for
insulators, calculated by Density Functional Perturbation Theory (DFPT)[41].
$\varepsilon_{s}$ effectively takes into account the interband contribution in
the dielectric screening. $\varepsilon^{\mathrm{intra}}\left({\bf q}\right)$
is the intraband contribution which involves only states with free carriers
and is critical for doped semiconductors (with free carriers). We computed it
at RPA without local field effect,

$\displaystyle\varepsilon^{\mathrm{intra}}\left({\bf q}\right)=$
$\displaystyle 1-V^{\mathrm{bare}}\left({\bf q}\right)\sum_{{\bf
k}mn}\left(\begin{array}[]{c}\frac{f_{{\bf k-q},m}-f_{{\bf
k}n}}{\epsilon_{{\bf k-q},m}-\epsilon_{{\bf k}n}}\times\\\ |\left\langle
u_{{\bf k-q},m}|u_{{\bf k}n}\right\rangle|^{2}\end{array}\right),$ (37)

where the sum runs over only states having free carriers, e.g., for a n-doped
semiconductor, $m$ and $n$ are conduction band indices. In the formula above,
$f$ is time-dependent non-equilibrium occupation instead of the equilibrium
one $f^{\mathrm{eq}}$. Therefore, if the optical field ${\bf
A}_{\mathrm{laser}}\left(t\right)$ is activated,
$\varepsilon^{\mathrm{intra}}\left({\bf q}\right)$ will be updated in every
time step, as $f$ will differ from $f^{\mathrm{eq}}$ and the difference
depends on the excitation density. $V^{\mathrm{bare}}\left({\bf
q}\right)=e^{2}/\left(V_{\mathrm{cell}}\varepsilon_{0}|q|^{2}\right)$ is the
bare Coulomb potential with $\varepsilon_{0}$ vacuum permittivity.

We then have the matrix elements in reciprocal space,

$\displaystyle g_{1234}^{\mathrm{e-e}}=$ $\displaystyle
V^{\mathrm{scr}}\left({\bf q}_{13}\right)\delta_{{\bf k}_{1}+{\bf k}_{2},{\bf
k}_{3}+{\bf k}_{4}}\left\langle u_{1}|u_{3}\right\rangle\left\langle
u_{2}|u_{4}\right\rangle,$ (38) $\displaystyle V^{\mathrm{scr}}\left({\bf
q}_{13}\right)=$ $\displaystyle V^{\mathrm{bare}}\left({\bf
q}_{13}\right)/\varepsilon\left({\bf q}_{13}\right),$ (39)

where $V^{\mathrm{scr}}\left({\bf q}\right)$ is the screened Coulomb potential
and ${\bf q}_{13}={\bf k}_{1}-{\bf k}_{3}$. $\left\langle
u_{1}|u_{3}\right\rangle$ is the overlap matrix element between two periodic
parts of Bloch wave functions.

##### (iv) Spontaneous emission.

As discussed above, in this work, the light-matter interaction consists of a
semiclassical part (absorption and simulated emission by a laser field) and a
quantum part. The semiclassical part has been discussed above (Sec. II.1.2).
Here we discuss the quantum part, which describes the electron-photon
interaction in the vacuum. Similar to the electron-phonon scattering, we write
a similar form for electron-photon interaction (the underlying consideration
on positive definition of density matrix is similar) [37].

$\displaystyle P_{1234}^{\mathrm{\mathrm{sp-em}}}=$
$\displaystyle\sum_{q\lambda\pm}A_{13}^{\mathrm{photon,}q\lambda\pm}A_{24}^{\mathrm{photon,}q\lambda\pm,*},$
(40) $\displaystyle A_{13}^{\mathrm{photon,}q\lambda\pm}=$
$\displaystyle\sqrt{\frac{2\pi}{\hbar}}g_{13}^{\mathrm{photon,}q\lambda\pm}\sqrt{\delta_{\sigma}^{G}\left(\omega_{13}\pm\omega_{q\lambda}^{\mathrm{photon}}\right)}$
(41) $\displaystyle\times\sqrt{n_{q\lambda}^{\mathrm{photon,}\pm}},$ (42)

where “ +” and “ -” correspond to photon emission and absorption respectively,
$\lambda$ is the photon mode, $g_{13}^{\mathrm{photon,}q\lambda\pm}$ is the
electron-photon matrix element proportional to ${\bf p}$ matrix element. The
differences from $P^{\mathrm{e-ph}}$ are: (i) As $\omega^{\mathrm{photon}}$ is
comparable to the band gap which is much greater than $k_{B}T$,
$n_{q\lambda}^{\mathrm{photon,}-}\approx 0$ and
$n_{q\lambda}^{\mathrm{photon,}+}\approx 1$. This represents that photon
absorption is not allowed in vacuum, and every photon emission process emits
one photon. (ii) As photon momentum $q$ is tiny (at long wavelength limit),
for $A_{13}^{\mathrm{photon}}$, we have ${\bf k}_{1}\approx\text{{\bf
k}}_{3}$. Therefore,

$\displaystyle A_{13}^{\mathrm{photon,}q\lambda+}=$
$\displaystyle\sqrt{\frac{2\pi}{\hbar}}g_{13}^{\mathrm{photon,}q\lambda+}\sqrt{\delta_{\sigma}^{G}\left(\omega_{13}+\omega_{q\lambda}^{\mathrm{photon}}\right)}\delta_{{\bf
k}_{1}{\bf k}_{3}}.$ (43)

The detailed form of $P_{1234}^{\mathrm{\mathrm{sp-em}}}$ will be given in our
future work.

###  II.2 Spin lifetime and spin diffusion length

####  II.2.1 Spin lifetime: relaxation and dephasing

In most spin lifetime $\tau_{s}$ experiments, $\tau_{s}$ of ensemble spins are
measured. The spin observable of an ensemble is

$\displaystyle S^{\mathrm{tot}}_{i}\left(t\right)=$
$\displaystyle\mathrm{Tr}\left[s_{i}\rho\left(t\right)\right]=\sum_{k}\sum_{mn}s_{i,kmn}\rho_{knm}\left(t\right),$
(44)

where $s_{i}$ is spin Pauli matrix in Blöch basis along direction $i$. For
spin dynamics, the total excited or excess spin observable is often more
relevant and is defined as,

$\displaystyle S^{\mathrm{ex}}_{i}\left(t\right)=$ $\displaystyle
S^{\mathrm{tot}}_{i}\left(t\right)-S^{\mathrm{eq}}_{i}\left(t\right),$ (45)

where “eq” corresponds to the final equilibrium state.

The time evolution must start at an initial state (at $t\,=\,t_{0}$) with a
net spin i.e. $\delta\rho(t_{0})=\rho(t_{0})-\rho^{\mathrm{eq}}\neq 0$ such
that $S^{\mathrm{ex}}_{i}\left(t\right)\neq 0$. We evolve $\rho\left(t\right)$
through the quantum master equation Eq. 5 for a long enough simulation time,
typically from a few ps to a few hundred ns, until the evolution of
$S^{\mathrm{ex}}_{i}\left(t\right)$ can be reliably fitted by

$\displaystyle S^{\mathrm{ex}}_{i}\left(t\right)=$ $\displaystyle
S^{\mathrm{ex}}_{i}\left(t_{0}\right)\mathrm{exp}\left[-\frac{t-t_{0}}{\tau_{s,i}}\right]$
(46)
$\displaystyle\times\mathrm{cos}\left[\omega\left(t-t_{0}\right)+\phi\right]$

to extract the relaxation time, $\tau_{s,i}$. Above, $\omega$ is oscillation
frequency due to Larmor precession.

In Ref. 16, we have shown that one can generate the initial spin imbalance by
applying a test magnetic field at $t=-\infty$, allowing the system to
equilibrate with a net spin and then turning it off suddenly at $t_{0}$, in
order to measure spin relaxation.

Historically, spin relaxation time (or longitudinal time) $T_{1}$ and ensemble
spin dephasing time (or transverse time) $T_{2}^{*}$ were used to characterize
the decay of spin ensemble[8, 42]. Suppose the spins of the system are
polarized along direction ${\bf r}_{0}$, possibly due to applying a constant
external field ${\bf B}_{0}$ along ${\bf r}_{0}$, and suppose the total
excess/excited spin is along direction ${\bf r}_{1}$, if ${\bf
r}_{1}||\text{{\bf r}}_{0}$ (or ${\bf r}_{1}\perp\text{{\bf r}}_{0}$),
$\tau_{s}$ is called $T_{1}$ (or $T_{2}^{*}$). For $T_{2}^{*}$, the excess
spin ${\bf S}^{\mathrm{ex}}\left(t\right)$ precesses with a frequency
proportional to $\left|{\bf B}_{0}\right|$.

The ensemble spin dephasing rate $1/T_{2}^{*}$ consists of reversible part and
irreversible part. The reversible part may be removed by the spin echo
technique. The irreversible part is called spin dephasing rate $1/T_{2}$,
which must be smaller than $1/T_{2}^{*}$. According to Ref. 42, $T_{2}$ may be
also defined using Eq. II.2.1 without the need of spin echo but instead of
$S^{\mathrm{tot}}_{i}\left(t\right)$, we need another quantity - the sum of
individual spin amplitudes

$\displaystyle S_{i}^{\mathrm{indiv}}=$
$\displaystyle\sum_{k}\sum_{mn}\left|s_{i,kmn}\rho_{knm}\left(t\right)\right|.$
(47)

####  II.2.2 Optical measurements

Experimentally, $\tau_{s}$ is often studied through ultrafast magneto-optical
pump-probe measurements[43, 44] where excess spin is generated by a
circularly-polarized pump pulse and its evolution is detected by probe pulses.

The probe pulse interacts with the material similarly to the pump pulse, and
could be described in exactly the same way in principle. However, this would
require repeating the simulation for several values of the pump-probe delay.
Instead, since the probe is typically chosen to be of sufficiently low
intensity, we use second-order time-dependent perturbation theory to capture
its interaction with the system [16],

$\displaystyle\Delta\rho^{\mathrm{probe}}=$
$\displaystyle\frac{1}{2}\sum_{345}\left\\{\begin{array}[]{c}\left[I-\rho\left(t\right)\right]_{13}P_{32,45}^{\mathrm{probe}}\rho\left(t\right)_{45}\\\
-\left[I-\rho\left(t\right)\right]_{45}P_{45,13}^{\mathrm{probe}}\rho\left(t\right)_{32}\end{array}\right\\}+H.C.,$
(50) $\displaystyle P_{1234}^{\mathrm{probe}}=$
$\displaystyle\sum_{\pm}A_{13}^{\mathrm{probe,\pm}}A_{24}^{\mathrm{probe,\pm,*}},$
(51) $\displaystyle A_{13}^{\mathrm{probe,\pm}}=$
$\displaystyle\sqrt{\frac{2\pi}{\hbar}}\frac{e}{m_{e}}\left({\bf
A}^{\mathrm{probe}}\cdot{\bf
p}\right)\sqrt{\delta_{\sigma}^{G}\left(\omega_{13}\pm\omega^{\mathrm{probe}}\right)}$
(52)

The change of dielectric function $\Delta\epsilon$ between the excited state
and ground state absorption detected by the probe is then

$\displaystyle\mathrm{Im}\Delta\epsilon=$
$\displaystyle\frac{2\pi}{\left(\omega^{\mathrm{probe}}\right)^{3}|{\bf
A}^{\mathrm{probe}}|^{2}}\mathrm{Tr}\left(H_{0}\Delta\rho^{\mathrm{probe}}\right).$
(53)

Note that $\Delta\rho^{\mathrm{probe}}$ contains $|{\bf
A}^{\mathrm{probe}}|^{2}$ so that $\mathrm{Im}\Delta\epsilon$ is independent
of ${\bf A}^{\mathrm{probe}}$. The $\mathrm{Im}\Delta\epsilon$ above is a
functional of $\rho$ according to Eq. 50 and is an extension of the usual
independent-particle $\mathrm{Im}\epsilon$ depending on just occupation
function $f$[45]. After obtaining $\mathrm{Im}\Delta\epsilon$, the real part
$\mathrm{Re}\Delta\epsilon$ can be obtained from the Krames-Kronig relation.

By summing up $\Delta\epsilon$ with the dielectric function for ground state
absorption, we can obtain the excited-state $\epsilon$ as inputs for Kerr and
Faraday rotation calculations[46]. These correspond to the rotations of the
polarization plane of a linearly-polarized light, reflected by (Kerr rotation)
and transmitted (Faraday rotation) through the material, after a pump
excitation with a circularly-polarized light. Time-Resolved Kerr/Faraday
Rotation (TRKR/TRFR) has been widely used to study spin dynamics of solids[47,
48]. In a TRKR experiment, a circularly-polarized pump pulse is used to excite
valence electrons to conduction bands. The transitions approximately satisfy
the selection rule of $\Delta m_{j}=\pm 1$ for left and right circularly-
polarized pulses, respectively, where $m_{j}$ is secondary total angular
momentum. TRKR works by measuring the polarization changes of reflected light,
which qualitatively is proportional to the small population imbalance of
electronic states with different $m_{j}$.

Specifically, the Kerr rotation angle $\theta_{K}$ is computed with dielectric
functions by

$\displaystyle\theta_{K}=$
$\displaystyle\mathrm{Im\frac{\sqrt{\epsilon_{+}}-\sqrt{\epsilon_{-}}}{1-\sqrt{\epsilon_{+}}\sqrt{\epsilon_{-}}}},$
(54)

where $\pm$ denotes the left and right circular polarization, respectively.

####  II.2.3 The work flow of a spin dynamics simulation

Figure 1: The work flow of a spin dynamics simulation. The codes implemented
for different steps are listed and explained in the main text.

As shown in Fig. 1, a spin dynamics simulation has three main steps:

(i) DFT step. The ground-state electronic structure, phonons, as well as the
e-ph and e-i matrix elements are firstly calculated using density functional
theory (DFT) with relatively coarse $k$ and $q$ meshes in the DFT plane-wave
code JDFTx[49]. The phonon calculation uses the finite difference method with
supercells. Alternatively, we can compute the same quantities in
QuantumEspresso, with phonon and e-ph couplings calculated by DFPT at the
coarse mesh as well.

(ii) Wannier step. We then transform all quantities from plane wave basis to
maximally localized Wannier function basis[50], and interpolate them[51, 38,
52, 53, 54, 55] to substantially finer k and q meshes. The Wannier
interpolation approach fully accounts for long-range polar terms (originally
from LO-TO splitting) in the e-ph matrix elements and phonon dispersion
relations using the approaches of Ref. 56 and 57 for the 3D and 2D systems.
This part is performed in JDFTx code.

(iii) Dynamics step. Starting from an initial state with a spin imbalance, we
evolve $\rho\left(t\right)$ through the quantum master equation of Eq. 5.
After obtaining spin observable $S\left(t\right)$ from $\rho\left(t\right)$
(Eq. 44) and fitting $S\left(t\right)$ to an exponentially oscillating decay
curve (Eq. II.2.1), the decay constant $\tau_{s}$ is obtained. This part is
performed in the DMD code interfacing with the JDFTx code.

####  II.2.4 Carrier lifetime and spin diffusion length

##### 4.1 Semiclassical limit and carrier lifetime.

At the semiclassical limit, $\rho$ is replaced by (non-equilibrium) occupation
$f$, then the scattering term originally with a full quantum description in
Eq. 19 required by DM dynamics becomes [16]:

$\displaystyle\frac{df_{1}}{dt}|_{c}=$ $\displaystyle\mathop{\sum_{2\neq
1}}\left[\left(1-f_{1}\right)P_{11,22}^{c}f_{2}-\left(1-f_{2}\right)P_{22,11}^{c}f_{1}\right],$
(55)

using the facts that $P_{11,22}$ is real and “2=1” term is zero. “c” represent
a scattering channel. Note that the weights of k points must be considered
when summing over k points.

Suppose $f$ is perturbed from its equilibrium value by $\delta f$, i.e.,
$f=f^{\mathrm{eq}}+\delta f$, then insert $f$ after perturbation into Eq. 55
and linearize it,

$\displaystyle\frac{df_{1}}{dt}|_{c}=$ $\displaystyle-\mathop{\sum_{2\neq
1}}\left[P_{11,22}^{c}f_{2}^{\mathrm{eq}}+\left(1-f_{2}^{\mathrm{eq}}\right)P_{22,11}^{c}\right]\delta
f_{1},$ (56)

using the fact that $\delta P_{11,22}$ is always zero, even for the e-e
scattering.

Define carrier/particle lifetime of state “1” $\tau_{p,1}^{c}$ by
$\frac{df_{1}}{dt}|_{c}=-\frac{\delta f_{1}}{\tau_{p,1}^{c}}$, we have

$\displaystyle\frac{1}{\tau_{p,1}^{c}}=$ $\displaystyle\mathop{\sum_{2\neq
1}}\left[P_{11,22}^{c}f_{2}^{\mathrm{eq}}+\left(1-f_{2}^{\mathrm{eq}}\right)P_{22,11}^{c}\right].$
(57)

The linewidth or the imaginary part of self-energy for the scattering channel
$c$ is related to $\tau_{p}^{c}$ by
$\mathrm{Im}\Sigma_{1}^{c}=\hbar/\left(2\tau_{p,1}^{c}\right)$.

##### 4.2 Spin diffusion length.

The DM approach has been widely employed to simulate transport properties[27,
58, 59] including spin diffusion length $l_{s}$.[8] The implementation of
simulating $l_{s}$ based on FPDM approach for general solid-state materials is
however still under development.

In this work, $l_{s}$ are obtained from first-principles using the commonly-
used relation[7, 8] based on the drift-diffusion model,

$\displaystyle l_{s}=$ $\displaystyle\sqrt{D\tau_{s}},$ (58) $\displaystyle
D=$ $\displaystyle-\mu_{c}n_{c}/\frac{dn_{c}}{d\varepsilon_{\mu_{F}}},$ (59)

where $D$ is the diffusion coefficient of carriers calculated using the
generalized Einstein relation[60]. The underlying assumption here is the
diffusion coefficient of carrier being the same as spin, which was shown to be
true in graphene-derived systems experimentally[61]. $\varepsilon_{\mu_{F}}$
is the electron chemical potential. $n_{c}$ is the carrier density. $\mu_{c}$
is the carrier mobility calculated by solving the linearized Boltzmann
equation in momentum-relaxation-time approximation[62, 63, 64],

$\displaystyle\mu_{c,i}=$
$\displaystyle\frac{e}{n_{c}V_{u}N_{k}}\sum_{1}\left[f^{\mathrm{eq}}\right]^{\prime}_{1}v_{1,i}^{2}\tau_{m,1},$
(60)

where $i=x,y,z$. $N_{k}$ is the number of k points. $V_{u}$ is the unit cell
volume. $n_{c}$ is electron density. $\left[f^{\mathrm{eq}}\right]^{\prime}$
is the derivative of the Fermi-Dirac distribution function. $v$ is the band
velocity. $\tau_{m}$ is the momentum relaxation time and is approximated
as[62, 63, 64]

$\displaystyle\tau_{m,1}^{-1}=$ $\displaystyle\mathop{\sum_{2\neq
1}}\left\\{\left[P_{11,22}f_{2}^{\mathrm{eq}}+\left(1-f_{2}^{\mathrm{eq}}\right)P_{22,11}\right]\mathrm{cos}\theta_{12}\right\\},$
(61) $\displaystyle\mathrm{cos}\theta_{12}=$ $\displaystyle\frac{{\bf
v}_{1}\cdot{\bf v}_{2}}{v_{1}v_{2}},$ (62)

where ${\bf v}$ is the velocity vector. We find that $\tau_{m}$ is similar to
the carrier lifetime $\tau_{p}$ except the angle factor
$\mathrm{cos}\theta_{12}$.

##  III Spin relaxation mechanism and analysis methods for spin relaxation

Through studying the relations between $\tau_{s}$ and important electronic
quantities, switching on/off certain dynamic processes, or tuning the key
factors affecting spin dynamics, our FPDM method is a powerful technique to
determine spin relaxation mechanism and quantitatively predict spin lifetime
in different materials at various conditions.

There are several mechanisms causing spin relaxation and dephasing of electron
carriers in non-magnetic systems[7, 8]. Among them, the most important
mechanisms are Elliot-Yafet (EY) and Dyakonov-Perel (DP) mechanisms. EY
represents spin relaxation due to spin-flip scattering. DP is activated when
inversion symmetry is broken, which results in random spin precession between
adjacent scattering events. Similar to DP mechanism, free induction decay
(FID) mechanism is caused by random spin precession but happens when
scattering is weak enough. Below we show how to determine spin relaxation
mechanism and analyse spin relaxation through _ab-initio_ simulations. We note
that our FPDM approach provides unified treatment for different mechanisms. We
also acknowledge that some approximate models were proposed in the past
specific for each mechanism, which was helpful for providing mechanistic
insights and qualitative understanding. We will introduce them in the
following sections. By comparing FPDM and these approximate models, we can
further grasp the detailed physical picture of spin relaxation.

###  III.1 Spin expectation value, spin texture and internal magnetic field

We first define some important spin-related properties in spin dynamics.

The spin expectation value along $i$ direction is defined as the diagonal
element of spin matrix $s_{i}$, i.e.,

$\displaystyle S_{i,1}^{\mathrm{exp}}=$ $\displaystyle s_{i,11}.$ (63)

Note that for degenerate bands, the matrix elements of $s_{i}$ are arbitrary,
thus, we need to diagonalize $s_{i}$ matrix in degenerate subspaces.

The spin texture presents when the spin up and down degeneracy (Kramers pair)
is lifted (e.g., due to broken inversion symmetry), and is the distribution of
the spin expectation value vector ${\bf
S}^{\mathrm{exp}}\equiv\left(S_{x}^{\mathrm{exp}},S_{y}^{\mathrm{exp}},S_{z}^{\mathrm{exp}}\right)$.

Suppose originally a system has time-reversal and inversion symmetries, so
that every two bands form a Kramers degenerate pair. Suppose the ${\bf
k}$-dependent spin matrix vectors in Bloch basis of the Kramers degenerate
pairs are ${\bf s}_{k}^{0}$ with ${\bf
s}\equiv\left(s_{x},s_{y},s_{z}\right)$. The inversion symmetry broken induces
${\bf k}$-dependent Hamiltonian terms[7]

$\displaystyle H_{k}^{\mathrm{ISB}}=$ $\displaystyle\mu_{B}g_{0}{\bf
B}_{k}^{\mathrm{in}}\cdot{\bf s}_{k}^{0},$ (64)

where ${\bf B}_{k}^{\mathrm{in}}$ is the so-called internal magnetic fields at
presence of SOC. ${\bf B}^{\mathrm{in}}$ splits the degenerate pair and
polarizes the spin along its direction. The definition of ${\bf
B}_{k}^{\mathrm{in}}$ is

$\displaystyle{\bf B}_{k}^{\mathrm{in}}\equiv$ $\displaystyle 2\Delta_{k}{\bf
S}_{k}^{\mathrm{exp}}/\left(\mu_{B}g_{0}\right),$ (65)

where $\Delta$ is the band splitting energy. From Eq. 65, ${\bf
S}_{k}^{\mathrm{exp}}||{\bf B}_{k}^{\mathrm{in}}$ (internal magnetic field
${\bf B}_{k}^{\mathrm{in}}$ is along the spin texture direction ${\bf
S}_{k}^{\mathrm{exp}}$).

###  III.2 Elliot-Yafet (EY) mechanism

EY mechanism dominates spin relaxation when spin-up and spin-down are well
defined (in absence of spin precession), so that it dominates in two types of
systems: (i) Materials with both time-reversal and spatial inversion
symmetries, e.g., silicon, free-standing graphene in absence of external
fields. With such symmetries, every two bands form a Kramers degenerate
pair[7]. Therefore, spin-up/down is well defined along an axis
${\bf\widehat{r}}$ by diagonalizing the corresponding spin matrix $s_{{\bf
r}}={\bf s}\cdot{\bf\widehat{r}}$ in degenerate subspaces. (ii) Systems with
high spin polarization e.g. due to intrinsic magnetization or spin-valley
lockng, e.g., ferromagnets, transition metal dichalcogenides (TMDs).

####  III.2.1 Fermi’s golden rule (FGR)

##### 1.1 FGR with spin-flip transitions

According to Ref. [17], if a solid-state system is close to equilibrium (but
not at equilibrium) and its spin relaxation is dominated by EY mechanism, its
free carriers’ $\tau_{s}$ due to spin-orbit coupling and e-ph scattering
approximately satisfies (for simplicity the band indices are dropped)

$\displaystyle\tau_{s}^{-1}\propto$
$\displaystyle\frac{N_{k}^{-2}}{\chi}\sum_{kq\lambda}\left\\{\begin{array}[]{c}|g_{k,k-q}^{\uparrow\downarrow,q\lambda}|^{2}n_{q\lambda}f_{k-q}^{\mathrm{eq}}\left(1-f_{k}^{\mathrm{eq}}\right)\\\
\delta\left(\epsilon_{k}-\epsilon_{k-q}-\omega_{q\lambda}\right)\end{array}\right\\},$
(68) $\displaystyle\chi=$ $\displaystyle
N_{k}^{-1}\sum_{k}f_{k}^{\mathrm{eq}}\left(1-f_{k}^{\mathrm{eq}}\right),$ (69)

where $g^{\uparrow\downarrow}$ is the spin-flip e-ph matrix element between
two electronic states of opposite spins, $n_{q,\lambda}$ is phonon occupation
at momentum $q$ and mode $\lambda$, and $f^{\mathrm{eq}}$ is Fermi-Dirac
distribution function. We will further discuss $g^{\uparrow\downarrow}$ in the
next subsection.

According to Eq. 68 and 69, $\tau_{s}^{-1}$ is proportional to
$|g_{q}^{\uparrow\downarrow}|^{2}$ and also the density of the spin-flip
transitions. Therefore we propose a temperature ($T$) and chemical potential
($\varepsilon_{\mu_{F}}$) dependent effective modulus square of the spin-flip
e-ph matrix element $\overline{|\widetilde{g}^{\uparrow\downarrow}|^{2}}$ and
a scattering density of states $D^{\mathrm{S}}$ as

$\displaystyle\overline{|\widetilde{g}^{\uparrow\downarrow}|^{2}}=$
$\displaystyle\frac{\sum_{kq}\mathrm{w}_{k,k-q}\sum_{\lambda}|g_{k,k-q}^{\uparrow\downarrow,q\lambda}|^{2}n_{q\lambda}}{\sum_{kq}\mathrm{w}_{k,k-q}},$
(70) $\displaystyle D^{\mathrm{S}}=$
$\displaystyle\frac{N_{k}^{-2}\sum_{kq}\mathrm{w}_{k,k-q}}{N_{k}^{-1}\sum_{k}f_{k}^{\mathrm{eq}}\left(1-f_{k}^{\mathrm{eq}}\right)},$
(71) $\displaystyle\mathrm{w}_{k,k-q}=$ $\displaystyle
f_{k-q}^{\mathrm{eq}}\left(1-f_{k}^{\mathrm{eq}}\right)\delta\left(\epsilon_{k}-\epsilon_{k-q}-\omega_{c}\right),$
(72)

where $\omega_{c}$ is the characteristic phonon energy chosen as the averaged
energy of the phonons contributing to spin relaxation, and
$\mathrm{w}_{k,k-q}$ is the weight function. The matrix element modulus square
is weighted by $n_{q\lambda}$ according to Eq. 68 and 69. This rules out high-
frequency phonons at low $T$ which are not excited. $\mathrm{w}_{k,k-q}$
selects transitions between states separated by $\omega_{c}$ and around the
band edge or $\varepsilon_{\mu_{F}}$, which are “more relevant” transitions to
spin relaxation.

$D^{\mathrm{S}}$ can be regarded as an effective density of spin-flip e-ph
transitions satisfying energy conservation between one state and its pairs.
When $\omega_{c}=0$, we have $D^{\mathrm{S}}=\int
d\epsilon\left(-\frac{df^{\mathrm{eq}}}{d\epsilon}\right)D^{2}\left(\epsilon\right)/\int
d\epsilon\left(-\frac{df^{\mathrm{eq}}}{d\epsilon}\right)D\left(\epsilon\right)$
with $D\left(\epsilon\right)$ density of electronic states (DOS). So
$D^{\mathrm{S}}$ can be roughly regarded as a weighted-averaged DOS with
weight
$\left(-\frac{df^{\mathrm{eq}}}{d\epsilon}\right)D\left(\epsilon\right)$.

With $\overline{|\widetilde{g}^{\uparrow\downarrow}|^{2}}$ and
$D^{\mathrm{S}}$, we have the approximate relation for spin relaxation rate,

$\displaystyle\tau_{s}^{-1}\propto$
$\displaystyle\overline{|\widetilde{g}^{\uparrow\downarrow}|^{2}}D^{\mathrm{S}}.$
(73)

We can similarly define an effective spin conserving matrix element
$\overline{|\widetilde{g}^{\uparrow\uparrow}|^{2}}$ by replacing
$g_{k,k-q}^{\uparrow\downarrow,q\lambda}$ to
$g_{k,k-q}^{\uparrow\uparrow,q\lambda}$ in Eq. 70. Then we have the
approximate relation for carrier relaxation rate due to e-ph scattering,

$\displaystyle\left\langle\tau_{p}^{-1}\right\rangle$
$\displaystyle\propto\overline{|\widetilde{g}^{\uparrow\uparrow}|^{2}}D^{\mathrm{S}}.$
(74)

Although such formula are not as general and accurate as our FPDM formulation,
Eq. 73 and 74 can be used to understand the spin flip and spin conserving e-ph
matrix element’s relation to spin or carrier relaxation, and one can develop
intuition on what type of phonons contribute more to
$\overline{|\widetilde{g}^{\uparrow\uparrow}|^{2}}$ or
$\overline{|\widetilde{g}^{\uparrow\downarrow}|^{2}}$. As our recent studies
show Frölich LO phonon strongly contributes to carrier relaxation in CsPbBr3,
but much less important in its spin relaxation due to the spin-conserving
nature of long-ranged Frölich electron-phonon coupling [18].

##### 1.2 Generalized FGR

In centrosymmetric systems with strong spin-mixing and band degeneracy, ${\bf
S}_{k}^{\mathrm{exp}}$ may have multiple values. For example, valence bands of
silicon at $\Gamma$ are four-fold degenerate and their ${\bf
S}_{k}^{\mathrm{exp}}$ are (more precisely, very close to) $\pm 1/2$ and $\pm
1/6$. In such cases, spin relaxation may be still driven by EY mechanism but
needs a generalized FGR formula[17] beyond spin-flip transition,

$\displaystyle\tau_{s,i}^{-1}\propto$
$\displaystyle\frac{N_{k}^{-2}}{\chi_{i}}\sum_{12\lambda}\left\\{\begin{array}[]{c}\left|\Delta
S_{i,12}^{\mathrm{exp}}g_{12}^{q,\lambda}\right|^{2}n_{q\lambda}\\\
f_{2}^{\mathrm{eq}}\left(1-f_{1}^{\mathrm{eq}}\right)\delta\left(\epsilon_{1}-\epsilon_{2}-\omega_{q\lambda}\right)\end{array}\right\\},$
(77) $\displaystyle\chi_{i}=$ $\displaystyle
N_{k}^{-1}\sum_{1}f_{1}^{\mathrm{eq}}\left(1-f_{1}^{\mathrm{eq}}\right)S_{i,1}^{\mathrm{exp},2},$
(78) $\displaystyle\Delta S_{i,12}^{\mathrm{exp}}=$ $\displaystyle
S_{i,1}^{\mathrm{exp}}-S_{i,2}^{\mathrm{exp}},$ (79)

where $\Delta S^{\mathrm{exp}}$ is the change of spin expectation value for a
pair of electronic states. Therefore, the transitions with non-negligible
$\Delta S^{\mathrm{exp}}$ can contribute to spin relaxation. Again, unlike our
FPDM method, this is approximate formula, which can be used to understand how
e-ph scattering between a pair of states changes spin expectation value as we
showcase bulk Si hole spin relaxation in Ref. 17.

####  III.2.2 Spin-mixing parameter $b^{2}$

##### 2.1 State-resolved $b^{2}$.

Suppose the spin of a state “1” is highly polarized along $i$ direction. Then
in general, the wavefunction of state “1” can be written as
$\Psi_{1}\left({\bf r}\right)=a_{i,1}\left({\bf
r}\right)\alpha+b_{i,1}\left({\bf r}\right)\beta$, where $a$ and $b$ are the
coefficients of the large and small components of the wavefunction, and
$\alpha$ and $\beta$ are spinors (one up and one down for direction $i$).
Define $a_{i,1}^{2}=\int|a_{i,1}\left({\bf r}\right)|^{2}d{\bf r}$ and
$b_{i,1}^{2}=\int|b_{i,1}\left({\bf r}\right)|^{2}d{\bf r}$, then
$a_{i,1}^{2}>b_{i,1}^{2}$ and $b_{i,1}^{2}$ is just spin-mixing parameter of
state “1” along direction $i$ [19].

$\displaystyle a_{i,1}^{2}+b_{i,1}^{2}=$ $\displaystyle 1,$ (80)
$\displaystyle 0.5\left(a_{i,1}^{2}-b_{i,1}^{2}\right)=$ $\displaystyle
S_{i,1}^{\mathrm{exp}},$ (81)

Therefore,

$\displaystyle b_{i,1}^{2}=$ $\displaystyle 0.5-S_{i,1}^{\mathrm{exp}}.$ (82)

##### 2.2 EY relation.

According to Eq. 68 and 74, we have

$\displaystyle\tau_{s}^{-1}/\tau_{p}^{-1}\propto$
$\displaystyle\overline{|\widetilde{g}^{\uparrow\downarrow}|^{2}}/\overline{|\widetilde{g}^{\uparrow\uparrow}|^{2}},$
(83)

where $\tau_{p}=1/\left\langle\tau_{p}^{-1}\right\rangle$.

As thermal averaging frequently appears in spin relaxation analysis, we define
$\left\langle A\right\rangle$ as the thermal average of electronic quantity
$A$,

$\displaystyle\left\langle A\right\rangle=$
$\displaystyle\frac{\sum_{kn}\left(-\left[f^{\mathrm{eq}}\right]^{\prime}_{kn}\right)A_{kn}}{\sum_{kn}\left(-\left[f^{\mathrm{eq}}\right]^{\prime}_{kn}\right)},$
(84)

where $\left[f^{\mathrm{eq}}\right]^{\prime}$ is the derivative of the Fermi-
Dirac distribution function.

According to Refs. 65, 66, 9,
$\overline{|\widetilde{g}^{\uparrow\downarrow}|^{2}}/\overline{|\widetilde{g}^{\uparrow\uparrow}|^{2}}\sim\left\langle
b^{2}\right\rangle$, so that

$\displaystyle\tau_{s}^{-1}/\tau_{p}^{-1}\sim$ $\displaystyle\left\langle
b^{2}\right\rangle.$ (85)

The above is a rough relation and cannot be used to predict $\tau_{s}$ (the
error may be several or even ten times). Practically, we use the following
approximate relation

$\displaystyle\tau_{s}^{-1}/\tau_{p}^{-1}=$ $\displaystyle 4\left\langle
b^{2}\right\rangle.$ (86)

This is called EY relation in this work.

####  III.2.3 Spin-flip angle $\theta^{\uparrow\downarrow}$

In our previous paper[21], we proposed an important electronic quantity for
intervalley spin-flip scattering - the spin-flip angle
$\theta^{\uparrow\downarrow}$ between two electronic states. For two states
$\left(k,n\right)$ and $\left(k^{\prime},n^{\prime}\right)$ with opposite spin
directions, $\theta^{\uparrow\downarrow}$ is the angle between $-{\bf
S}_{kn}^{\mathrm{exp}}$ and ${\bf S}_{k^{\prime}n^{\prime}}^{\mathrm{exp}}$.
The motivation of examining $\theta^{\uparrow\downarrow}$ is that: According
to Ref. [67], due to time-reversal symmetry, the spin-flip matrix element of
the same band between ${\bf k}$ and $-{\bf k}$ is exactly zero, so that
$g^{\uparrow\downarrow}$ is zero at lowest order for intervalley transitions
between two opposite valleys (e.g., ${\bf K}$ and ${\bf-K}$). In the first-
order perturbation level, $\left|g^{\uparrow\downarrow}\right|$ between two
states is determined by $\theta^{\uparrow\downarrow}$ between these two states
and proportional to
$\left|\sin\left(\theta^{\uparrow\downarrow}/2\right)\right|$.

Suppose (i) the inversion symmetry broken induces ${\bf B}_{k}^{\mathrm{in}}$
(Eq. 65) for a Kramers degenerate pair; (ii) there are two valleys centered at
wavevectors ${\bf Q}$ and $-{\bf Q}$ and (iii) there are two wavevectors ${\bf
k}_{1}$ and ${\bf k}_{2}$ near ${\bf Q}$ and $-{\bf Q}$ respectively. Due to
time-reversal symmetry, the directions of ${\bf B}_{k_{1}}^{\mathrm{in}}$ and
${\bf B}_{k_{2}}^{\mathrm{in}}$ are almost opposite.

We can prove that for a general operator $\widehat{A}$,

$\displaystyle\left|A_{k_{1}k_{2}}^{\uparrow\downarrow}\right|^{2}\approx$
$\displaystyle\mathrm{sin}^{2}\left(\theta_{k_{1}k_{2}}^{\uparrow\downarrow}/2\right)\left|A_{k_{1}k_{2}}^{\downarrow\downarrow}\right|^{2},$
(87)

where $A_{k_{1}k_{2}}^{\uparrow\downarrow}$ and
$A_{k_{1}k_{2}}^{\downarrow\downarrow}$ are the spin-flip and spin-conserving
matrix elements between ${\bf k}_{1}$ and ${\bf k}_{2}$ respectively. We
present the detailed derivation in SI Sec. SIII of Ref. 21.

From Eq. 87, for the intervalley e-ph matrix elements, we have

$\displaystyle\left|g_{k_{1}k_{2}}^{\uparrow\downarrow}\right|^{2}\approx$
$\displaystyle\mathrm{sin}^{2}\left(\theta_{k_{1}k_{2}}^{\uparrow\downarrow}/2\right)\left|g_{k_{1}k_{2}}^{\downarrow\downarrow}\right|^{2}.$
(88)

Finally, similar to Eq. 70, we propose an effective modulus square of
$\mathrm{sin}^{2}\left(\theta_{k_{1}k_{2}}^{\uparrow\downarrow}/2\right)$,

$\displaystyle\overline{\mathrm{sin}^{2}\left(\theta^{\uparrow\downarrow}/2\right)}=$
$\displaystyle\frac{\sum_{kq}\mathrm{w}_{k,k-q}\mathrm{sin}^{2}\left(\theta_{k,k-q}^{\uparrow\downarrow}/2\right)}{\sum_{kq}\mathrm{w}_{k,k-q}}.$
(89)

We found spin relaxation time (obtained from our FPDM method) linearly
proportional to this quantity
($\overline{\mathrm{sin}^{2}\left(\theta^{\uparrow\downarrow}/2\right)}$) when
intervalley spin relaxation process dominates, as an example of substrate
effects on spin relaxation of strong SOC Dirac materials in Ref. 21. Such
angle relates to the spin expectation value direction between initial
$\uparrow$ and final $\downarrow$ states.

###  III.3 Dyakonov-Perel (DP) and Free induction decay (FID) mechanisms

####  III.3.1 Model relations

For non-magnetic materials, with nonzero k-dependent internal magnetic field
${\bf B}_{k}^{\mathrm{in}}$ induced by inversion symmetry broken and spin-
orbit coupling, the spins at ${\bf k}$ precess about ${\bf
B}_{k}^{\mathrm{in}}$. The Larmor precession frequency vector can be defined
as:

$\displaystyle{\bf\Omega}_{k}=$ $\displaystyle\Delta_{k}\widehat{{\bf
S}}_{k}^{\mathrm{exp}},$ (90)

where $\widehat{{\bf S}}_{k}^{\mathrm{exp}}$ is the normalized ${\bf
S}_{k}^{\mathrm{exp}}$.

We define ${\bf\Omega}_{\perp\widehat{{\bf r}}}$ as the component of
${\bf\Omega}$ perpendicular to direction $\widehat{{\bf r}}$. Suppose the
fluctuation amplitude among different k-points of
${\bf\Omega}_{\perp\widehat{{\bf r}}}$ is
$\Delta{\bf\Omega}_{\perp\widehat{{\bf r}}}$ and numerically we define it as
(using Eq. 84)

$\displaystyle\Delta{\bf\Omega}_{\perp\widehat{{\bf r}}}=$
$\displaystyle\sqrt{\left\langle\left|{\bf\Omega}_{\perp\widehat{{\bf
r}}}-\left\langle{\bf\Omega}_{\perp\widehat{{\bf
r}}}\right\rangle\right|^{2}\right\rangle}.$ (91)

According to Refs. 7 and 8, a nonzero $\Delta{\bf\Omega}$ leads to finite spin
lifetime $\tau_{s}^{\Delta{\bf\Omega}}$ along ${\bf\widehat{r}}$ and the spin
relaxation mechanism depends on the magnitude of
$\tau_{p}\Delta{\bf\Omega}$[8, 7] (the subindex “$\perp\widehat{{\bf r}}$” is
dropped for simplicity):

(i) DP mechanism if $\tau_{p}\Delta{\bf\Omega}\ll 1$ (strong scattering
limit). We have the DP relation

$\displaystyle\left(\tau_{s}^{\Delta\Omega}\right)^{-1}\sim$
$\displaystyle\left(\tau_{s}^{\mathrm{DP}}\right)^{-1}\sim\tau_{p}\left(\Delta{\bf\Omega}\right)^{2}.$
(92)

(ii) FID mechanism if $\tau_{p}\Delta\Omega\gtrsim 1$ (weak scattering limit).
We have

$\displaystyle\left(\tau_{s}^{\Delta\Omega}\right)^{-1}\sim$
$\displaystyle\left(\tau_{s}^{\mathrm{FID}}\right)^{-1}\sim\Delta{\bf\Omega}.$
(93)

(iii) Between (i) and (ii) regimes, there isn’t a good approximate relation
for $\left(\tau_{s}^{\Delta\Omega}\right)^{-1}$, but we may expect that[8]

$\left(\tau_{s}^{\mathrm{DP}}\right)^{-1}<\left(\tau_{s}^{\Delta\Omega}\right)^{-1}<\left(\tau_{s}^{\mathrm{FID}}\right)^{-1}.$
(94)

####  III.3.2 Landé g-factor of free carriers

Besides intrinsic spin-orbit coupling, one channel to induce nonzero
$\Delta{\bf\Omega}$ is through $g$ factor fluctuations under magnetic field.
In Sec. II.1.2, we have described how ${\bf B}^{\mathrm{ext}}$ is considered
in FPDM approach, therefore, with $H_{z}\left({\bf B}^{\mathrm{ext}}\right)$,
the magnetic-field effects on spin dynamics are straightforwardly included in
FPDM simulations.

For the description of the magnetic-field effects, Landé $g$-factor has been
widely used[18]. For a single band and a pair of bands, $g$-factor is well
defined and relates to ${\bf B}^{\mathrm{ext}}$-induced change of energy,
energy splitting and/or Larmor precession frequency. In this work, we limit
our discussions to two Kramers degenerate bands under a transverse ${\bf
B}^{\mathrm{ext}}$ (perpendicular to the direction of the excess/excited spin
${\bf S}\equiv\left(S_{x},S_{y},S_{z}\right)$). In general, in the two-band
case, $g$-factor is a tensor and may be defined by the relation
${\bf\Omega}_{k}\left({\bf B}^{\mathrm{ext}}\right)=\mu_{B}{\bf
B}^{\mathrm{ext}}g_{k}^{S}$. Here, for simplicity, we assume
${\bf\Omega}_{k}\left({\bf B}^{\mathrm{ext}}\right)||{\bf B}^{\mathrm{ext}}$,
which is valid in many materials including CsPbBr3. With this assumption and
from Eq. 90,

$\displaystyle g_{k}^{S}=$ $\displaystyle\frac{{\bf\Omega}_{k}\left({\bf
B}^{\mathrm{ext}}\right)\cdot\widehat{{\bf
B}}^{\mathrm{ext}}}{\mu_{B}B^{\mathrm{ext}}}.$ (95)

${\bf\Omega}_{k}\left({\bf B}^{\mathrm{ext}}\right)$ is computed using Eq. 90.

However, in many previous theoretical studies[68, 69], $g$-factors were
defined based on pseudo-spins related to the total magnetic momenta
$J^{\mathrm{at}}$, which are determined from the atomic-orbital models. The
pseudo-spins can have opposite directions to the actual spins ($S$). Most
previous experimental studies adopted the same convention for the signs of
carrier $g$-factors. Therefore, to compare with $g$-factors obtained in those
previous studies, we introduce a correction factor $C^{\mathrm{S\rightarrow
J}}$ and define a new $g$-factor:

$\displaystyle\widetilde{g}_{k}\left(\widehat{{\bf B}}^{\mathrm{ext}}\right)=$
$\displaystyle C^{S\rightarrow J}g_{k}^{S}.$ (96)

$C^{\mathrm{S\rightarrow J}}=m_{S}^{\mathrm{at}}/m_{J}^{\mathrm{at}}$ with
$m_{J}^{\mathrm{at}}$ and $m_{S}^{\mathrm{at}}$ the total and spin magnetic
momenta respectively, obtained from the atomic-orbital model[69].
$C^{\mathrm{S\rightarrow J}}$ is independent from k-point, and is $\mp$1 for
electrons and holes respectively for CsPbBr3 [18].

As $\widetilde{g}_{k}$ is different at different ${\bf k}$, we can define its
fluctuation amplitude as

$\displaystyle\Delta\widetilde{g}=$
$\displaystyle\left\langle\left(\widetilde{g}-\left\langle\widetilde{g}\right\rangle\right)^{2}\right\rangle.$
(97)

From the above equation, we have the fluctuation amplitude of
${\bf\Omega}_{k}\left({\bf B}^{\mathrm{ext}}\right)$

$\displaystyle\Delta\Omega\left({\bf B}^{\mathrm{ext}}\right)=$
$\displaystyle\mu_{B}B^{\mathrm{ext}}\Delta\widetilde{g}.$ (98)

A nonzero $\Delta\Omega\left({\bf B}^{\mathrm{ext}}\right)$ leads to spin
dephasing under external magnetic field ${\bf B}^{\mathrm{ext}}$ and the
mechanism may be DP or FID depending on the magnitude of
$\tau_{p}\Delta{\bf\Omega}$ as discussed above.

###  III.4 Determination of spin relaxation mechanism

In general, the applicability of our FPDM approach does not depend on specific
spin relaxation mechanism. To determine the dominant relaxation or dephasing
mechanism, we have the following two approaches by utilizing our FPDM
calculations.

####  III.4.1 Method 1: Comparing FPDM calculations with FGR or model
relations

The first approach is to compare FPDM results directly with simple models
designed for various mechanisms as introduced in previous two sections. If
FPDM agrees with one of the model relations, it is a good indicator of
dominant mechanism.

##### 1.1 EY mechanism.

Since spin precession is suppressed when EY mechanism dominates, the DM master
equation in semiclassical limit from Eq. 55 and FGR formula with spin-flip
scattering Eq. 68 should describe well spin relaxation (suppose spin matrix is
diagonalized in degenerate subspaces along the direction of the excess spin).
Therefore, if the values and trends of $\tau_{s}$ by Eq. 55 and 68 are similar
to those obtained from FPDM calculations, the dominating mechanism is likely
the EY.

##### 1.2 DP or FID mechanism.

With nonzero $\Delta{\bf\Omega}$, whether spin relaxation is dominated by DP
or FID mechanism may be determined by comparing the values and trends of
$\tau_{s}$ by Eq. 92 or 93 with the FPDM results. We have shown the success of
such analysis in our previous work on spin relaxation and transport in
germanene and silicene under electric field [19] as well as CsPbBr3 under
magnetic field [18].

####  III.4.2 Method 2: Tuning the scattering strength or precession frequency
in FPDM calculations

Conventionally, spin relaxation mechanism is determined from the relation
between $\tau_{s}$ and $\tau_{p}$ \- EY mechanism leads to
$\tau_{s}\propto\tau_{p}$ while DP mechanism leads to
$\tau_{s}\propto\tau_{p}^{-1}$. Another way to look at this is through tuning
the strength of scatterings (e.g. physically, increasing impurity
concentration increases scattering). If $\tau_{s}$ decreases with increasing
scattering, it’s likely EY; otherwise, it’s likely DP mechanism, because
$\tau_{p}$ always decreases with increasing scatterings.

Practically in our FPDM calculations, this approach can be implemented by
introducing a scaling factor $F^{\mathrm{sc}}$ to tune the scattering
strength, i.e. multiplying $F^{\mathrm{sc}}$ to the scattering term
$\frac{d\rho}{dt}|_{\mathrm{scatt}}$ of the DM master equation (Eq. 5) but
keep the coherent term (the first term of Eq. 5) unchanged. This is equivalent
to multiply $F^{\mathrm{sc}}$ to all elements of the generalized scattering-
rate matrix $P$ (Eq. 20, 23 and 27).

To avoid confusions , we name carrier and spin lifetimes after introducing
$F^{\mathrm{sc}}$ as $\tau_{s}^{\prime}$ and $\tau_{p}^{\prime}$ respectively.
According to Eq. 57, we always have
$\tau_{p}^{\prime}/\tau_{p}=\left(F^{\mathrm{sc}}\right)^{-1}$. On the other
hand, $\tau_{s}^{\prime}$ depends on the spin relaxation mechanism. For EY
mechanism,
$\tau_{s}^{\prime}/\tau_{p}^{\prime}=\tau_{s}/\tau_{p}=C^{\mathrm{EY}}$, where
$C^{\mathrm{EY}}$ is a constant unrelated to $F^{\mathrm{sc}}$. Therefore, it
can be proven that
$\tau_{s}^{\prime}/\tau_{s}=\left(F^{\mathrm{sc}}\right)^{-1}=\tau_{p}^{\prime}/\tau_{p}$.
Similarly, for DP mechanism, we can prove that
$\tau_{s}^{\prime}/\tau_{s}=F^{\mathrm{sc}}=\left(\tau_{p}^{\prime}/\tau_{p}\right)^{-1}$.
For FID mechanism, as $\tau_{s}^{\prime}$ is irrelevant to
$\tau_{p}^{\prime}$, we have $\tau_{s}^{\prime}/\tau_{s}\equiv 1$. Therefore,
the relation between $\tau_{s}^{\prime}/\tau_{s}$ and
$\tau_{p}^{\prime}/\tau_{p}$ by tuning $F^{\mathrm{sc}}$ is useful to
understand spin relaxation mechanism.

Moreover, the relation between $\tau_{s}$ and $\Delta\Omega$ is also useful to
understand spin relaxation mechanism. $\Delta\Omega$ can be tuned easily by
tuning the energy splitting.

##  IV Applications

###  IV.1 Non-magnetic materials with inversion symmetry

Figure 2: Spin relaxation in Silicon (Si). Spin ($\tau_{s}$) and carrier
($\tau_{p}$) lifetimes of (a) electrons and (b) holes of Si. Exp. A and B are
from Ref. 70 and 71. (c) Cumulative contributions to spin relaxation by change
in spin, $\Delta$s, per scattering event defined based on Eqs. 77-79:
electrons in Si exhibit spin flips with all contributions at $\Delta$s = 1,
whereas holes in Si and electrons in iron exhibit a broad distribution in
$\Delta$s [17].

We first present results for systems with inversion symmetry traditionally
described by EY mechanism. Figure 2(a) shows that our theoretical electron
$\tau_{s}$ of Si as a function of $T$ are in excellent agreement with
experimental measurements[70, 71]. Note that previous first-principles
calculations[9] approximated spin-flip electron-phonon matrix elements from
pseudospin wavefunction overlap and spin-conserving electron-phonon matrix
element, effectively assuming that the scattering potential varies slowly on
the scale of a unit cell; we make no such approximation in our FPDM approach.
In contrast, holes in Si exhibit strong spin mixing with spin-2/3 character
and spin expectation values no longer close to $\hbar/2$. Figure 2(b) shows
our predictions for the hole spin lifetime $\tau_{s,h}$ which is much shorter
than the electron spin as a result of the strong mixing (450 fs for holes
compared to 7 ns for electrons at 300 K) and is much closer to carrier
relaxation time $\tau_{p}$. Additionally, Fig. 2(c) shows that the change in
spin expectation values ($\Delta S^{\mathrm{exp}}$) per scattering event
(evaluated using Eqs. 77 -79) has a broad distribution for holes in Si,
indicating that they cannot be described purely by spin-flip transitions,
while conduction electrons in Si predominantly exhibit spin-flip transitions
with $\Delta S^{\mathrm{exp}}=1$.

###  IV.2 Materials with high spin polarization

In ferromagnets and antiferromagnets, spins are highly polarized, so that spin
relaxation is likely dominated by EY mechanism. In Ref. 17, we simulated
$\tau_{s}$ of bcc iron, where we found good agreement with experiments, with a
dominant EY mechanism (magnon was not considered here). In nonmagnetic
materials, when inversion symmetry is broken and SOC strength is moderate or
large, it has been found that some of them such as transition metal
dichalcogenides also have highly polarized carrier spins and exhibit EY spin
relaxation.

####  IV.2.1 Transition metal dichalcogenides (TMDs)

Monolayer TMDs exhibit exciting features including spin-polarized bands,
valley-specific optical selection rules and spin-valley locking. In Ref. 72,
73, it has been shown that by introducing doping in monolayer TMDs, ultraslow
decays of Kerr rotations, which correspond to ultralong spin/valley lifetimes
of resident carriers especially resident holes, can be observed at low
temperatures. Those features establish TMDs’ advantages for spin-valleytronics
and (quantum) information processing. Although spin/valley relaxation time of
resident carriers in monolayer TMDs, most relevant to spinvalleytronics’
application, has been extensively examined, the underlying relaxation
mechanisms especially the effects of different types of impurities have not
been addressed through _ab-initio_ simulations.

In Ref. 16, we clarified the above problem by conducting _ab initio_ real-time
dynamics simulations with relevant scattering mechanisms. We focused on WSe2
due to its larger valence band SOC splitting and weaker interlayer coupling
compared with other TMDs and focused on dynamics of resident holes as
$\tau_{s}$ of holes seem longer than electrons.

For holes of monolayer WSe2, spin/valley relaxation is completely determined
by intervalley spin-flip scattering between $K$ and $K^{\prime}$ valleys
because of spin-valley locking. Previously, we reported spin/valley lifetimes
of resident holes of monolayer TMDs at $T\geq$50 K with e-ph scattering[17].
At very low temperatures, e.g., 10 K, intervalley e-ph scattering is however
not activated as the corresponding phonon occupation is negligible, so that
other scattering mechanisms need to be included. Note that e-e scattering
should not play an important role in spin relaxation of holes of TMDs. The
reason is: The e-e scattering is a two-particle process where a transition is
accompanied by another transition with energy and momentum being conserved.
Considering the fact that only the highest occupied band is involved in
dynamics of TMD holes, for a e-e process, a $K$$\rightarrow$$K^{\prime}$
($K^{\prime}$$\rightarrow$$K$) spin-flip transition must be accompanied by an
opposite $K^{\prime}$$\rightarrow$$K$ ($K$$\rightarrow$$K^{\prime}$) spin-flip
transition. Overall, an e-e scattering process cannot change the spin of the
system. Therefore, we will include only e-ph and e-i scatterings.

Figure 3: (a) The schematics of four types of impurities in WSe2. (b) Spin
lifetimes of holes of monolayer WSe2 with a relatively low hole density
$10^{11}$ cm-2 with impurities compared with experimental data[73, 74, 75,
76]. The choices of impurity concentration $n_{i}$ of different impurities are
given in the main text [16].

Experimentally there exists many types of impurities/defects in TMD samples.
Here we pick four types of impurities with different symmetries and chemical
bonds (see Fig. 3(a)) - Se vacancy (Se vac.), two neighboring Se vacancies
(2N-Se vac.), W vacancy (W vac.) and two Se vacancies with the same in-plane
position (2S-Se vac.). According to Refs. 77, 78, 79, $n_{i}$ ranges from
8$\times$1010 to 1014 cm-2 depending on samples. Considering that Se vac. is
often regarded as the most abundant impurity, we choose a reasonable $n_{i}$
of Se vac. - 7$\times$1011 cm-2, which is taken relatively low for better
comparison with experimental $\tau_{s}$ shown in Fig. 3(b). Since the
formation energy of Se vac. seems lowest and $n_{i}$ of larger impurities are
found lower than smaller impurities[79], $n_{i}$ of 2N-Se vac. is chosen as
8$\times$109 cm-2 lower than Se vac. and also for better comparison with
experimental $\tau_{s}$. $n_{i}$ of W vac. and 2S-Se vac. are chosen
arbitrarily as we found they have rather weak effects on spin relaxation and
are 7$\times$1011 and 3.5$\times$1011 cm-2, respectively.

From Fig. 3, we first find that assuming $n_{i}$ not so high, at $T$>20 K,
spin relaxation is almost driven by e-ph scattering and impurities can only
affect spin relaxation at $T$$\leq$20 K. For the effects on spin relaxation of
different impurities, we have 2N-Se vac. $\gg$ Se vac. $\gg$ W vac. $\sim$
2S-Se vac.. Moreover, the temperature dependence of $\tau_{s}$ with 2N-Se vac.
is much weaker and in better agreement with experiments than that with Se
vac.. Therefore, the observed weak temperature dependence in some experiments
is probably related the existence of larger impurities with lower symmetries.
Our observations suggest that the local symmetry and chemical bonds surround a
impurity have large impact on spin relaxation. To understand why the effects
of different impurities on spin relaxation significantly differ, further
theoretical investigations on how impurities affect impurity potentials and
their SOC corrections are required.

####  IV.2.2 Germanene under an electric field or on a substrate

In Ref. 19, through FPDM simulations, we predicted that monolayer germanene
(ML-Ge) has giant spin lifetime anisotropy, spin-valley-locking (SVL) effect
under nonzero perpendicular electric field $E_{z}$ and long $\tau_{s}$
($\sim$100 ns at 50 K), which makes it advantageous for spin-valleytronic
applications. In 2D-material-based spintronic devices, the materials are
usually supported on a substrate. Therefore, for the design of those devices,
it is crucial to understand substrate effects on spin relaxation. We further
examined $\tau_{s}$ of ML-Ge on different substrates in Ref. 21.

Figure 4: Band structures and spin textures around the Dirac cones of ML-Ge
systems with and without substrates. (a)-(d) show band structures of ML-Ge
under $E_{z}=0$ and under -7 V/nm and ML-Ge on silicane (SiH) and on InSe
substrates respectively. The red and blue bands correspond to spin-up and
spin-down states. Due to time-reversal symmetry, band structures around
another Dirac cone at $K^{\prime}=-K$ are the same except that the spin-up and
spin-down bands are reversed. The grey, white, blue, pink and green balls
correspond to Ge, H, Si, In and Se atoms, respectively. Band structures of ML-
Ge on germanane (GeH) and GaTe are shown in Supplementary Figure 4 in the
Supporting Information, and are similar to those of ML-Ge on SiH and InSe
substrates, respectively. (e) and (f) show spin textures in the
$k_{x}$-$k_{y}$ plane and 3D plots of the spin vectors ${\bf
S}_{k_{1}}^{\mathrm{exp}}$ on the circle $|\overrightarrow{k}|=0.005$ bohr-1
of the band at the band edge around $K$ of ML-Ge on SiH and InSe substrates
respectively. The color scales $S_{z}^{\mathrm{exp}}$ and the arrow length
scales the vector length of in-plane spin expectation value [19].

We show band structures and spin textures of free-standing and supported ML-Ge
in Fig. 4, which are essential for understanding spin relaxation mechanisms.
At $E_{z}=0$, ML-Ge has time-reversal and inversion symmetries, so that its
bands are Kramers degenerate[7]. A finite $E_{z}$ or a substrate breaks the
inversion symmetry and induces a strong out-of-plane ${\bf B}^{\mathrm{in}}$
(and also ${\bf S}^{\mathrm{exp}}$, Eq. 65), which splits the Kramers pairs
into spin-up and spin-down bands[19]. Interestingly, we find that band
structures of ML-Ge-SiH (Fig. 4c) and ML-Ge-GeH (not shown) are quite similar
to free-standing ML-Ge under $E_{z}$=-7 V/nm (ML-Ge@-7V/nm, Fig. 4b), which
indicates that the impact of the SiH/GeH substrate on band structure (and also
${\bf B}^{\mathrm{in}}$) may be similar to a finite $E_{z}$. This similarity
is frequently assumed in model Hamiltonian studies[80, 81]. On the other hand,
the band structures of ML-Ge-InSe (Fig. 4d) and ML-Ge-GaTe (not shown) have
more differences from the free-standing one under $E_{z}$, with larger band
gaps, smaller band curvatures at Dirac Cones, and larger electron-hole
asymmetry of band splittings. This implies that the impact of the InSe/GaTe
substrates can not be approximated by applying an $E_{z}$ to the free-standing
ML-Ge, unlike SiH/GeH substrates. We further examine ${\bf S}^{\mathrm{exp}}$
of substrate-supported ML-Ge. Importantly, from Fig. 4e and 4f, although ${\bf
S}^{\mathrm{exp}}$ of ML-Ge on substrates are highly polarized along $z$ (out-
of-plane) direction, the in-plane components of ${\bf S}^{\mathrm{exp}}$ of
ML-Ge-InSe (and ML-Ge-GaTe) are much more pronounced than ML-Ge-SiH (and ML-
Ge-GeH). Such differences are crucial to the out-of-plane spin relaxation as
discussed later.

Figure 5: The out-of-plane $\tau_{s}$ of ML-Ge under $E_{z}=0$, -7 V/nm and
substrate-supported ML-Ge as a function of temperature without impurities.
Here we show electron $\tau_{s}$ for intrinsic ML-Ge systems except that hole
$\tau_{s}$ is shown for ML-Ge-InSe, since electron $\tau_{s}$ are longer than
hole $\tau_{s}$ at low $T$ except ML-Ge-InSe [21].

We compare out-of-plane $\tau_{s}$ due to e-ph scattering between the free-
standing ML-Ge (with/without an electric field) and ML-Ge on different
substrates in Fig. 5. From Fig. 5, we find that $\tau_{s}$ of ML-Ge under
$E_{z}=0$ and -7 V/nm are at the same order of magnitude for a wide range of
temperatures. The differences are only considerable at low $T$, e.g, by 3-4
times at 20 K. On the other hand, $\tau_{s}$ of supported ML-Ge are very
sensitive to the specific substrates. While $\tau_{s}$ of ML-Ge-GeH and ML-Ge-
SiH have the same order of magnitude as the free-standing ML-Ge, in particular
very close between ML-Ge-GeH and ML-Ge@-7 V/nm, $\tau_{s}$ of ML-Ge-GaTe and
ML-Ge-InSe are shorter by at least 1-2 orders of magnitude in the whole
temperature range.

Since spins of ML-Ge at $E_{z}\neq 0$ and on a substrate are highly polarized,
spin relaxation in ML-Ge systems is dominated by EY mechanism caused by spin-
flip scattering. According to Eq. 73 and 85, $\tau_{s}^{-1}$ is roughly
proportional to density of states (DOS) and spin-mixing parameter
$\left\langle b_{z}^{2}\right\rangle$. Indeed, in Ref. 21, we find that at 300
K, the differences of $\tau_{s}$ of ML-Ge on different substrates are well
explained by the differences of the products of DOS and $\left\langle
b_{z}^{2}\right\rangle$. However, at $T\leq$50 K, the differences of the
products of DOS and $\left\langle b_{z}^{2}\right\rangle$ for different
substrates are about 3-7 times, while the differences of $\tau_{s}$ for
different substrates can be as large as two orders of magnitude. Therefore,
the substrate effects on $\tau_{s}$ can not be fully explained by the changes
of DOS and $\left\langle b_{z}^{2}\right\rangle$, in particular at relatively
low $T$.

Figure 6: The relation between $\tau_{s}^{-1}$ and
$\overline{\mathrm{sin}^{2}\left(\theta^{\uparrow\downarrow}/2\right)}$
multiplied by the scattering density of states $D^{S}$ at 20 K. $D^{S}$ is
defined in 71. $\theta^{\uparrow\downarrow}$ is the spin-flip angle between
two electronic states. For two states $\left(k,n\right)$ and
$\left(k^{\prime},n^{\prime}\right)$ with opposite spin directions,
$\theta^{\uparrow\downarrow}$ is the angle between $-{\bf
S}_{kn}^{\mathrm{exp}}$ and ${\bf S}_{k^{\prime}n^{\prime}}^{\mathrm{exp}}$.
$\overline{\mathrm{sin}^{2}\left(\theta^{\uparrow\downarrow}/2\right)}$ is
defined in Eq. 89. The variation of $D^{S}$ among different substrates is at
most three times, much weaker than the variations of $\tau_{s}^{-1}$ and other
quantities shown here [21].

Since
$\tau_{s}^{-1}\propto\overline{|\widetilde{g}^{\uparrow\downarrow}|^{2}}D^{\mathrm{S}}$
(Eq. 73), to understand substrate effects on $\tau_{s}$ at low $T$, we first
compare $\tau_{s}^{-1}$ and
$\overline{|\widetilde{g}^{\uparrow\downarrow}|^{2}}D^{\mathrm{S}}$ of
different ML-Ge systems. In Ref. 21, we found that $\tau_{s}^{-1}$ is almost
linearly proportional to
$\overline{|\widetilde{g}^{\uparrow\downarrow}|^{2}}D^{S}$ at 20 K. As the
variation of $D^{S}$ among ML-Ge on different substrates is at most three
times, which is much weaker than the large variation of $\tau_{s}^{-1}$, this
indicates that the substrate-induced change of $\tau_{s}$ is mostly due to the
substrate-induced change of spin-flip matrix elements.

To have deeper intuitive understanding, we then propose an important
electronic quantity for intervalley spin-flip scattering - the spin-flip angle
$\theta^{\uparrow\downarrow}$ between two electronic states. For two states
$\left(k_{1},n_{1}\right)$ and $\left(k_{2},n_{2}\right)$ with opposite spin
directions, $\theta^{\uparrow\downarrow}$ is the angle between $-{\bf
S}_{k_{1}n_{1}}^{\mathrm{exp}}$ and ${\bf S}_{k_{2}n_{2}}^{\mathrm{exp}}$ or
equivalently the angle between $-{\bf B}_{k_{1}}^{\mathrm{in}}$ and ${\bf
B}_{k_{2}}^{\mathrm{in}}$. At low $T$, due to large SOC splittings of
conduction and valence bands of supported ML-Ge, intravalley spin-flip e-ph
scattering processes are forbidden, because the corresponding phonons have too
high energies and are not occupied. So spin relaxation in supported ML-Ge is
dominated by intervalley spin-flip e-ph scattering. Therefore,
$\theta^{\uparrow\downarrow}$ is helpful for understanding spin relaxation in
supported ML-Ge.

The motivation of examining $\theta^{\uparrow\downarrow}$ is that: Suppose two
wavevectors ${\bf k}_{1}$ and ${\bf k}_{2}={\bf-k}_{1}$ are in two opposite
valleys $Q$ and -$Q$ respectively and there is a pair of bands, which are
originally Kramers degenerate but splitted by ${\bf B}^{\mathrm{in}}$. Due to
time-reversal symmetry, we have ${\bf B}_{k_{1}}^{\mathrm{in}}=-{\bf
B}_{k_{2}}^{\mathrm{in}}$, which means the two states at the same band $n$ at
${\bf k_{1}}$ and ${\bf k_{2}}$ have opposite spins and
$\theta^{\uparrow\downarrow}$ between them is zero. Therefore, the matrix
element of operator $\widehat{A}$ between states $\left(k_{1},n\right)$ and
$\left(k_{2},n\right)$ \- $A_{k_{1}n,k_{2}n}$ is a spin-flip one and we name
it as $A_{k_{1}k_{2}}^{\uparrow\downarrow}$. According to Ref. 67, with time-
reversal symmetry, $A_{k_{1}k_{2}}^{\uparrow\downarrow}$ is exactly zero. In
general, for another wavevector ${\bf k_{3}}$ within valley -$Q$ but not
${\bf-k}_{1}$, $A_{k_{1}k_{3}}^{\uparrow\downarrow}$ is usually non-zero. One
critical quantity that determines the intervalley spin-flip matrix element
$A_{k_{1}k_{3}}^{\uparrow\downarrow}$ for a band within the pair introduced
above is $\theta_{k_{1}k_{3}}^{\uparrow\downarrow}$. Based on time-independent
perturbation theory, we can prove that $\left|A^{\uparrow\downarrow}\right|$
between two states is approximately proportional to
$\left|\sin\left(\theta^{\uparrow\downarrow}/2\right)\right|$. The derivation
is given in Sec. III.2.3.

As shown in Figure 6c, $\tau_{s}^{-1}$ of ML-Ge on different substrates at 20
K is almost linearly proportional to
$\overline{\sin^{2}\left(\theta^{\uparrow\downarrow}/2\right)}D^{S}$, where
$\overline{\sin^{2}\left(\theta^{\uparrow\downarrow}/2\right)}$ is the
statistically-averaged modulus square of
$\sin\left(\theta^{\uparrow\downarrow}/2\right)$. This indicates that the
relation
$\overline{|\widetilde{g}^{\uparrow\downarrow}|^{2}}\propto\overline{\sin^{2}\left(\theta^{\uparrow\downarrow}/2\right)}$
is nearly perfectly satisfied at low $T$, where intervalley processes dominate
spin relaxation. We additionally examined the relation between $\tau_{s}^{-1}$
and $\overline{\sin^{2}\left(\theta^{\uparrow\downarrow}/2\right)}D^{S}$ at
300 K and found that the trend of $\tau_{s}^{-1}$ is still approximately
captured by the trend of
$\overline{\sin^{2}\left(\theta^{\uparrow\downarrow}/2\right)}D^{S}$.

Since $\theta^{\uparrow\downarrow}$ is defined by ${\bf S}^{\mathrm{exp}}$ at
different states, $\tau_{s}$ is highly correlated with ${\bf
S}^{\mathrm{exp}}$ and more specifically with the anisotropy of ${\bf
S}^{\mathrm{exp}}$ (equivalent to the anisotropy of ${\bf B}^{\mathrm{in}}$).
Qualitatively, the larger anisotropy of ${\bf S}^{\mathrm{exp}}$ leads to
smaller $\theta^{\uparrow\downarrow}$ (consistent with Fig. 4e and f) and
longer $\tau_{s}$ along the high-spin-polarization direction. This finding may
be applicable to spin relaxation in other materials whenever intervalley spin-
flip scattering dominates or spin-valley locking exists, e.g., in TMDs[72],
Stanene[82], 2D hybrid perovskites with persistent spin helix[6], etc.

###  IV.3 DP (Dyakonov-Perel) systems

In many non-magnetic materials with broken inversion symmetry, spin relaxation
is dominated by DP mechanism.

####  IV.3.1 GaAs

Spin dynamics in GaAs has broad interest in spintronics over past decades[7,
47, 83, 84, 85] and more recently[86, 87, 88], partly due to its long spin
lifetime in n-type GaAs at low temperatures[47]. Despite various experimental
[47, 83, 89, 48, 90] and theoretical [7, 91, 84, 92, 85, 93] (mostly using
parameterized model Hamiltonian) studies previously, the dominant spin
relaxation mechanism in bulk GaAs at various $T$ and doping concentrations
$n_{i}$ remains unclear.

Through FPDM simulations of $n$-GaAs at various $T$ and $n_{i}$ with different
scattering mechanisms in Ref. 16, we pointed out that although at low
temperatures and moderate doping concentrations e-i scattering dominates
carrier relaxation, e-e scattering is the most dominant process in spin
relaxation. Moreover, we found that the relative contributions of phonon modes
vary considerably between spin and carrier relaxation.

Figure 7: (a) $\tau_{s}$ and (b) $\tau_{p}$ of $n$-GaAs with different doping
concentrations $n_{i}$ at 30 K with different scattering mechanisms. “All”
represents all the e-ph, e-i and e-e scattering mechanisms being considered.
(c) $\Delta\Omega_{\perp z}^{2}$ as a function of carrier density $n$, where
${\bf\Delta\Omega_{\perp z}}$ is the fluctuation amplitude of Larmor frequency
defined in Eq. 91 [16].

Figure 7 shows $\tau_{s}$ and $\tau_{p}$ with different $n_{i}$ at 30 K with
individual and total scattering pathways, respectively. It is found that the
roles of different scattering mechanism differ considerably between spin and
carrier relaxation processes. Specifically, for the carrier relaxation in Fig.
7b, except when $n_{i}$ is very low (e.g. at $10^{14}$ cm-3), the e-i
scattering dominates. On the other hand, for the spin relaxation in Fig. 7a,
the e-e scattering dominates except at very high concentration (above
$10^{17}$ cm-3), while e-i scattering is only important in the very high
doping region (close to or above $10^{17}$ cm-3).

Figure 7 shows the calculated $\tau_{s}$ has a maximum at
$n_{i}\,=\,1\text{-}2\times 10^{16}$ cm-3, and $\tau_{s}$ decreases fast with
$n_{i}$ going away from its peak position. This is in good agreement with the
experimental finding in Ref. 47, which also reported $\tau_{s}$ at
$n_{i}=10^{16}$ cm-3 is longer than $\tau_{s}$ at other lower and higher
$n_{i}$ at a low temperature (a few Kelvin). The $n_{i}$ dependence of
$\tau_{s}$ may be qualitatively interpreted from the DP relation[7]
$\tau_{s,z}\sim\tau_{p}\Delta\Omega_{\perp z}^{2}$ (Eq. 92), where
${\bf\Delta\Omega_{\perp z}}$ is the fluctuation amplitude of Larmor frequency
defined in Eq. 91. From Fig. 7, we find that with $n_{i}$ from $10^{14}$ cm-3
to $5\times 10^{15}$ cm-3, $\tau_{p}$ decreases rapidly (black curve in Fig.
7b) and $\Delta\Omega_{\perp z}^{2}$ remains flat in Fig. 7c, which may
explain why $\tau_{s}$ increases in Fig. 7a based on the DP relation; however,
when $n_{i}>10^{16}$ cm-3, $\tau_{p}$ decreases with a similar speed but
$\Delta\Omega_{\perp z}^{2}$ experiences a sharp increase, which may explain
why $\tau_{s}$ decreases in Fig. 7b and owns a maximum at $10^{16}$ cm-3.

Note that although the DP relation is intuitive to understand the cause of
doping-level dependence of spin lifetime, it may break down when we evaluate
individual scattering processes. For example, when $n_{i}$ increases from
$10^{14}$ cm-3 to $10^{15}$ cm-3, both carrier lifetime $\tau_{s}$ and
$\tau_{s,z}$ due to e-i scattering decrease while the internal magnetic field
remains unchanged. Moreover, the simple empirical relation cannot possibly
explain our first-principles results that the e-e and e-i scatterings have
largely different contributions in carrier and spin relaxation. First-
principles calculations are critical to provide unbiased mechanistic insights
to spin and carrier relaxation of general systems.

####  IV.3.2 Graphene on hBN

Graphene samples exhibit exciting spintronic properties such as long
$\tau_{s}$ and $l_{s}$ at room temperature [1, 94]. In practice, graphene is
usually supported by a substrate[95, 96, 97, 94, 98]. Actually $\tau_{s}$ of
the free-standing graphene sample was found relatively low[99] $\sim$150 ps at
300 K (compared with the longest reported value $\sim$12 ns), probably due to
free-standing graphene samples often have more imperfections. Therefore, it is
important to understand spin relaxation in supported graphene.

Figure 8: Theoretical spin and carrier relaxation in graphene (red lines),
graphene at $E_{\perp}$=0.4 V/nm (Gr+$E_{\perp}$, green lines) and graphene on
hBN (Gr-hBN, blue lines) as a function of $n$ (positive or negative for
electron and hole doping respectively) at room temperature. (a) In-plane spin
lifetime $\tau_{s\parallel}$, (b) out-of-plane spin lifetime $\tau_{s\perp}$
and (c) their ratios $\tau_{s\perp}/\tau_{s\parallel}$ calculated using FPDM
approach (solid lines) compared with those estimated from EY/DP models with
first-principles inputs (dashed for conventional DP and dotted for modified DP
with intervalley scattering contribution lines). (d) Total $\tau_{p}$ (three
dashed lines) and intervalley only contribution for Gr-hBN (dotted blue line),
(e) $\left\langle b_{i}^{2}\right\rangle$ and (f)
$\left\langle\Omega_{i}^{2}\right\rangle$, all predicted from first-principles
for the DP and EY models estimates of $\tau_{s}$ in (a-c) [20].

From Fig. 8, we find that the EY model for graphene, and the DP model for
Gr+$E_{\perp}$ and Gr-hBN, agree qualitatively with FPDM predictions, but with
some important quantitative differences discussed next. First, the EY model
for graphene is more accurate for electrons than for holes, for both
$\tau_{s\parallel}$ and $\tau_{s\perp}$ (Fig. 8(a-b)). The conventional DP
model matches FPDM predictions quantitatively for both $\tau_{s\parallel}$ and
$\tau_{s\perp}$ of Gr+$E_{\perp}$, but only for $\tau_{s\perp}$ of Gr-hBN.

The discrepancy of the conventional DP model for $\tau_{s\parallel}$ of Gr-hBN
can be rectified by modifying the model. Briefly, the DP model assumes that
${\bf B}^{\mathrm{in}}$ effectively changes randomly each time the electron
scatters. The in-plane magnetic field ${\bf B}_{||}^{\mathrm{in}}$ rotates
over the Fermi circle and covers all in-plane directions, satisfying this
condition, in both Gr+$E_{\perp}$ and Gr-hBN. However, the out-of-plane
magnetic field ${\bf B}_{\perp}^{\mathrm{in}}$, which matters only for
$\tau_{s\parallel}$ and is present only for Gr-hBN, has the same direction
within each valley. Consequently, only intervalley scattering will change the
${\bf B}_{\perp}^{\mathrm{in}}$ for a given electron spin. As proposed in Ref.
100, this can be captured by changing the DP model from
$(\tau_{s,x}^{\mathrm{DP}})^{-1}\sim\tau_{p}\left(\Delta{\bf\Omega}_{\perp
x}\right)^{2}=\tau_{p}\left\langle\Omega_{y}^{2}+\Omega_{z}^{2}\right\rangle$
as given by Eq. 92 for in-plane $x$ spins, to

$\left(\tau_{s,x}^{\mathrm{mDP}}\right)^{-1}\approx\tau_{p}\left\langle\Omega_{y}^{2}\right\rangle+\tau_{p}^{\mathrm{Inter}}\left\langle\Omega_{z}^{2}\right\rangle,$
(99)

where $\tau_{p}^{\mathrm{Inter}}$ is the intervalley scattering time (dotted
line in Fig. 8(d)). This modified DP model agrees with FPDM predictions for
$\tau_{s\parallel}$ of Gr-hBN (Fig. 8(a)).

The ratio $\tau_{s\perp}/\tau_{s\parallel}$ (Fig. 8(c)) is nearly 1/2 for
graphene, consistent with the EY relation Eq. 85 and the fact that
$\left\langle b_{||}^{2}\right\rangle/\left\langle b_{\perp}^{2}\right\rangle$
is also 1/2 (Fig. 8(e)). This ratio remains unchanged for Gr+$E_{\perp}$, but
now because $\left(\Delta{\bf\Omega}_{\perp
z}\right)^{2}=\langle\Omega_{x}^{2}+\Omega_{y}^{2}\rangle=2\langle\Omega_{y}^{2}\rangle$,
while $\left(\Delta{\bf\Omega}_{\perp
x}\right)^{2}=\langle\Omega_{y}^{2}\rangle$ since $\Omega_{z}=0$ (Fig. 8(f)),
leading to $\tau_{s,x}^{\mathrm{DP}}=2\tau_{s,z}^{\mathrm{DP}}$ using Eq. 92.
This ratio deviates substantially from 1/2 only for Gr-hBN (Fig. 8(c)) due to
the substrate-induced $\Omega_{z}\neq 0$. The conventional DP model (Eq. 92)
only captures part of this dramatic effect seen in the FPDM calculations,
while the modifications in Eq. 99 account for $\Omega_{z}$ correctly and agree
with the FPDM results in Fig. 8(c).

$\tau_{s}$ decrease with increasing carrier density magnitude in graphene
(Fig. 8(a,b)), but this trend reverses for both inversion-symmetry-broken
cases, in agreement with some experiments.[94, 101] The overall $\tau_{s}$ are
reduced by one-two orders of magnitude in free-standing graphene by
$E_{\perp}$ of 0.4 V/nm, down from $\mu s$ to tens of ns.
$\left\langle\Omega_{i}^{2}\right\rangle$ of Gr-hBN in Fig. 8(f) is about 100
times larger than that of Gr+$E_{\perp}$, further reducing $\tau_{s}$ of Gr-
hBN to the ns scale, comparable to experimental measurements.[97, 102, 96, 94]

Finally, $\tau_{s}$ is mostly symmetric between electrons and holes for
Gr+$E_{\perp}$ (Fig. 8(a, b)). However, we find hole $\tau_{s}$ to be
typically 2-3 times smaller than electrons for both free-standing graphene and
Gr-hBN. On hBN, this asymmetry is captured by the (modified) DP model and
stems primarily from the larger spin-splitting and hence
$\left\langle\Omega_{i}^{2}\right\rangle$ in the valence band compared to the
conduction band (not shown), consistent with previous calculations.[103]
Importantly, this effect depends sensitively on the substrate, and even on
hBN, could reverse for a different layer stacking.[103] Consequently,
experiments may find electron-hole asymmetries of either sign depending on the
substrate and precise structure,[96, 101, 97, 102, 94, 98, 103] and we focus
here on the comparison between FPDM predictions and DP model for the specific
lowest-energy stacking of Gr-hBN.

###  IV.4 Spin dephasing under magnetic field

As introduced in Sec. II.2.1, the dephasing time (or transverse time) of a
spin ensemble is called $T_{2}^{*}$ and describes the decay of the total
excess spin ${\bf S}-{\bf S}^{\mathrm{eq}}$ at an nonzero ${\bf
B}^{\mathrm{ext}}$ perpendicular to ${\bf S}-{\bf S}^{\mathrm{eq}}$. As
discussed in Sec. II.1.2, $T_{2}^{*}$ can be simulated straightforwardly
through FPDM simulations with Zeeman Hamiltonian considering both spin and
orbital angular momenta (Eq. 13 and 14). For the purpose of analysing and
understanding spin dephasing, one key parameter is the Landé $g$-factor
$\widetilde{g}$ (Eq. 96). It value relates to ${\bf B}^{\mathrm{ext}}$-induced
energy splitting (Zeeman effect) $\Delta E\left({\bf B}^{\mathrm{ext}}\right)$
and Larmor precession frequency $\Omega$, satisfying $\Omega\approx\Delta
E=\mu_{B}B^{\mathrm{ext}}\widetilde{g}$. More importantly, the $g$-factor
fluctuation (near Fermi surface or $\mu_{F,c}$) $\Delta\widetilde{g}$ (Eq. 97)
leads to a nonzero $\Delta\Omega$ (Eq. 98) and then $T_{2}^{*}$ due to DP or
FID mechanism.

Spintronics in halide perovskites has drawn significant attention in recent
years, due to highly tunable spin-orbit fields and intriguing interplay with
lattice symmetry. Recently, we simulate[18] $\widetilde{g}$,
$\Delta\widetilde{g}$ and $T_{2}^{*}$ of a typical halide perovskite -
CsPbBr3. We find that $T_{2}^{*}$ is sensitive to $B^{\mathrm{ext}}$ at $T<$20
K but not at $T\geq$20 K. At 4 K, we predict that
$\left(T_{2}^{*}\right)^{-1}$ is linear to $B^{\mathrm{ext}}$ at
$B^{\mathrm{ext}}\geq$0.4 Tesla and the predicted slope of
$\left(T_{2}^{*}\right)^{-1}$ is in good agreement with experimental data.
Moreover, we predict strong $n$-dependence of $\widetilde{g}$,
$\Delta\widetilde{g}$ and $T_{2}^{*}$. Together with FPDM simulations of spin
relaxation time $T_{1}$ of CsPbBr3 at various conditions, our work provides
fundamental insights on how to control and manipulate spin
relaxation/dephasing in halide perovskites, which are vital for their
spintronics and quantum information applications.

##  V Outlooks

The first-principles density-matrix dynamics (FPDM) approach is an important
technique for studying spin and electron dynamics and transport, with
outstanding advantages: i) it can accurately describe various interactions,
scattering processes, and spin precession simultaneously; ii) it can describe
processes far from equilibrium; iii) it can simulate dynamical processes on
various time scales from femtoseconds to milliseconds. We have shown its
success in simulating ultrafast spin dynamics, spin and charge transport,
relaxation and dephasing. In the future study, the FPDM approach still has a
broad open area for new theory development.

####  V.0.1 Exciton dynamics

Exciton is a bound electron-hole pair generated by optical excitation.
Excitons play an important role in optical properties of semiconductors, in
particular in low-dimensional systems. Exciton dynamics has been extensively
studied using different methods, in particular non-equilibrium Greens’
function theory (NEGF) within Kadanoff-Baym equations[104, 105, 45]. We focus
on theoretical simulations of exciton spin relaxation based on Lindbladian DM
master equation with quantum description of the exciton-phonon scattering.
Understanding the exciton spin relaxation is also important for understanding
optical measurements of spin dynamics. By replacing electrons to excitons, the
existing FPDM approach can be generalized to simulation spin dynamics of
excitons. The exciton dynamics can involve many types of processes including
exciton-light interaction, exciton-phonon scattering, phonon-mediated exciton
recombination and dissociation, exciton–exciton annihilation, etc. Initially,
we will focus on exciton-light interaction and exciton-phonon scattering,
which are two most important processes leading to exciton spin relaxation.

The exciton-light interaction is responsible to exciton spin generation by
light absorption with a circularly-polarized pump pulse and exciton spin
relaxation through exciton-phonon scattering and exciton radiative
recombination. The light absorption for excitons is similar to that for
electrons (see Sec. II.1.2) except that the momentum matrix elements are now
between two excitonic states [106, 107, 108].

The exciton-phonon scattering is responsible to exciton spin relaxation.
According to Refs. 109, 110, the exciton-phonon scattering matrix elements can
be approximately obtained from the excition wavefunction at finite momentum
and electron-phonon scattering matrix elements. Recently we implemented the
exciton-phonon scattering matrix elements and used them to simulate phonon-
assisted indirect exciton radiative recombination. In the density-matrix
master equation, the term of the exciton-phonon scattering is rather similar
to the electron-phonon scattering, except that the excitons are bosons while
the electrons are fermions. Therefore, the implementation of the exciton-
phonon scattering in the frame of FPDM approach can be straightforward.

####  V.0.2 Circular photogalvanic effect (CPGE)

CPGE is the effect that under the circularly-polarized light a DC current may
be induced in a solid-state material, and widely presents in materials without
inversion symmetry, in absence of p-n junction and applied electric field. In
recent years, CPGE has attracted growing interests in the fields of
topological physics and spin-optotronics [111, 112]. In order to deeply
understand CPGE, theoretical methods have been developed based on perturbation
theory. However, several critical issues remain: (i) the scattering mechanisms
are highly simplified with a single relaxation time approximation; (ii)
spontaneous recombination was rarely considered; (iii) available theory
requires different formulation for each contribution i.e. interband or
intraband contributions; (iv) many-body effects were not considered. These
issues may be resolved based on our FPDM approach.

If a constant laser field is applied to an inversion symmetry broken system,
the quantum master equation for describing density matrix of the electronic
system reads:

$\displaystyle\frac{d\rho}{dt}=$
$\displaystyle\frac{d\rho}{dt}|_{\mathrm{laser}}+\frac{d\rho}{dt}|_{\mathrm{sp-
em}}+\frac{d\rho}{dt}|_{\mathrm{e-ph}},$ (100)

where $\frac{d\rho}{dt}|_{\mathrm{laser}}$ is coherent dynamics due to a laser
field (Eq. 12), $\frac{d\rho}{dt}|_{\mathrm{sp-em}}$ and
$\frac{d\rho}{dt}|_{\mathrm{e-ph}}$ are the spontaneous emission (Eq. 40) and
e-ph scattering terms (Eq. 19 and 20), respectively.

CPGE is measuring the steady state current when a nonmagnetic system is under
a constant circularly-polarized laser field for a long enough time. In the
velocity gauge, the current density is computed using[32]

$\displaystyle{\bf J}=$ $\displaystyle\mathrm{Tr}\left(\rho{\bf j}\right),$
(101) $\displaystyle{\bf j}_{kmn}\left(t\right)=$ $\displaystyle-e\left({\bf
v}_{kmn}+\delta_{mn}\frac{e}{m_{e}}{\bf A}\left(t\right)\right),$ (102)

where ${\bf v}$ is the velocity operator matrix and ${\bf v}={\bf p}/m_{e}$.
The average value of the second term of ${\bf j}$ is usually zero so that this
term does not contribute to CPGE.

####  V.0.3 More on scattering terms

Although in our FPDM approach, the _ab-initio_ treatment of quantum scattering
is rather general, the following theoretical development can further improve
the applicability of our methods to different systems:

##### Anharmonic phonons.

In current implementation, phonons are harmonic and it is required that the
material must be dynamically stable at zero temperature. However, for some
soft materials (e.g., hybrid halide perovskites), significant anharmonic
effects appear at high-$T$ phase, which requires extracting phonon properties
from finite-temperature simulations, e.g. using _ab-initio_ molecular
dynamics. Such approach has been implemented [113] for studying phonon,
electron-phonon, and carrier transport properties, which can be also applied
to study the anharmonic effect on spin relaxation under our FPDM framework.

##### Fröhlich interaction and LO-TO splitting in doped semiconductor.

The intraband dielectric screening is not considered for Fröhlich interaction
and LO-TO splitting in current implementation. This may be problematic in
doped semiconductor and may lead to significant errors at moderate or high
carrier density. Since we have already implemented the intraband dielectric
function (Eq. 37), its effect can be straightforwardly included in FPDM
approach similar to Ref. 114.

##### Short-range contribution to the electron-ionized-impurity scattering.

For the electron-ionized-impurity scattering, we currently only consider the
long-range contribution. This is appropriate for certain systems like GaAs,
but may be problematic when intravalley processes dominate spin relaxation or
the short-range spin-flip electron-ionized-impurity scattering is unimportant
to spin relaxation. Therefore, it is important to include both long-range and
short-range contributions to the electron-ionized-impurity scattering. Similar
to Ref. 56, the long-range and short-range parts are treated separately - the
short-range part is treated from first principles similar to neutral impurity,
and the long-range part is simulated using the screened Coulomb potential.

## Acknowledgements

This work is supported by National Science Foundation under grant No.
DMR-1956015. This research used resources of the Center for Functional
Nanomaterials, which is a US DOE Office of Science Facility, and the
Scientific Data and Computing center, a component of the Computational Science
Initiative, at Brookhaven National Laboratory under Contract No. DE-SC0012704,
the lux supercomputer at UC Santa Cruz, funded by NSF MRI grant AST 1828315,
the National Energy Research Scientific Computing Center (NERSC) a U.S.
Department of Energy Office of Science User Facility operated under Contract
No. DE-AC02-05CH11231, the Extreme Science and Engineering Discovery
Environment (XSEDE) which is supported by National Science Foundation Grant
No. ACI-1548562 [115].

## References

  * Avsar  _et al._ [2020] A. Avsar, H. Ochoa, F. Guinea, B. Özyilmaz, B. J. van Wees, and I. J. Vera-Marun, Rev. Mod. Phys. 92, 021003 (2020).
  * Sierra  _et al._ [2021] J. F. Sierra, J. Fabian, R. K. Kawakami, S. Roche, and S. O. Valenzuela, Nat. Nanotechnol 16, 856 (2021).
  * Garcia  _et al._ [2020] J. H. Garcia, M. Vila, C.-H. Hsu, X. Waintal, V. M. Pereira, and S. Roche, Phys. Rev. Lett. 125, 256603 (2020).
  * Šmejkal  _et al._ [2018] L. Šmejkal, Y. Mokrousov, B. Yan, and A. H. MacDonald, Nat. Phys. 14, 242 (2018).
  * Kim  _et al._ [2021] Y.-H. Kim, Y. Zhai, H. Lu, X. Pan, C. Xiao, E. A. Gaulding, S. P. Harvey, J. J. Berry, Z. V. Vardeny, J. M. Luther, _et al._ , Science 371, 1129 (2021).
  * Zhang  _et al._ [2022] L. Zhang, J. Jiang, C. Multunas, C. Ming, Z. Chen, Y. Hu, Z. Lu, S. Pendse, R. Jia, M. Chandra, _et al._ , Nat. Photon. , 1 (2022).
  * Žutić  _et al._ [2004] I. Žutić, J. Fabian, and S. D. Sarma, Rev. Mod. Phys. 76, 323 (2004).
  * Wu  _et al._ [2010] M. Wu, J. Jiang, and M. Weng, Phys. Rep. 493, 61 (2010).
  * Restrepo and Windl [2012] O. D. Restrepo and W. Windl, Phys. Rev. Lett. 109, 166604 (2012).
  * Fedorov  _et al._ [2013] D. V. Fedorov, M. Gradhand, S. Ostanin, I. V. Maznichenko, A. Ernst, J. Fabian, and I. Mertig, Phys. Rev. Lett. 110, 156602 (2013).
  * Marques  _et al._ [2012] M. A. L. Marques, N. T. Maitra, F. M. S. Nogueira, E. K. U. Gross, and A. Rubio,  _Fundamentals of Time-Dependent Density Functional Theory_ , Vol. 837 (Springer Science & Business Media, 2012).
  * Chen and Wang [2019] Z. Chen and L. Wang, Sci. Adv. 5, eaau8000 (2019).
  * Acharya  _et al._ [2020] S. R. Acharya, V. Turkowski, G. P. Zhang, and T. S. Rahman, Phys. Rev. Lett. 125, 017202 (2020).
  * Krieger  _et al._ [2015] K. Krieger, J. K. Dewhurst, P. Elliott, S. Sharma, and E. K. U. Gross, J. Chem. Theory Comput. 11, 4870 (2015).
  * Park  _et al._ [2022] J. Park, J.-J. Zhou, Y. Luo, and M. Bernardi, Phys. Rev. Lett. 129, 197201 (2022).
  * Xu  _et al._ [2021a] J. Xu, A. Habib, R. Sundararaman, and Y. Ping, Phys. Rev. B 104, 184418 (2021a).
  * Xu  _et al._ [2020] J. Xu, A. Habib, S. Kumar, F. Wu, R. Sundararaman, and Y. Ping, Nat. Commun. 11, 2780 (2020).
  * Xu  _et al._ [2023] J. Xu, K. Li, U. N. Huynh, J. Huang, R. Sundararaman, V. Vardeny, and Y. Ping, Preprint at https://arxiv.org/abs/2210.17074 (2023).
  * Xu  _et al._ [2021b] J. Xu, H. Takenaka, A. Habib, R. Sundararaman, and Y. Ping, Nano Lett. 21, 9594 (2021b).
  * Habib  _et al._ [2022] A. Habib, J. Xu, Y. Ping, and R. Sundararaman, Phys. Rev. B 105, 115122 (2022).
  * Xu and Ping [2023] J. Xu and Y. Ping, npj Comput. Mater. 9, 47 (2023).
  * Rossi and Kuhn [2002] F. Rossi and T. Kuhn, Rev. Mod. Phys. 74, 895 (2002).
  * Bonitz [2016] M. Bonitz,  _Quantum Kinetic Theory_ (Springer, Cham, 2016).
  * Axt and Stahl [1994] V. M. Axt and A. Stahl, Z. Physik B - Condensed Matter 93, 195 (1994).
  * Iotti and Rossi [2017] R. C. Iotti and F. Rossi, Eur. Phys. J. B 90, 1 (2017).
  * Marini [2013] A. Marini, in  _Journal of Physics: Conference Series_ , Vol. 427 (IOP Publishing, 2013) p. 012003.
  * Kohn and Luttinger [1957] W. Kohn and J. Luttinger, Phys. Rev. 108, 590 (1957).
  * Krieger and Iafrate [1987] J. Krieger and G. Iafrate, Phys. Rev. B 35, 9644 (1987).
  * Ciancio  _et al._ [2004] E. Ciancio, R. C. Iotti, and F. Rossi, Phys. Rev. B 69, 165319 (2004).
  * Kané  _et al._ [2012] G. Kané, M. Lazzeri, and F. Mauri, Phys. Rev. B 86, 155433 (2012).
  * Sekine  _et al._ [2017] A. Sekine, D. Culcer, and A. H. MacDonald, Phys. Rev. B 96, 235134 (2017).
  * Ventura  _et al._ [2017] G. Ventura, D. Passos, J. L. Dos Santos, J. V. P. Lopes, and N. Peres, Phys. Rev. B 96, 035431 (2017).
  * Iafrate  _et al._ [2017] G. Iafrate, V. Sokolov, and J. Krieger, Phys. Rev. B 96, 144303 (2017).
  * King-Smith and Vanderbilt [1993] R. D. King-Smith and D. Vanderbilt, Phys. Rev. B 47, 1651 (1993).
  * Resta [1994] R. Resta, Rev. Mod. Phys. 66, 899 (1994).
  * Multunas  _et al._ [2022] C. Multunas, A. Grieder, J. Xu, Y. Ping, and S. R., Under Preparation (2022).
  * Rosati  _et al._ [2014] R. Rosati, R. C. Iotti, F. Dolcini, and F. Rossi, Phys, Rev. B 90, 125140 (2014).
  * Giustino [2017] F. Giustino, Rev. Mod. Phys. 89, 015003 (2017).
  * Sundararaman and Ping [2017] R. Sundararaman and Y. Ping, J. Chem. Phys. 146, 104109 (2017).
  * Jacoboni [2010] C. Jacoboni,  _Theory of Electron Transport in Semiconductors: A Pathway from Elementary Physics to Nonequilibrium Green Functions_ , Vol. 165 (Springer Science & Business Media, 2010).
  * Wu  _et al._ [2005] X. Wu, D. Vanderbilt, and D. R. Hamann, Phys, Rev. B 72, 035105 (2005).
  * Lü  _et al._ [2007] C. Lü, J. Cheng, M. Wu, and I. da Cunha Lima, Phys. Lett. A 365, 501 (2007).
  * Meier and Zakharchenya [2012] F. Meier and B. P. Zakharchenya,  _Optical Orientation_ (Elsevier, 2012).
  * Dean  _et al._ [2016] M. P. M. Dean, Y. Cao, X. Liu, S. Wall, D. Zhu, R. Mankowsky, V. Thampy, X. M. Chen, J. G. Vale, D. Casa, _et al._ , Nat. Mater. 15, 601 (2016).
  * Molina-Sánchez  _et al._ [2017] A. Molina-Sánchez, D. Sangalli, L. Wirtz, and A. Marini, Nano Lett. 17, 4549 (2017).
  * Mainkar  _et al._ [1996] N. Mainkar, D. A. Browne, and J. Callaway, Phys, Rev. B 53, 3692 (1996).
  * Kikkawa and Awschalom [1998] J. M. Kikkawa and D. D. Awschalom, Phys. Rev. Lett. 80, 4313 (1998).
  * Kimel  _et al._ [2001] A. V. Kimel, F. Bentivegna, V. N. Gridnev, V. V. Pavlov, R. V. Pisarev, and T. h. Rasing, Phys, Rev. B 63, 235201 (2001).
  * Sundararaman  _et al._ [2017] R. Sundararaman, K. Letchworth-Weaver, K. A. Schwarz, D. Gunceler, Y. Ozhabes, and T. A. Arias, SoftwareX 6, 278 (2017).
  * Marzari and Vanderbilt [1997] N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 (1997).
  * Brown  _et al._ [2016a] A. M. Brown, R. Sundararaman, P. Narang, W. A. Goddard, and H. A. Atwater, ACS Nano 10, 957 (2016a).
  * Narang  _et al._ [2017] P. Narang, L. Zhao, S. Claybrook, and R. Sundararaman, Adv. Opt. Mater. 5, 1600914 (2017).
  * Brown  _et al._ [2017] A. M. Brown, R. Sundararaman, P. Narang, A. M. Schwartzberg, W. A. Goddard III, and H. A. Atwater, Phys. Rev. Lett. 118, 087401 (2017).
  * Habib  _et al._ [2018] A. Habib, R. Florio, and R. Sundararaman, J. Opt. 20, 064001 (2018).
  * Brown  _et al._ [2016b] A. M. Brown, R. Sundararaman, P. Narang, W. A. Goddard III, and H. A. Atwater, Phys. Rev. B 94, 075120 (2016b).
  * Verdi and Giustino [2015] C. Verdi and F. Giustino, Phys. Rev. Lett. 115, 176401 (2015).
  * Sohier  _et al._ [2017] T. Sohier, M. Gibertini, M. Calandra, F. Mauri, and N. Marzari, Nano Lett. 17, 3758 (2017).
  * Iotti and Rossi [2001] R. C. Iotti and F. Rossi, Phys. Rev. Lett. 87, 146603 (2001).
  * Watanabe and Yanase [2021a] H. Watanabe and Y. Yanase, Phys. Rev. X 11, 011001 (2021a).
  * Kubo [1966] R. Kubo, Rep. Prog. Phys. 29, 255 (1966).
  * Maassen  _et al._ [2011] T. Maassen, F. Dejene, M. Guimarães, C. Józsa, and B. Van Wees, Physical Review B 83, 115410 (2011).
  * Ciccarino  _et al._ [2018] C. J. Ciccarino, T. Christensen, R. Sundararaman, and P. Narang, Nano Lett. 18, 5709 (2018).
  * Gunst  _et al._ [2016] T. Gunst, T. Markussen, K. Stokbro, and M. Brandbyge, Phys. Rev. B 93, 035414 (2016).
  * Mahan [2000] G. D. Mahan,  _Many-Particle Physics_ (Springer US, 2000).
  * Fabian and Sarma [1998] J. Fabian and S. D. Sarma, Phys. Rev. Lett. 81, 5624 (1998).
  * Leyland  _et al._ [2007] W. Leyland, R. Harley, M. Henini, A. Shields, I. Farrer, and D. Ritchie, Phys. Rev. B 76, 195305 (2007).
  * Yafet [1963] Y. Yafet, in  _Solid state physics_ , Vol. 14 (Elsevier, 1963) pp. 1–98.
  * Yu [2016] Z. Yu, Sci. Rep. 6, 1 (2016).
  * Kirstein  _et al._ [2022] E. Kirstein, D. Yakovlev, M. Glazov, E. Zhukov, D. Kudlacik, I. Kalitukha, V. Sapega, G. Dimitriev, M. Semina, M. Nestoklon, _et al._ , Nat. Commun. 13, 1 (2022).
  * Lepine [1972] D. J. Lepine, Phys. Rev. B 6, 436 (1972).
  * Huang  _et al._ [2007] B. Huang, D. J. Monsma, and I. Appelbaum, Phys. Rev. Lett. 99, 177209 (2007).
  * Dey  _et al._ [2017] P. Dey, L. Yang, C. Robert, G. Wang, B. Urbaszek, X. Marie, and S. A. Crooker, Phys. Rev. Lett. 119, 137401 (2017).
  * Li  _et al._ [2021] J. Li, M. Goryca, K. Yumigeta, H. Li, S. Tongay, and S. A. Crooker, Phys. Rev. Mater. 5, 044001 (2021).
  * Goryca  _et al._ [2019] M. Goryca, N. P. Wilson, P. Dey, X. Xu, and S. A. Crooker, Sci. Adv. 5, eaau4899 (2019).
  * Song  _et al._ [2016] X. Song, S. Xie, K. Kang, J. Park, and V. Sih, Nano Lett. 16, 5010 (2016).
  * Yan  _et al._ [2017] T. Yan, S. Yang, D. Li, and X. Cui, Phys. Rev. B 95, 241406 (2017).
  * Edelberg  _et al._ [2019] D. Edelberg, D. Rhodes, A. Kerelsky, B. Kim, J. Wang, A. Zangiabadi, C. Kim, A. Abhinandan, J. Ardelean, M. Scully, _et al._ , Nano Lett. 19, 4371 (2019).
  * Rhodes  _et al._ [2019] D. Rhodes, S. H. Chae, R. Ribeiro P., and J. Hone, Nat. Mater. 18, 541 (2019).
  * Yankowitz  _et al._ [2015] M. Yankowitz, D. McKenzie, and B. J. LeRoy, Phys. Rev. Lett. 115, 136803 (2015).
  * Van Tuan  _et al._ [2016] D. Van Tuan, F. Ortmann, A. W. Cummings, D. Soriano, and S. Roche, Sci. Rep. 6, 1 (2016).
  * Ertler  _et al._ [2009] C. Ertler, S. Konschuh, M. Gmitra, and J. Fabian, Phys. Rev. B 80, 041405 (2009).
  * Tao and Tsymbal [2019] L. Tao and E. Y. Tsymbal, Phys. Rev. B 100, 161110 (2019).
  * Hilton and Tang [2002] D. J. Hilton and C. L. Tang, Phys. Rev. Lett. 89, 146601 (2002).
  * Jiang and Wu [2009] J. H. Jiang and M. W. Wu, Phys, Rev. B 79, 125206 (2009).
  * Kamra and Ghosh [2011] A. Kamra and B. Ghosh, J. Appl. Phys 109, 024501 (2011).
  * Dettwiler  _et al._ [2017] F. Dettwiler, J. Fu, S. Mack, P. J. Weigele, J. C. Egues, D. D. Awschalom, and D. M. Zumbühl, Phys. Rev. X 7, 031010 (2017).
  * Huber  _et al._ [2018] D. Huber, M. Reindl, S. F. C. Da Silva, C. Schimpf, J. Martín-Sánchez, H. Huang, G. Piredda, J. Edlinger, A. Rastelli, and R. Trotta, Phys. Rev. Lett. 121, 033902 (2018).
  * Belykh  _et al._ [2018] V. V. Belykh, A. Y. Kuntsevich, M. M. Glazov, K. V. Kavokin, D. R. Yakovlev, and M. Bayer, Phys. Rev. X 8, 031021 (2018).
  * Ohno  _et al._ [1999] Y. Ohno, R. Terauchi, T. Adachi, F. Matsukura, and H. Ohno, Phys. Rev. Lett. 83, 4196 (1999).
  * Hohage  _et al._ [2006] P. E. Hohage, G. Bacher, D. Reuter, and A. D. Wieck, Appl. Phys. Lett. 89, 231101 (2006).
  * Yu  _et al._ [2005] Z. G. Yu, S. Krishnamurthy, M. Van Schilfgaarde, and N. Newman, Phys, Rev. B 71, 245312 (2005).
  * Mower  _et al._ [2011] M. D. Mower, G. Vignale, and I. V. Tokatly, Phys, Rev. B 83, 155205 (2011).
  * Marchetti  _et al._ [2014] G. Marchetti, M. Hodgson, J. McHugh, R. Chantrell, and I. D’Amico, Materials 7, 2795 (2014).
  * Drogeler  _et al._ [2016] M. Drogeler, C. Franzen, F. Volmer, T. Pohlmann, L. Banszerus, M. Wolter, K. Watanabe, T. Taniguchi, C. Stampfer, and B. Beschoten, Nano Lett. 16, 3533 (2016).
  * Dean  _et al._ [2010] C. R. Dean, A. F. Young, I. Meric, C. Lee, L. Wang, S. Sorgenfrei, K. Watanabe, T. Taniguchi, P. Kim, K. L. Shepard, and J. Hone, Nat. Nanotechnol. 5, 722 (2010).
  * Han and Kawakami [2011] W. Han and R. K. Kawakami, Phys. Rev. Lett. 107, 047207 (2011).
  * Guimarães  _et al._ [2014] M. H. Guimarães, P. J. Zomer, J. Ingla-Aynés, J. C. Brant, N. Tombros, and B. J. van Wees, Phys. Rev. Lett. 113, 086602 (2014).
  * Raes  _et al._ [2016] B. Raes, J. E. Scheerder, M. V. Costache, F. Bonell, J. F. Sierra, J. Cuppens, J. Van de Vondel, and S. O. Valenzuela, Nat. Commun. 7, 11444 (2016).
  * Guimarães  _et al._ [2012] M. H. Guimarães, A. Veligura, P. Zomer, T. Maassen, I. Vera-Marun, N. Tombros, and B. Van Wees, Nano Letters 12, 3512 (2012).
  * Cummings  _et al._ [2017] A. W. Cummings, J. H. Garcia, J. Fabian, and S. Roche, Phys. Rev. Lett. 119, 206601 (2017).
  * Zomer  _et al._ [2012] P. J. Zomer, M. H. D. Guimarães, N. Tombros, and B. J. van Wees, Phys. Rev. B 86, 161416(R) (2012).
  * Kamalakar  _et al._ [2015] M. V. Kamalakar, C. Groenveld, A. Dankert, and S. P. Dash, Nat. Commun. 6, 1 (2015).
  * Zollner  _et al._ [2019] K. Zollner, M. Gmitra, and J. Fabian, Phys. Rev. B 99, 125151 (2019).
  * Paleari and Marini [2022] F. Paleari and A. Marini, Physical Review B 106, 125403 (2022).
  * Chan  _et al._ [2021] Y.-H. Chan, D. Y. Qiu, F. H. da Jornada, and S. G. Louie, Proceedings of the National Academy of Sciences 118, e1906938118 (2021).
  * Guo  _et al._ [2021] C. Guo, J. Xu, and Y. Ping, Journal of Physics: Condensed Matter 33, 234001 (2021).
  * Wu  _et al._ [2019] F. Wu, D. Rocca, and Y. Ping, J. Mater. Chem. C 7, 12891 (2019).
  * Smart  _et al._ [2021] T. J. Smart, K. Li, J. Xu, and Y. Ping, npj Computational Materials (2021).
  * Antonius and Louie [2022] G. Antonius and S. G. Louie, Phys. Rev. B 105, 085111 (2022).
  * Chen  _et al._ [2020] H.-Y. Chen, D. Sangalli, and M. Bernardi, Physical Review Letters 125, 107401 (2020).
  * Watanabe and Yanase [2021b] H. Watanabe and Y. Yanase, Phys. Rev. X 11, 011001 (2021b).
  * de Juan  _et al._ [2017] F. de Juan, A. G. Grushin, T. Morimoto, and J. E. Moore, Nat. Commun. 8, 15995 (2017).
  * Hellman  _et al._ [2013] O. Hellman, P. Steneteg, I. A. Abrikosov, and S. I. Simak, Physical Review B 87, 104111 (2013).
  * Macheda  _et al._ [2022] F. Macheda, P. Barone, and F. Mauri, Physical Review Letters 129, 185902 (2022).
  * Towns  _et al._ [2014] J. Towns, T. Cockerill, M. Dahan, I. Foster, K. Gaither, A. Grimshaw, V. Hazlewood, S. Lathrop, D. Lifka, G. D. Peterson, R. Roskies, J. R. Scott, and N. Wilkins-Diehr, Comput. Sci. Eng. 16, 62 (2014).