Many-body theory of phonon-induced spin relaxation and decoherence

Jinsoo Park Department of Applied Physics and Materials Science, California Institute of Technology, Pasadena, CA 91125, USA. Yao Luo Department of Applied Physics and Materials Science, California Institute of Technology, Pasadena, CA 91125, USA. Jin-Jian Zhou School of Physics, Beijing Institute of Technology, Beijing 100081, China. Marco Bernardi [email protected] Department of Applied Physics and Materials Science, California Institute of Technology, Pasadena, CA 91125, USA.

Abstract

First-principles calculations enable accurate predictions of electronic interactions and dynamics. However, computing the electron spin dynamics remains challenging. The spin-orbit interaction causes various dynamical phenomena that couple with phonons, such as spin precession and spin-flip $e$ -ph scattering, which are difficult to describe with current first-principles calculations. In this work, we show a rigorous framework to study phonon-induced spin relaxation and decoherence, by computing the spin-spin correlation function and its vertex corrections due to $e$ -ph interactions. We apply this approach to a model system and develop corresponding first-principles calculations of spin relaxation in GaAs. Our vertex-correction formalism is shown to capture the Elliott-Yafet, Dyakonov-Perel, and strong-precession mechanisms $-$ three independent spin decoherence regimes with distinct physical origins $–$ thereby unifying their theoretical treatment and calculation. Our method is general and enables quantitative studies of spin relaxation, decoherence, and transport in a wide range of materials and devices.

I Introduction

Linear response theory provides a microscopic understanding of the response of a system to external perturbations and computes the associated correlation functions [1, 2, 3, 4, 5, 6]. First-principles calculations of electronic interactions [7, 8, 9, 10, 11, 12, 13, 14] complement this formalism, enabling precise predictions of materials properties and transport coefficients without resorting to empirical models or fitting parameters. In this context, electron-phonon ( $e$ -ph) interactions are particularly important as they govern a wide range of phenomena such as charge transport [15], superconductivity [16], spin transport [17, 18, 19] and spin decoherence [20, 21, 22].
The Boltzmann transport equation (BTE) is widely used to study the response to an external electric field [23, 24]. The field drives the electronic populations $f_{n\bm{k}}$ , for states with band $n$ and crystal momentum $\bm{k}$ , away from the equilibrium Fermi-Dirac distribution $f^{0}_{n\bm{k}}$ , while the $e$ -ph interactions dissipate electron energy and act to restore equilibrium, resulting in a steady-state current $e(f_{n\bm{k}}-f^{0}_{n\bm{k}})v_{n\bm{k}}$ , where $e$ is the electron charge and $v_{n\bm{k}}$ is the band velocity [23]. In the many-body formalism, the BTE at low electric field is formally equivalent to the ladder vertex-correction to the dc conductivity [4]. In that framework, one determines the current-current correlation function, with vertex corrections from the $e$ -ph interactions obtained by summing over ladder diagrams, and computes the conductivity from the dissipative part of the susceptibility. A key factor making this approach equivalent to the BTE is that the electron velocity is band-diagonal in the Bloch basis, $\bra{m\bm{k}}\hat{v}\ket{n\bm{k}}=\delta_{nm}\partial_{\bm{k}}E_{n\bm{k}}/\hbar$ [4].
However, studying the response to an external field of an arbitrary operator that couples with phonons is more difficult. The matrix representation of an operator $\hat{A}$ is in general nondiagonal in the Bloch basis, $\hat{A}_{nm\bm{k}}\!=\!\bra{m\bm{k}}\hat{A}\ket{n\bm{k}}$ , and can mix states in different bands. The BTE cannot be applied in this case because due to its population-based formalism it neglects such off-diagonal (inter-band) components. A framework treating the response of non-diagonal operators coupled with $e$ -ph interactions is still missing.
An important example is spin relaxation and decoherence, where spin-orbit coupling (SOC) makes the spin operators non-diagonal in the band index, and phonons can change the electron spin through $e$ -ph interactions [19]. Theories of spin decoherence focus on two distinct models $-$ the Elliott-Yafet (EY) mechanism [25, 26], where $e$ -ph collisions rotate the spin direction, and the Dyakonov-Perel (DP) mechanism [27], where spin precession in the SOC field induces a motional narrowing of the spin. The dominant mechanism depends on the system $-$ typically, EY dominates in centrosymmetric and DP in non-centrosymmetric materials. Spin relaxation exhibits opposite trends in these two mechanisms, with spin relaxation times proportional to the $e$ -ph relaxation times in EY, and inversely proportional in DP. We have recently shown that EY spin relaxation can be computed from first-principles in the spin relaxation time approximation (sRTA) [28] $-$ the spin counterpart of the transport RTA for charge transport [1, 4] $-$ but spin precession and the DP mechanism are neglected in the sRTA.
Here we show a many-body approach to compute the susceptibility for an arbitrary non-diagonal operator coupled to $e$ -ph interactions. Our diagrammatic approach, based on the Kubo formula with vertex corrections to the susceptibility in an external injection field, calculates an effective phonon-dressed operator and its renormalized dynamics. We derive a Bethe-Salpeter equation (BSE) for the vertex corrections, and specializing to the spin operator, we use the vertex corrections to compute spin relaxation and precession. We show that the vertex corrections can capture spin decoherence due to both the EY and DP mechanisms and can also model the strong-precession regime, a third mechanism distinct from EY and DP. We find these three mechanisms in the exact solution of a two-level system, and also identify them in a real material, GaAs, using first-principles calculations. Combined with first-principles $e$ -ph calculations, our method is poised to advance microscopic understanding of phonon-induced spin decoherence [29], with applications ranging from solid-state qubits to quantum materials with spin Hall effect, valley-dependent spin physics, and Rashba effect.
The paper is organized as follows: In Sec. II, we derive the BSE for the phonon-dressed vertex, discuss its physical interpretation, and calculate the susceptibility in response to an injection field. In Sec. III $-$ IV, we apply this formalism to study spin dynamics in a model two-level system and in a real material, GaAs, discussing spin relaxation due to the EY, DP, and strong-precession mechanisms.

II Theory

We derive a self-consistent BSE for the vertex correction to the susceptibility due to $e$ -ph interactions, focusing on a general vector observable $\mathbf{\hat{A}}$ . We then present a physical interpretation of the vertex corrections and the renormalized dynamics of the operator. We employ atomic units and set $\hbar=1$ .

II.1 Interacting Green’s function

We consider an unperturbed Hamiltonian $H_{0}$ diagonal in the Bloch basis, $\bra{n^{\prime}\bm{k}}H_{0}\ket{n\bm{k}}=\varepsilon_{n\bm{k}}\delta_{nn^{\prime}}$ . The interacting imaginary-time Green’s function $\mathcal{G}(i\omega_{a})$ is written using the Dyson equation as [1]

\mathcal{G}(i\omega_{a})^{-1}=\mathcal{G}^{(0)}(i\omega_{a})^{-1}-\Sigma(i\omega_{a}),

(1)

where $\omega_{a}$ are fermionic Matsubara frequencies, $\mathcal{G}^{(0)}(i\omega_{a})$ is the non-interacting Green’s function, and $\Sigma(i\omega_{a})$ is the lowest order (Fan-Migdal) $e$ -ph self-energy [1, 15, 30], whose band- and $\bm{k}$ -dependent expression is

\begin{split}\Sigma_{nn^{\prime}\bm{k}}&(i\omega_{a})=-\frac{1}{\beta N_{q}V_{\text{uc}}}\sum_{mm^{\prime}\bm{q}\nu,iq_{c}}\left[g_{n^{\prime}m^{\prime}\nu}(\bm{k},\bm{q})\right]^{*}g_{nm\nu}(\bm{k},\bm{q})\\ &\times\mathcal{D}_{\nu\bm{q}}(iq_{c})\mathcal{G}_{mm^{\prime}\bm{k}+\bm{q}}(i\omega_{a}+iq_{c}).\end{split}

(2)

Here, $\beta=1/{k_{B}T}$ at temperature $T$ , $N_{q}$ is the number of $\bm{q}$ -points in the summation, $V_{\text{uc}}$ is the unit cell volume, $q_{c}$ is the bosonic Matsubara frequency of the phonon, and $\mathcal{D}_{\nu\bm{q}}(iq_{c})=2\omega_{\nu\bm{q}}/((iq_{c})^{2}-\omega_{\nu\bm{q}}^{2})$ is the non-interacting phonon Green’s function for a phonon with mode index $\nu$ , wave-vector $\bm{q}$ , and energy $\omega_{\nu\bm{q}}$ . The $e$ -ph matrix elements $g_{nm\nu}(\bm{k},\bm{q})$ quantify the probability amplitude for an electron in a Bloch state $\ket{\psi_{n\bm{k}}}$ , with band index $n$ and crystal momentum $\bm{k}$ , to scatter into a final state $\ket{\psi_{m\bm{k}+\bm{q}}}$ by emitting or absorbing a phonon [15, 31],

g_{nm\nu}(\bm{k},\bm{q})\!=\!\bra{\psi_{m\bm{k}+\bm{q}}}\partial_{\nu\bm{q}}\hat{V}\!\ket{\psi_{n\bm{k}}}\!,

(3)

where $\partial_{\nu\bm{q}}\hat{V}$ is the perturbation to the potential acting on an electron due to a given phonon mode $(\nu,\bm{q})$ .

Refer to caption — Figure 1: (a) Bare bubble diagram without the vertex correction. (b) Bubble diagram including the vertex correction. (c) Bethe-Salpeter equation for the vertex corrections $\Lambda$ from electron-phonon interactions within the ladder approximation. The wavy line is the phonon propagator and the red dots are the $e$ -ph matrix elements $g_{nm\nu}(\bm{k},\bm{q})$ .

II.2 Kubo formula and correlation function

We consider a complex vector operator $\mathbf{\hat{A}}$ , with matrix elements in the direction $\alpha$ written as $A^{\alpha}_{nm\bm{k}}=\bra{m\bm{k}}\hat{A}^{\alpha}\ket{n\bm{k}}$ . We derive the $\mathbf{\hat{A}}-\mathbf{\hat{A}}$ correlation function with a procedure analogous to the derivation of the dc conductivity in the ladder approximation [4]. Here, the operator $\mathbf{\hat{A}}$ is in general non-diagonal in the band index, leading to matrix elements $A^{\alpha}_{nm\bm{k}}$ , so the derivation for the diagonal case given in Ref. [4] needs to be extended to non-diagonal operators and vertex corrections.
We first derive the correlation function in imaginary time and frequency, and then extend it to real frequencies via analytic continuation. The retarded correlation function for the operator $\mathbf{\hat{A}}$ can be obtained from the Kubo formula [1]

\chi_{\alpha\beta}(\bm{p},i\nu_{b})=\int_{0}^{\beta}d\tau e^{i\nu_{b}\tau}\left<T_{\tau}\hat{A}^{\alpha}(\bm{p},\tau)\hat{A}^{\beta}(-\bm{p},0)\right>,

(4)

where $\bm{p}$ is a wave-vector, $\nu_{b}$ is a bosonic Matsubara frequency, $\tau$ is imaginary time ranging from $0$ to $\beta=1/{k_{B}T}$ at temperature $T$ , and $T_{\tau}$ is the imaginary time-ordering operator. Here we focus on the $\bm{p}\rightarrow 0$ limit, so we drop $\bm{p}$ from the equations. This correlation function can be expressed as a sum of bubble diagrams $P$ as [1]

\chi_{\alpha\beta}(i\nu_{b})=\frac{1}{\beta}\sum_{i\omega_{a}}P(i\omega_{a},i\omega_{a}+i\nu_{b}).

(5)

Let us consider the bare bubble diagram that includes the electron self-energy only in the electron propagator $\mathcal{G}$ , as shown in Fig. 1(a):

\chi_{\alpha\beta}(i\nu_{b})=\frac{1}{\beta V_{\text{uc}}}\!\sum_{i\omega_{a}}\Tr[\mathcal{G}(i\omega_{a})\hat{A}^{\alpha}\mathcal{G}(i\omega_{a}+i\nu_{b})\hat{A}^{\beta}\big{]},

(6)

where the trace is evaluated over the band and momentum indices. In this expression, the operator $\mathbf{\hat{A}}$ can be regarded as the bare vertex of the correlation function. For the velocity operator, Eq. (6) leads to the well-known Drude conductivity [4, 1].
In this work, the corrections to the vertex originate from the $e$ -ph interactions, which couple electronic states with different band and crystal momenta. Figure 1(b) shows the correlation function including the vertex correction $\Lambda$ ,

\begin{split}\chi_{\alpha\beta}(i\nu_{b})=\frac{1}{\beta V_{\text{uc}}}\!\sum_{i\omega_{a}}\Tr&\big{[}\mathcal{G}(i\omega_{a})\hat{A}^{\alpha}\mathcal{G}(i\omega_{a}+i\nu_{b})\\ &\times\hat{A}^{\beta}\Lambda^{\beta}(i\omega_{a},i\omega_{a}+i\nu_{b})\big{]},\end{split}

(7)

where $\hat{A}^{\beta}\Lambda^{\beta}(i\omega_{a},i\omega_{a}+i\nu_{b})$ is the phonon-dressed vertex for the operator $\hat{A}$ in the Cartesian direction $\beta$ . Note that the vertex correction ${\Lambda}^{\beta}(i\omega_{a},i\omega_{a}+i\nu_{b})$ is a complex-valued vector that contains information about the operator dynamics renormalized by the $e$ -ph interactions.

II.3 Bethe-Salpeter equation for the phonon-dressed vertex

The leading correction to the vertex is obtained by summing over ladder diagrams, which can be viewed as an abstract form of charge conservation in the presence of $e$ -ph scattering [1, 4]. The vertex correction $\Lambda^{\alpha}_{nn^{\prime}\bm{k}}$ satisfies the self-consistent BSE, shown diagrammatically in Fig. 1(c) and written as

\begin{split}A^{\alpha}_{nn^{\prime}\bm{k}}&\Lambda^{\alpha}_{nn^{\prime}\bm{k}}(i\omega_{a},i\omega_{a}+i\nu_{b})=A^{\alpha}_{nn^{\prime}\bm{k}}\\ -&\frac{1}{\beta N_{q}V_{\text{uc}}}\sum_{mm^{\prime}ll^{\prime}\bm{q}\nu,iq_{c}}\left[g_{n^{\prime}m^{\prime}\nu}(\bm{k},\bm{q})\right]^{*}g_{nm\nu}(\bm{k},\bm{q})\mathcal{D}_{\nu\bm{q}}(iq_{c})\\ &\times\mathcal{G}_{ml\bm{k}+\bm{q}}(i\omega_{a}+iq_{c})\mathcal{G}_{l^{\prime}m^{\prime}\bm{k}+\bm{q}}(i\omega_{a}+i\nu_{b}+iq_{c})\\ &\times A^{\alpha}_{ll^{\prime}\bm{k}}\Lambda^{\alpha}_{ll^{\prime}\bm{k}+\bm{q}}(i\omega_{a}+iq_{c},i\omega_{a}+i\nu_{b}+iq_{c}).\end{split}

(8)

The kernel of this BSE [32] is the $e$ -ph interaction $\left[g_{n^{\prime}m^{\prime}\nu}(\bm{k},\bm{q})\right]^{*}g_{nm\nu}(\bm{k},\bm{q})\mathcal{D}_{\nu\bm{q}}(iq_{c})$ .
Following Mahan [1] and Ref. [4], we first sum over the bosonic Matsubara frequency $iq_{c}$ in Eq. (8). This summation, defined as $S(i\omega_{a},i\omega_{a}+i\nu_{b})$ , reads:

\begin{split}S&(i\omega_{a},i\omega_{a}+i\nu_{b})=\sum_{ll^{\prime}}S_{ll^{\prime}}(i\omega_{a},i\omega_{a}+i\nu_{b})\\ =&\frac{1}{\beta}\sum_{ll^{\prime}iq_{c}}\,\mathcal{D}_{\nu\bm{q}}(iq_{c})\,\Lambda^{\alpha}_{ll^{\prime}\bm{k}+\bm{q}}(i\omega_{a}+iq_{c},i\omega_{a}+i\nu_{b}+iq_{c})\\ &\times\mathcal{G}_{ml\bm{k}+\bm{q}}(i\omega_{a}+iq_{c})\,\mathcal{G}_{l^{\prime}m^{\prime}\bm{k}+\bm{q}}(i\omega_{a}+i\nu_{b}+iq_{c}).\end{split}

(9)

As usual, the summation is done by constructing a contour integral along a circle at infinity:

\begin{split}\oint\frac{dz}{2\pi i}&n_{\rm B}(z)\mathcal{D}_{\nu\bm{q}}(z)\,\Lambda^{\alpha}_{ll^{\prime}\bm{k}+\bm{q}}(i\omega_{a}+z,i\omega_{a}+i\nu_{b}+z)\\ \times\,&\mathcal{G}_{ml\bm{k}+\bm{q}}(i\omega_{a}+z)\,\mathcal{G}_{l^{\prime}m^{\prime}\bm{k}+\bm{q}}(i\omega_{a}+i\nu_{b}+z),\end{split}

(10)

where $n_{\rm B}$ are Bose-Einstein occupations. The integrand has poles at $z=iq_{c}$ , $z=\pm\omega_{\nu\bm{q}}$ , and branch cuts along $z=-i\omega_{a}$ and $z=-i\omega_{a}-i\nu_{b}$ [1, 4]. Employing Cauchy’s residue theorem, we obtain

\begin{split}S_{ll^{\prime}}(i&\omega_{a},i\omega_{a}+i\nu_{b})\\ =-&N_{\nu\bm{q}}\Lambda^{\alpha}_{ll^{\prime}\bm{k}+\bm{q}}(i\omega_{a}+\omega_{\nu\bm{q}},i\omega_{a}+i\nu_{b}+\omega_{\nu\bm{q}})\\ &\times\mathcal{G}_{ml\bm{k}+\bm{q}}(i\omega_{a}+\omega_{\nu\bm{q}})\mathcal{G}_{l^{\prime}m^{\prime}\bm{k}+\bm{q}}(i\omega_{a}+i\nu_{b}+\omega_{\nu\bm{q}})\\ -&[N_{\nu\bm{q}}+1]\Lambda^{\alpha}_{ll^{\prime}\bm{k}+\bm{q}}(i\omega_{a}-\omega_{\nu\bm{q}},i\omega_{a}+i\nu_{b}-\omega_{\nu\bm{q}})\\ &\times\mathcal{G}_{ml\bm{k}+\bm{q}}(i\omega_{a}-\omega_{\nu\bm{q}})\mathcal{G}_{l^{\prime}m^{\prime}\bm{k}+\bm{q}}(i\omega_{a}+i\nu_{b}-\omega_{\nu\bm{q}})\\ -&\int\frac{d\varepsilon^{\prime}}{2\pi i}f(\varepsilon^{\prime})\frac{2\omega_{\nu\bm{q}}}{(\varepsilon^{\prime}-i\omega_{a})^{2}-\omega_{\nu\bm{q}}^{2}}\mathcal{G}_{l^{\prime}m^{\prime}\bm{k}+\bm{q}}(\varepsilon^{\prime}+i\nu_{b})\\ &\times[\Lambda^{\alpha}_{ll^{\prime}\bm{k}+\bm{q}}(\varepsilon^{\prime}+i\eta,\varepsilon^{\prime}+i\nu_{b})\mathcal{G}_{ml\bm{k}+\bm{q}}(\varepsilon^{\prime}+i\eta)\\ &-\Lambda^{\alpha}_{ll^{\prime}\bm{k}+\bm{q}}(\varepsilon^{\prime}-i\eta,\varepsilon^{\prime}+i\nu_{b})\mathcal{G}_{ml\bm{k}+\bm{q}}(\varepsilon^{\prime}-i\eta)]\\ -&\int\frac{d\varepsilon^{\prime}}{2\pi i}f(\varepsilon^{\prime})\frac{2\omega_{\nu\bm{q}}}{(\varepsilon^{\prime}-i\omega_{a}-i\nu_{b})^{2}-\omega_{\nu\bm{q}}^{2}}\mathcal{G}_{ml\bm{k}+\bm{q}}(\varepsilon^{\prime}-i\nu_{b})\\ &\times[\Lambda^{\alpha}_{ll^{\prime}\bm{k}+\bm{q}}(\varepsilon^{\prime}-i\nu_{b},\varepsilon^{\prime}+i\eta)\mathcal{G}_{l^{\prime}m^{\prime}\bm{k}+\bm{q}}(\varepsilon^{\prime}+i\eta)\\ &-\Lambda^{\alpha}_{ll^{\prime}\bm{k}+\bm{q}}(\varepsilon^{\prime}-i\nu_{b},\varepsilon^{\prime}-i\eta)\mathcal{G}_{l^{\prime}m^{\prime}\bm{k}+\bm{q}}(\varepsilon^{\prime}-i\eta)]\,,\end{split}

(11)

where $N_{\nu\bm{q}}\!=\!n_{\rm B}(\omega_{\nu\bm{q}})$ are temperature dependent phonon occupations, $f(\varepsilon)$ is the Fermi-Dirac distribution function, and $\eta$ is a positive infinitesimal.
The leading contribution to $S_{ll^{\prime}}(i\omega_{a},i\omega_{a}+i\nu_{b})$ comes from the combination of retarded and advanced Green’s functions, $G^{R}$ and $G^{A}$ , while terms of $O([G^{R}]^{2},[G^{A}]^{2})$ can be neglected at low electron density [1, 4]. Therefore, after the analytic continuations $i\omega_{a}\rightarrow\varepsilon-i\eta$ and $i\omega_{a}+i\nu_{b}\rightarrow\varepsilon+\nu+i\eta$ , and using the identity $\frac{1}{x+i\eta}=P\frac{1}{x}-i\pi\delta(x)$ , we obtain $S_{ll^{\prime}}(\varepsilon-i\eta,\varepsilon+i\eta)$ in limit of $\nu\rightarrow 0$ ,

\begin{split}S_{ll^{\prime}}(\varepsilon-&i\eta,\varepsilon+i\eta)\\ =-[&N_{\nu\bm{q}}+f(\varepsilon+\omega_{\nu\bm{q}})]\Lambda^{\alpha}_{ll^{\prime}\bm{k}+\bm{q}}(\varepsilon+\omega_{\nu\bm{q}})\\ &\times G_{ml\bm{k}+\bm{q}}^{R}(\varepsilon+\omega_{\nu\bm{q}})G_{l^{\prime}m^{\prime}\bm{k}+\bm{q}}^{A}(\varepsilon+\omega_{\nu\bm{q}})\\ -[&N_{\nu\bm{q}}+1-f(\varepsilon-\omega_{\nu\bm{q}})]\Lambda^{\alpha}_{ll^{\prime}\bm{k}+\bm{q}}(\varepsilon-\omega_{\nu\bm{q}})\\ &\times G_{ml\bm{k}+\bm{q}}^{R}(\varepsilon-\omega_{\nu\bm{q}})G_{l^{\prime}m^{\prime}\bm{k}+\bm{q}}^{A}(\varepsilon-\omega_{\nu\bm{q}}),\end{split}

(12)

where the index A (R) stands for advanced (retarded) function, and $\Lambda^{\alpha}(\varepsilon)\equiv\Lambda^{\alpha}(\varepsilon-i\eta,\varepsilon+i\eta)$ .
Using this result, we write the self-consistent BSE for the phonon-dressed vertex $\mathbf{\hat{A}}\Lambda$ at energy $\varepsilon$ as:

\begin{split}&A^{\alpha}_{nn^{\prime}\bm{k}}\Lambda^{\alpha}_{nn^{\prime}\bm{k}}(\varepsilon)=A^{\alpha}_{nn^{\prime}\bm{k}}\\ &+\frac{1}{N_{q}V_{\text{uc}}}\sum_{mm^{\prime}ll^{\prime}\bm{q}\nu}\left[g_{n^{\prime}m^{\prime}\nu}(\bm{k},\bm{q})\right]^{*}g_{nm\nu}(\bm{k},\bm{q})A^{\alpha}_{ll^{\prime}\bm{k}+\bm{q}}\\ &\times\bigg{[}(N_{\nu\bm{q}}+f(\varepsilon+\omega_{\nu\bm{q}}))\Lambda^{\alpha}_{ll^{\prime}\bm{k}+\bm{q}}(\varepsilon+\omega_{\nu\bm{q}})\\ &\times G_{ml\bm{k}+\bm{q}}^{R}(\varepsilon+\omega_{\nu\bm{q}})G_{l^{\prime}m^{\prime}\bm{k}+\bm{q}}^{A}(\varepsilon+\omega_{\nu\bm{q}})\\ &+(N_{\nu\bm{q}}+1-f(\varepsilon-\omega_{\nu\bm{q}}))\Lambda^{\alpha}_{ll^{\prime}\bm{k}+\bm{q}}(\varepsilon-\omega_{\nu\bm{q}})\\ &\times G_{ml\bm{k}+\bm{q}}^{R}(\varepsilon-\omega_{\nu\bm{q}})G_{l^{\prime}m^{\prime}\bm{k}+\bm{q}}^{A}(\varepsilon-\omega_{\nu\bm{q}})\bigg{]}.\end{split}

(13)

By solving Eq. (13), we obtain the phonon-dressed vertex $A^{\alpha}_{nn^{\prime}\bm{k}}\Lambda^{\alpha}_{nn^{\prime}\bm{k}}(\varepsilon)$ and its dependence on band, crystal momentum and energy.
In the weak scattering regime, where the electron spectral function has a well-defined quasiparticle peak [5] and the off-diagonal self-energy can be neglected [33, 34], the Green’s function becomes band-diagonal and the self-energies can be evaluated on-shell. Then the product of the retarded and advanced Green’s functions, $G^{R}G^{A}$ , can be approximated as [35]

\begin{split}&G^{R}_{m\bm{k}+\bm{q}}(\varepsilon)G^{A}_{m^{\prime}\bm{k}+\bm{q}}(\varepsilon)\\ &=\frac{G_{m^{\prime}\bm{k}+\bm{q}}^{A}(\varepsilon)-G_{m\bm{k}+\bm{q}}^{R}(\varepsilon)}{G_{m\bm{k}+\bm{q}}^{R}(\varepsilon)^{-1}-G_{m^{\prime}\bm{k}+\bm{q}}^{A}(\varepsilon)^{-1}}\\ &\approx\frac{\pi\delta(\varepsilon\!-\!\varepsilon_{m^{\prime}\bm{k}+\bm{q}})\!+\!\pi\delta(\varepsilon\!-\!\varepsilon_{m\bm{k}+\bm{q}})\!-\!iP\frac{1}{\varepsilon-\varepsilon_{m^{\prime}\bm{k}+\bm{q}}}\!+\!iP\frac{1}{\varepsilon-\varepsilon_{m\bm{k}+\bm{q}}}}{i(\Sigma^{R}_{m\bm{k}+\bm{q}}-\Sigma^{A}_{m^{\prime}\bm{k}+\bm{q}})+i(\varepsilon_{m\bm{k}+\bm{q}}-\varepsilon_{m^{\prime}\bm{k}+\bm{q}})},\end{split}

(14)

a function that is strongly peaked at electron energies $\varepsilon\!=\!\varepsilon_{m\bm{k}+\bm{q}}$ and $\varepsilon\!=\!\varepsilon_{m^{\prime}\bm{k}+\bm{q}}$ . Therefore, we can further simplify the full-frequency BSE in Eq. (13) to a double-pole ansatz, which evaluates the vertex corrections only at these two energies:

\begin{split}&A^{\alpha}_{nn^{\prime}\bm{k}}\Lambda^{\alpha}_{nn^{\prime}\bm{k}}(\varepsilon)=A^{\alpha}_{nn^{\prime}\bm{k}}+\frac{2\pi}{N_{q}V_{\text{uc}}}\sum_{mm^{\prime}\bm{q}\nu}\left[g_{n^{\prime}m^{\prime}\nu}(\bm{k},\bm{q})\right]^{*}g_{nm\nu}(\bm{k},\bm{q})\\ &\times\frac{1}{2}\bigg{[}\{(N_{\nu\bm{q}}+f_{m\bm{k}+\bm{q}})(\delta(\varepsilon+\omega_{\nu\bm{q}}-\varepsilon_{m\bm{k}+\bm{q}})-\frac{i}{\pi}P\frac{1}{\varepsilon+\omega_{\nu\bm{q}}-\varepsilon_{m^{\prime}\bm{k}+\bm{q}}})\\ &+(N_{\nu\bm{q}}+1-f_{m\bm{k}+\bm{q}})(\delta(\varepsilon-\omega_{\nu\bm{q}}-\varepsilon_{m\bm{k}+\bm{q}})-\frac{i}{\pi}P\frac{1}{\varepsilon-\omega_{\nu\bm{q}}-\varepsilon_{m^{\prime}\bm{k}+\bm{q}}})\}\times\frac{A^{\alpha}_{mm^{\prime}\bm{k}+\bm{q}}\Lambda^{\alpha}_{mm^{\prime}\bm{k}+\bm{q}}(\varepsilon_{m\bm{k}+\bm{q}})}{i(\Sigma^{R}_{m\bm{k}+\bm{q}}-\Sigma^{A}_{m^{\prime}\bm{k}+\bm{q}})+i(\varepsilon_{m\bm{k}+\bm{q}}-\varepsilon_{m^{\prime}\bm{k}+\bm{q}})}\\ &+\{(N_{\nu\bm{q}}+f_{m^{\prime}\bm{k}+\bm{q}})(\delta(\varepsilon+\omega_{\nu\bm{q}}-\varepsilon_{m^{\prime}\bm{k}+\bm{q}})+\frac{i}{\pi}P\frac{1}{\varepsilon+\omega_{\nu\bm{q}}-\varepsilon_{m\bm{k}+\bm{q}}})\\ &+(N_{\nu\bm{q}}+1-f_{m^{\prime}\bm{k}+\bm{q}})(\delta(\varepsilon-\omega_{\nu\bm{q}}-\varepsilon_{m^{\prime}\bm{k}+\bm{q}})+\frac{i}{\pi}P\frac{1}{\varepsilon-\omega_{\nu\bm{q}}-\varepsilon_{m\bm{k}+\bm{q}}})\}\times\frac{A^{\alpha}_{mm^{\prime}\bm{k}+\bm{q}}\Lambda^{\alpha}_{mm^{\prime}\bm{k}+\bm{q}}(\varepsilon_{m^{\prime}\bm{k}+\bm{q}})}{i(\Sigma^{R}_{m\bm{k}+\bm{q}}-\Sigma^{A}_{m^{\prime}\bm{k}+\bm{q}})+i(\varepsilon_{m\bm{k}+\bm{q}}-\varepsilon_{m^{\prime}\bm{k}+\bm{q}})}\bigg{]},\end{split}

(15)

where $\varepsilon$ equals $\varepsilon_{n\bm{k}}$ or $\varepsilon_{n^{\prime}\bm{k}}$ , and $f_{m\bm{k}+\bm{q}}\equiv f(\varepsilon_{m\bm{k}+\bm{q}})$ .
We have tested the consistency of this theory by deriving a Ward identity [36, 4, 1] relating the self-energy and vertex corrections (see Appendix A). This result guarantees that $e$ -ph diagrams are taken into account consistently in the self-energy and in our BSE.

II.4 The dressed vertex and its interpretation

Table 1: Summary of the formalism for charge transport and spin decoherence.

	Charge transport (Ref. [4])	Spin decoherence
Operator	$v_{n\bm{k}}$ (diagonal)	$s_{nm\bm{k}}$ (non-diagonal)
External field $\mathcal{F}$	Vector potential ( $\mathbf{A}$ )	Magnetic field ( $\mathbf{B}$ )
Injection field $\dot{\mathcal{F}}$	$\mathbf{E}(\nu)=-i\nu\mathbf{A}(\nu)$	$\dot{\mathbf{B}}(\nu)=-i\nu\mathbf{B}(\nu)$
Vertex correction $\Lambda$	$\Lambda_{n\bm{k}}^{\alpha}(\varepsilon_{n\bm{k}})$	$\Lambda_{nn^{\prime}\bm{k}}^{\alpha}(\varepsilon_{n\bm{k}}),\Lambda_{nn^{\prime}\bm{k}}^{\alpha}(\varepsilon_{n^{\prime}\bm{k}})$
Renormalized dynamics $\tau$ , $\omega$	$\tau_{n\bm{k}}^{(\text{tr})\alpha}=\tau^{\text{e-ph}}_{n\bm{k}}\Lambda_{n\bm{k}}^{\alpha}(\epsilon_{n\bm{k}})$	$\frac{1}{\frac{1}{\tau_{nn^{\prime}\bm{k}}^{\alpha}(\varepsilon_{n\bm{k}})}+i{\omega}_{nn^{\prime}\bm{k}}^{\alpha}(\varepsilon_{n\bm{k}})}=\frac{\Lambda^{\alpha}_{nn^{\prime}\bm{k}}(\varepsilon_{n\bm{k}})}{i(\Sigma^{R}_{n\bm{k}}-\Sigma^{A}_{n^{\prime}\bm{k}})+i(\varepsilon_{n\bm{k}}-\varepsilon_{n^{\prime}\bm{k}})}$

We focus on the dressed operator divided by the band energy difference, a key term in Eq. (15):

\frac{A^{\alpha}_{mm^{\prime}\bm{k}+\bm{q}}\Lambda^{\alpha}_{mm^{\prime}\bm{k}+\bm{q}}(\varepsilon_{m^{\prime}\bm{k}+\bm{q}})}{i(\Sigma^{R}_{m\bm{k}+\bm{q}}\!-\!\Sigma^{A}_{m^{\prime}\bm{k}+\bm{q}})\!+\!i(\varepsilon_{m\bm{k}+\bm{q}}\!-\!\varepsilon_{m^{\prime}\bm{k}+\bm{q}})}.

(16)

This ratio describes the renormalized dynamics associated with the operator $\mathbf{\hat{A}}$ in the presence of $e$ -ph interactions. This dynamics is obtained by dividing Eq. (16) by the bare operator expectation value $A^{\alpha}_{mm^{\prime}\bm{k}+\bm{q}}$ , obtaining

\frac{\Lambda^{\alpha}_{mm^{\prime}\bm{k}+\bm{q}}(\varepsilon_{m^{\prime}\bm{k}+\bm{q}})}{i(\Sigma^{R}_{m\bm{k}+\bm{q}}\!-\!\Sigma^{A}_{m^{\prime}\bm{k}+\bm{q}})\!+\!i(\varepsilon_{m\bm{k}+\bm{q}}\!-\!\varepsilon_{m^{\prime}\bm{k}+\bm{q}})}.

(17)

The physical meaning of this ratio can be understood by analyzing the simple case of the velocity operator. As the velocity operator is band-diagonal and satisfies $v^{\alpha}_{mm^{\prime}\bm{k}+\bm{q}}=v^{\alpha}_{m\bm{k}+\bm{q}}\delta_{mm^{\prime}}$ , the band energy difference in the denominator vanishes, so the denominator is purely real because $\Sigma^{A}_{m^{\prime}\bm{k}+\bm{q}}=(\Sigma^{R}_{m^{\prime}\bm{k}+\bm{q}})^{*}$ . Thus Eq. (16) for the velocity operator becomes

\frac{v^{\alpha}_{m\bm{k}+\bm{q}}\Lambda^{\alpha}_{mm\bm{k}+\bm{q}}(\varepsilon_{m\bm{k}+\bm{q}})}{i(\Sigma^{R}_{m\bm{k}+\bm{q}}-\Sigma^{A}_{m\bm{k}+\bm{q}})}=v^{\alpha}_{m\bm{k}+\bm{q}}\tau^{\text{e-ph}}_{m\bm{k}+\bm{q}}\Lambda^{\alpha}_{mm\bm{k}+\bm{q}}(\varepsilon_{m\bm{k}+\bm{q}}),

(18)

where we used $\tau^{\text{e-ph}}_{m\bm{k}+\bm{q}}=1/|2\imaginary\Sigma_{m\bm{k}+\bm{q}}|$ for the $e$ -ph collision time. This equation gives the renormalized $e$ -ph mean free path, and dividing by the bare velocity we obtain the renormalized relaxation time, also known as the transport relaxation time [4],

\tau_{m\bm{k}+\bm{q}}^{\alpha(\text{tr})}\equiv\tau^{\text{e-ph}}_{m\bm{k}+\bm{q}}\,\Lambda^{\alpha}_{mm\bm{k}+\bm{q}}(\varepsilon_{m\bm{k}+\bm{q}})=\frac{\Lambda^{\alpha}_{mm\bm{k}+\bm{q}}(\varepsilon_{m\bm{k}+\bm{q}})}{{i(\Sigma^{R}_{m\bm{k}+\bm{q}}-\Sigma^{A}_{m\bm{k}+\bm{q}}})}.

(19)

For a non-diagonal operator, both the vertex correction and the operator expectation value are complex, so the ratio in Eq. (17) cannot be represented by a single real quantity with units of time as in Eq. (19). To extend the vertex correction to non-diagonal operators, we generalize this formalism by defining the renormalized microscopic relaxation times $\tau_{mm^{\prime}\bm{k}+\bm{q}}^{\alpha}(\varepsilon)$ and introducing the precession frequencies $\omega_{mm^{\prime}\bm{k}+\bm{q}}^{\alpha}(\varepsilon)$ :

\begin{split}&\frac{1}{\frac{1}{\tau_{mm^{\prime}\bm{k}+\bm{q}}^{\alpha}(\varepsilon)}+i{\omega}_{mm^{\prime}\bm{k}+\bm{q}}^{\alpha}(\varepsilon)}\\ &\equiv\frac{\Lambda^{\alpha}_{mm^{\prime}\bm{k}+\bm{q}}(\varepsilon)}{i(\Sigma^{R}_{m\bm{k}+\bm{q}}\!-\!\Sigma^{A}_{m^{\prime}\bm{k}+\bm{q}})\!+\!i(\varepsilon_{m\bm{k}+\bm{q}}\!-\!\varepsilon_{m^{\prime}\bm{k}+\bm{q}})},\end{split}

(20)

where $\varepsilon$ equals $\varepsilon_{m\bm{k}+\bm{q}}$ or $\varepsilon_{m^{\prime}\bm{k}+\bm{q}}$ . This way, without the vertex correction, the renormalized relaxation time reduces to the (non-diagonal) $e$ -ph collision time, $\tau^{\text{e-ph}}_{mm^{\prime}\bm{k}+\bm{q}}=1/|\imaginary\Sigma_{m\bm{k}+\bm{q}}+\imaginary\Sigma_{m^{\prime}\bm{k}+\bm{q}}|$ , and the renormalized precession frequency reduces to the bare operator rotation frequency, $\omega_{\textrm{B}}=(\varepsilon_{m\bm{k}+\bm{q}}+\real\Sigma_{m\bm{k}+\bm{q}})-(\varepsilon_{m^{\prime}\bm{k}+\bm{q}}+\real\Sigma_{m^{\prime}\bm{k}+\bm{q}})$ , with $A^{\alpha}_{mm^{\prime}\bm{k}+\bm{q}}(t)\propto e^{i\,{\omega}_{\textrm{B}}t}$ .

II.5 Vertex correction to the susceptibility

We derive the vertex-corrected susceptibility in response to an external field for the generic observable $\mathbf{\hat{A}}$ . Suppose that the complex operator $\hat{A}^{\alpha}$ couples to a vector field $\mathcal{F}^{\alpha}$ , with perturbation Hamiltonian $H^{\prime}=-\hat{A}^{\alpha}\mathcal{F}^{\alpha}$ . The susceptibility is defined as the response function in

\langle\hat{A}^{\alpha}(\nu)\rangle=\chi_{\alpha\beta}(\nu)\mathcal{F}^{\beta}(\nu),

(21)

where $\mathcal{F}$ is the external field along the direction $\beta$ , and $\langle\hat{A}^{\alpha}(\nu)\rangle$ is the response of the system along $\alpha$ at frequency $\nu$ due to the applied field.
To study relaxation and dissipation, we rewrite the response of the system as

\langle\hat{A}^{\alpha}(\nu)\rangle=\sigma_{\alpha\beta}(\nu)\dot{\mathcal{F}}^{\beta}(\nu),

(22)

thus expressing it in terms of the susceptibility $\sigma_{\alpha\beta}$ to the “injection field” at frequency $\nu$ , and $\dot{\mathcal{F}}^{\beta}(\nu)=-i\nu\mathcal{F}^{\beta}(\nu)$ . The injection field produces a nonequilibrium electron distribution with an injection rate equal to the inverse relaxation time of $\mathbf{\hat{A}}$ [37]. From Eqs. (21)-(22), we obtain

\sigma_{\alpha\beta}(\nu)=\frac{\chi_{\alpha\beta}(\nu)}{-i\nu}.

(23)

When $\mathcal{F}$ is the vector potential $\mathbf{A}$ , the injection field becomes the electric field $\mathbf{E}(\nu)=-i\nu\mathbf{A}(\nu)$ , the observable of interest is the current operator $A^{\alpha}_{nm\bm{k}}=e\delta_{nm}v^{\alpha}_{n\bm{k}}$ , and $\sigma_{\alpha\beta}(\nu)$ is the frequency-dependent conductivity tensor. When $\mathcal{F}$ is the magnetic field $\mathbf{B}$ , the injection field is its time derivative, $\dot{\mathbf{B}}(\nu)=-i\nu\mathbf{B}(\nu)$ , and the observable is the electron magnetic moment $A^{\alpha}_{nm\bm{k}}=g\mu_{B}s^{\alpha}_{nm\bm{k}}$ , which is proportional to the spin matrix $s^{\alpha}_{nm\bm{k}}$ [37, 38]. These results are summarized in Table 1.
We write the correlation function with vertex correction [see Eq. (7)] as a contour integral along a circle at infinity [1, 4],

\begin{split}\chi_{\alpha\beta}(i\nu_{b})=&-\frac{1}{V_{\text{uc}}}\oint\frac{dz}{2\pi i}f(z)\Tr[\\ \hfil&\mathcal{G}(z)\hat{A}^{\alpha}\mathcal{G}(z+i\nu_{b})\hat{A}^{\beta}\Lambda^{\beta}(z,z+i\nu_{b})\big{]},\end{split}

(24)

which has branch cuts along $z=-i\nu_{b}$ and $z=0$ , and poles at $z=i\omega_{a}$ , and thus

\begin{split}&\chi_{\alpha\beta}(i\nu_{b})=\frac{1}{V_{\text{uc}}}\int\frac{d\varepsilon}{2\pi i}f(\varepsilon)\Tr[\\ \hfil&-\mathcal{G}(\varepsilon+i\eta)\hat{A}^{\alpha}\mathcal{G}(\varepsilon+i\nu_{b})\hat{A}^{\beta}\Lambda^{\beta}(\varepsilon+i\eta,\varepsilon+i\nu_{b})\\ &+\mathcal{G}(\varepsilon-i\eta)\hat{A}^{\alpha}\mathcal{G}(\varepsilon+i\nu_{b})\hat{A}^{\beta}\Lambda^{\beta}(\varepsilon-i\eta,\varepsilon+i\nu_{b})\\ &-\mathcal{G}(\varepsilon-i\nu_{b})\hat{A}^{\alpha}\mathcal{G}(\varepsilon+i\eta)\hat{A}^{\beta}\Lambda^{\beta}(\varepsilon-i\nu_{b},\varepsilon+i\eta)\\ &+\mathcal{G}(\varepsilon-i\nu_{b})\hat{A}^{\alpha}\mathcal{G}(\varepsilon-i\eta)\hat{A}^{\beta}\Lambda^{\beta}(\varepsilon-i\nu_{b},\varepsilon-i\eta)\bigg{]}.\end{split}

(25)

After the analytic continuation $i\nu_{b}\rightarrow\nu+i\eta$ , we obtain the retarded correlation function to leading order by neglecting the terms $G^{R}G^{R}$ and $G^{A}G^{A}$ [1, 4]:

\begin{split}&\chi_{\alpha\beta}(\nu)\\ &=\frac{1}{V_{\text{uc}}}\int\frac{d\varepsilon}{2\pi i}(f(\varepsilon)-f(\varepsilon+\nu))\Tr[\\ \hfil&~{}~{}~{}~{}G^{R}(\varepsilon)\hat{A}^{\alpha}G^{A}(\varepsilon+\nu)\hat{A}^{\beta}\Lambda^{\beta}(i\omega_{a}-i\eta,i\omega_{a}+i\nu_{b}+i\eta)\big{]}\\ &\approx\frac{1}{N_{k}V_{\text{uc}}}\sum_{nm\bm{k}}\int\frac{d\varepsilon}{2\pi i}(f(\varepsilon)-f(\varepsilon+\nu))A^{\alpha}_{nm\bm{k}}A^{\beta}_{mn\bm{k}}\\ &\times\!\Lambda^{\beta}_{mn\bm{k}}(\varepsilon\!-\!i\eta,\varepsilon\!+\!\nu\!+\!i\eta)\frac{\pi\delta(\varepsilon\!-\!\varepsilon_{n\bm{k}})\!+\!\pi\delta(\varepsilon\!+\!\nu\!-\!\varepsilon_{m\bm{k}})}{i(\Sigma^{R}_{m\bm{k}}\!-\!\Sigma^{A}_{n\bm{k}})\!+\!i(\varepsilon_{m\bm{k}}\!+\!\nu\!-\!\varepsilon_{n\bm{k}})},\end{split}

(26)

where $N_{k}$ is the number of $\bm{k}$ -points, and we used Eq. (14) in the last equality. This equation characterizes the frequency-dependent response of the system [39].
We focus on the dc limit $\nu\rightarrow 0$ , where the driving field is static. The susceptibility with respect to the injection field becomes

\begin{split}&\lim_{\nu\rightarrow 0}\sigma_{\alpha\beta}(\nu)=-\lim_{\nu\rightarrow 0}\frac{1}{\nu}\imaginary\chi_{\alpha\beta}(\nu)\\ &=\frac{1}{N_{k}V_{\text{uc}}}\real\sum_{nm\bm{k}}A^{\alpha}_{nm\bm{k}}A^{\beta}_{mn\bm{k}}\\ &~{}~{}~{}~{}~{}\times\frac{\frac{1}{2}[(-\frac{df_{n\bm{k}}}{d\varepsilon})\Lambda^{\beta}_{mn\bm{k}}(\varepsilon_{n\bm{k}})+(-\frac{df_{m\bm{k}}}{d\varepsilon})\Lambda^{\beta}_{mn\bm{k}}(\varepsilon_{m\bm{k}})]}{{i(\Sigma^{R}_{m\bm{k}}-\Sigma^{A}_{n\bm{k}})+i(\varepsilon_{m\bm{k}}-\varepsilon_{n\bm{k}})}}.\end{split}

(27)

This static susceptibility has both band-diagonal ( $n=m$ ) and off-diagonal ( $n\neq m$ ) contributions. The band-diagonal contribution

\sigma_{\alpha\beta}^{\rm(d)}(0)=\frac{1}{N_{k}V_{\text{uc}}}\sum_{n\bm{k}}A^{\alpha}_{nn\bm{k}}A^{\beta}_{nn\bm{k}}\tau^{\text{e-ph}}_{n\bm{k}}\Lambda^{\beta}_{nn\bm{k}}(\varepsilon_{n\bm{k}})(-\frac{df_{n\bm{k}}}{d\varepsilon})

(28)

is the only contribution to the static susceptibility for a band-diagonal operator. For example, for the velocity operator, Eq. (28) becomes the well-known electrical conductivity tensor within the BTE, which has been studied extensively using both empirical and first-principles calculations [23, 24, 40, 41, 31] (note that solving exactly the BTE is equivalent to computing the velocity vertex correction [4]). For non-diagonal operators, our formalism introduces an off-diagonal contribution in Eq. (27):

\begin{split}&\sigma_{\alpha\beta}^{\rm(nd)}(0)=\frac{1}{N_{\bm{k}}V_{\text{uc}}}\real\sum_{n\neq m\bm{k}}A^{\alpha}_{nm\bm{k}}A^{\beta}_{mn\bm{k}}\\ &~{}~{}~{}~{}~{}\times\frac{\frac{1}{2}[(-\frac{df_{n\bm{k}}}{d\varepsilon})\Lambda^{\beta}_{mn\bm{k}}(\varepsilon_{n\bm{k}})+(-\frac{df_{m\bm{k}}}{d\varepsilon})\Lambda^{\beta}_{mn\bm{k}}(\varepsilon_{m\bm{k}})]}{{i(\Sigma^{R}_{m\bm{k}}-\Sigma^{A}_{n\bm{k}})+i(\varepsilon_{m\bm{k}}-\varepsilon_{n\bm{k}})}}.\end{split}

(29)

For the velocity operator, this term enables studies of charge transport in the presence of inter-band coherence [42]; for the spin operator, this contribution is essential to describe how $e$ -ph interactions modify spin precession.

II.6 Renormalized relaxation time

We derive an expression for the renormalized macroscopic relaxation time of an operator $\mathbf{\hat{A}}$ due to $e$ -ph interactions. Using Eq. (28), the average relaxation time for the band-diagonal components $A^{\alpha}_{nn\bm{k}}$ is

\begin{split}\tau_{\alpha\beta}=\frac{\sum_{n\bm{k}}A^{\alpha}_{nn\bm{k}}A^{\beta}_{nn\bm{k}}\tau^{\text{e-ph}}_{n\bm{k}}\Lambda^{\beta}_{nn\bm{k}}(\varepsilon_{n\bm{k}})(-\frac{df_{n\bm{k}}}{d\varepsilon})}{\sum_{n\bm{k}}A^{\alpha}_{nn\bm{k}}A^{\beta}_{nn\bm{k}}\big{(}-\frac{df_{n\bm{k}}}{d\varepsilon}\big{)}}.\end{split}

(30)

For the velocity operator, this equation gives the well-known Drude dc electrical conductivity, while for the spin operator one obtains the phonon-dressed macroscopic spin relaxation time, as discussed below. These results generalize the linear response treatment for band-diagonal operators presented in Ref. [4] and extend it to non-diagonal operators.

III Spin relaxation and decoherence

We now specialize to the non-diagonal spin operator, and apply our formalism to study phonon-induced spin relaxation and decoherence. The BSE for the phonon-dressed spin vertex $-$ called hereafter spin-phonon BSE $-$ is a key result obtained from Eq. (13) by replacing $A^{\alpha}_{nn^{\prime}\bm{k}}$ with the spin operator $s^{\alpha}_{nn^{\prime}\bm{k}}$ . In matrix form and using a compact notation, the spin-phonon BSE can be written in a way that clearly matches the diagram in Fig. 1(c):

\begin{split}\bm{s}\bm{\Lambda}_{\bm{k}}(\varepsilon)\!=\!\bm{s}_{\bm{k}}+\frac{1}{N_{q}V_{\text{uc}}}\!\sum_{\nu\bm{q}\pm}&\textbf{g}_{\nu\bm{k}\bm{q}}^{\dagger}\!\left[{G}^{A}\bm{s}\bm{\Lambda}{G}^{R}\right]_{\!\!\begin{smallmatrix}\bm{k}+\bm{q},~{}~{}\\ \varepsilon\pm\omega_{\nu\bm{q}}\end{smallmatrix}}\!\!\!\textbf{g}_{\nu\bm{k}\bm{q}}\,F_{\pm}(T),\end{split}

(31)

where $\bm{s}\bm{\Lambda}_{\bm{k}}(\varepsilon)=\bm{s}_{nn^{\prime}\bm{k}}\bm{\Lambda}_{nn^{\prime}\bm{k}}(\varepsilon)$ is the phonon-dressed spin vertex, $F_{\pm}(T)=N_{\nu\bm{q}}+\frac{1}{2}\pm[f(\varepsilon\pm\omega_{\nu\bm{q}})-\frac{1}{2}]$ is a thermal occupation factor at temperature $T$ , and $\left[\mathbf{g}_{\nu\bm{k}\bm{q}}\right]_{nm}=g_{nm\nu}(\bm{k},\bm{q})$ are $e$ -ph matrix elements [31].
In the weak scattering regime, this spin-phonon BSE can be rewritten using the double-pole ansatz discussed above (where $\varepsilon$ equals $\varepsilon_{n\bm{k}}$ or $\varepsilon_{n^{\prime}\bm{k}}$ ):

\begin{split}&s^{\alpha}_{nn^{\prime}\bm{k}}\Lambda^{\alpha}_{nn^{\prime}\bm{k}}(\varepsilon)=s^{\alpha}_{nn^{\prime}\bm{k}}+\frac{2\pi}{N_{q}V_{\text{uc}}}\sum_{mm^{\prime}\bm{q}\nu}\left[g_{n^{\prime}m^{\prime}\nu}(\bm{k},\bm{q})\right]^{*}g_{nm\nu}(\bm{k},\bm{q})\\ &\times\frac{1}{2}\bigg{[}\{(N_{\nu\bm{q}}+f_{m\bm{k}+\bm{q}})(\delta(\varepsilon+\omega_{\nu\bm{q}}-\varepsilon_{m\bm{k}+\bm{q}})-\frac{i}{\pi}P\frac{1}{\varepsilon+\omega_{\nu\bm{q}}-\varepsilon_{m^{\prime}\bm{k}+\bm{q}}})\\ &+(N_{\nu\bm{q}}+1-f_{m\bm{k}+\bm{q}})(\delta(\varepsilon-\omega_{\nu\bm{q}}-\varepsilon_{m\bm{k}+\bm{q}})-\frac{i}{\pi}P\frac{1}{\varepsilon-\omega_{\nu\bm{q}}-\varepsilon_{m^{\prime}\bm{k}+\bm{q}}})\}\times\frac{s^{\alpha}_{mm^{\prime}\bm{k}+\bm{q}}\Lambda^{\alpha}_{mm^{\prime}\bm{k}+\bm{q}}(\varepsilon_{m\bm{k}+\bm{q}})}{i(\Sigma^{R}_{m\bm{k}+\bm{q}}-\Sigma^{A}_{m^{\prime}\bm{k}+\bm{q}})+i(\varepsilon_{m\bm{k}+\bm{q}}-\varepsilon_{m^{\prime}\bm{k}+\bm{q}})}\\ &+\{(N_{\nu\bm{q}}+f_{m^{\prime}\bm{k}+\bm{q}})(\delta(\varepsilon+\omega_{\nu\bm{q}}-\varepsilon_{m^{\prime}\bm{k}+\bm{q}})+\frac{i}{\pi}P\frac{1}{\varepsilon+\omega_{\nu\bm{q}}-\varepsilon_{m\bm{k}+\bm{q}}})\\ &+(N_{\nu\bm{q}}+1-f_{m^{\prime}\bm{k}+\bm{q}})(\delta(\varepsilon-\omega_{\nu\bm{q}}-\varepsilon_{m^{\prime}\bm{k}+\bm{q}})+\frac{i}{\pi}P\frac{1}{\varepsilon-\omega_{\nu\bm{q}}-\varepsilon_{m\bm{k}+\bm{q}}})\}\times\frac{s^{\alpha}_{mm^{\prime}\bm{k}+\bm{q}}\Lambda^{\alpha}_{mm^{\prime}\bm{k}+\bm{q}}(\varepsilon_{m^{\prime}\bm{k}+\bm{q}})}{i(\Sigma^{R}_{m\bm{k}+\bm{q}}-\Sigma^{A}_{m^{\prime}\bm{k}+\bm{q}})+i(\varepsilon_{m\bm{k}+\bm{q}}-\varepsilon_{m^{\prime}\bm{k}+\bm{q}})}\bigg{]}.\end{split}

(32)

This BSE for the phonon-dressed spin vertex, used in this work to study spin dynamics, should not be confused with the widely used BSE for excitons and optical spectra [2], which is entirely unrelated.
The vertex corrections $\Lambda^{\alpha}_{nn^{\prime}\bm{k}}$ obtained by solving the BSE govern spin dynamics as they renormalize spin relaxation and precession [29]. The macroscopic spin relaxation times are obtained using the thermal average in Eq. (30),

\tau^{(s)}_{\alpha\beta}=\frac{\sum_{n\bm{k}}s^{\alpha}_{nn\bm{k}}s^{\beta}_{nn\bm{k}}\tau^{\text{e-ph}}_{n\bm{k}}\Lambda^{\beta}_{nn\bm{k}}(\varepsilon_{n\bm{k}})(-\frac{df_{n\bm{k}}}{d\varepsilon})}{\sum_{n\bm{k}}s^{\alpha}_{nn\bm{k}}s^{\beta}_{nn\bm{k}}\,\big{(}\!-\frac{df_{n\bm{k}}}{d\varepsilon}\big{)}}.

(33)

For $\alpha=\beta$ along the external magnetic field, Eq. (33) gives the longitudinal spin relaxation time, usually called $T_{1}$ , along the direction $\alpha$ , while for a perpendicular magnetic field one obtains the transverse spin relaxation time, $T_{2}$ [43]. The renormalized microscopic spin relaxation times ( $\tau_{nn^{\prime}\bm{k}}^{\alpha}$ ) and spin precession rates ( $\omega_{nn^{\prime}\bm{k}}^{\alpha}$ ), which are matrices in Bloch basis, are computed from the vertex corrections $\Lambda^{\alpha}_{nn^{\prime}\bm{k}}$ using

\frac{1}{\frac{1}{\tau_{nn^{\prime}\bm{k}}^{\alpha}(\varepsilon)}+i{\omega}_{nn^{\prime}\bm{k}}^{\alpha}(\varepsilon)}\equiv\frac{\Lambda^{\alpha}_{nn^{\prime}\bm{k}}(\varepsilon)}{i(\Sigma^{R}_{n\bm{k}}-\Sigma^{A}_{n^{\prime}\bm{k}})+i(\varepsilon_{n\bm{k}}-\varepsilon_{n^{\prime}\bm{k}})}.

(34)

The diagonal components with $n\!=\!n^{\prime}$ give the renormalized microscopic spin relaxation times, $\tau^{\beta}_{nn\bm{k}}=\tau^{\text{e-ph}}_{n\bm{k}}\,\Lambda^{\beta}_{nn\bm{k}}(\varepsilon_{n\bm{k}})$ , entering Eq. (33).

IV Results

We apply our formalism to study spin relaxation and decoherence. We first present analytic results for a two-level model system, and then focus on first-principles calculations on a real material, GaAs. Application to a wider range of materials is presented in our companion paper [29].

IV.1 Two-level system with optical phonon scattering

We study spin dynamics in a two-level system to understand different phonon-induced spin relaxation mechanisms. In our model, the electron spins undergo phonon-induced spin-flip transitions together with spin precession in the SOC field modified by the $e$ -ph interactions. We solve the spin-phonon BSE for this system and derive analytic expressions for the vertex corrections and spin relaxation times. Our analysis sheds light on phonon-dressed operators and their renormalized dynamics, providing a starting point to understand phonon-induced spin relaxation in real materials with complex band structures, phonon dispersions, and $e$ -ph interactions.
Consider a periodic two-level system where each level is spin degenerate (see Fig. 2). The Hilbert space consists of four Bloch states, which are eigenstates of the Hamiltonian:

H\!\ket{n}=\varepsilon_{n}\!\ket{n}

(35)

where $\ket{n}$ is the $n$ -th energy eigenstate. The two lowest-energy states $\ket{1}$ and $\ket{2}$ are degenerate ( $\varepsilon_{1}=\varepsilon_{2}$ ) and differ only in their spin part. The other two states, $\ket{3}$ and $\ket{4}$ , are higher in energy by $\omega_{O}$ and are perturbed by an internal magnetic field along $\hat{x}$ due to SOC. This field causes a small Zeeman splitting, $\Delta\ll\omega_{O}$ , such that $\varepsilon_{3,4}=\varepsilon_{1}+\omega_{O}\pm\frac{\Delta}{2}$ .
We separate the space-dependent part $\ket{\psi_{n}}$ and the spin-dependent part $\ket{\chi_{n}}$ of the two eigenstates as $\ket{n}=\ket{\psi_{n}}\otimes\ket{\chi_{n}}$ . The two lowest states have an identical space-dependent part $\ket{\psi}$ and are spin polarized along $\hat{z}$ :

\ket{1}=\ket{\psi}\otimes\left(\begin{matrix}1\\ 0\end{matrix}\right),~{}\ket{2}=\ket{\psi}\otimes\left(\begin{matrix}0\\ 1\end{matrix}\right),

(36)

with the following spin matrix elements along $z$ :

\begin{split}&\bra{1}\hat{s}_{z}\ket{1}=-\bra{2}\hat{s}_{z}\ket{2}=\frac{1}{2},\\ &\bra{1}\hat{s}_{z}\ket{2}=0.\end{split}

(37)

Above, $\hat{s}=\hat{\sigma}/2$ is the spin operator and $\hat{\sigma}$ are Pauli matrices. The two upper bands have an identical space-dependent part $\ket{\phi}$ and are spin polarized along $\hat{x}$ :

\ket{3}=\ket{\phi}\otimes\frac{1}{\sqrt{2}}\left(\begin{matrix}1\\ 1\end{matrix}\right),~{}\ket{4}=\ket{\phi}\otimes\frac{1}{\sqrt{2}}\left(\begin{matrix}1\\ -1\end{matrix}\right),

(38)

with spin matrix elements

\begin{split}&\bra{3}\hat{s}_{z}\ket{3}=\bra{4}\hat{s}_{z}\ket{4}=0,\\ &\bra{3}\hat{s}_{z}\ket{4}=\frac{1}{2}.\end{split}

(39)

The space-dependent part of the two lower states is orthogonal to that of the upper states, and the spin matrix elements between the two sets of states are zero.
In our model, an electron can scatter between the lower and upper levels by emitting or absorbing an optical phonon. These transitions are associated with $e$ -ph matrix elements $g_{nm}=\bra{m}\Delta\hat{V}\ket{n}$ , where $\Delta\hat{V}$ is the perturbation potential due to the optical phonon. We assume that this perturbation potential has the form

\Delta\hat{V}=\Delta\hat{V}(\bm{r})\otimes\left(\begin{matrix}a&b\\ b&a\end{matrix}\right),

(40)

where $\hat{V}(\bm{r})$ is the space-dependent part and $\left(\begin{smallmatrix}a&b\\ b&a\end{smallmatrix}\right)$ the spin-dependent part of the perturbation, with $a$ and $b$ real numbers. Due to the presence of SOC, the spin-dependent part $\left(\begin{smallmatrix}a&b\\ b&a\end{smallmatrix}\right)$ is different from the identity matrix. We consider a small spin-mixing $b$ , where $a^{2}+b^{2}=1$ , so that each phonon collision has a small probability $b^{2}\ll 1$ to flip the $z$ -component of the spin [25]. The $e$ -ph matrix elements become $g_{nm}=g_{0}\bra{\chi_{m}}\left(\begin{smallmatrix}a&b\\ b&a\end{smallmatrix}\right)\ket{\chi_{n}}$ , where $g_{0}=\bra{\phi}\Delta\hat{V}(\bm{r})\ket{\psi}$ . Thus the spin-dependent $e$ -ph matrix elements are

\begin{split}g_{13}=g_{23}=\frac{1}{\sqrt{2}}g_{0}(a+b),\\ g_{14}=-g_{24}=\frac{1}{\sqrt{2}}g_{0}(a-b).\end{split}

(41)

IV.2 Spin-phonon BSE and spin relaxation times

We construct the spin-phonon BSE for our two-level system. There are four non-zero spin matrix elements in our model: $s^{z}_{11}$ , $s^{z}_{22}$ , $s^{z}_{34}$ , and $s^{z}_{43}$ , where $s^{z}_{nm}=\bra{m}\hat{s}_{z}\ket{n}$ . Using Eq. (32), the diagonal matrix elements have only one vertex correction, while the off-diagonal matrix elements have two energy-dependent vertex corrections. Assuming a small Zeeman splitting $\Delta$ , we can neglect the energy dependence of the off-diagonal vertex corrections, and define $\Lambda^{z}_{nm}=\Lambda^{z}_{nm}(\varepsilon_{n})$ . We thus have four unknown vertex corrections: $\Lambda^{z}_{11}$ , $\Lambda^{z}_{22}$ , $\Lambda^{z}_{34}$ , and $\Lambda^{z}_{43}$ . Taking the $e$ -ph scattering rates of the two lower and upper states to be $\Gamma_{1}$ and $\Gamma_{3}$ , respectively, the spin-phonon BSE becomes

\begin{split}&\Lambda^{z}_{11}=1+\left(\frac{1}{2}-b^{2}\right)\left(\frac{\Gamma_{1}\Lambda^{z}_{34}}{-i\Delta+\Gamma_{3}}+\frac{\Gamma_{1}\Lambda^{z}_{43}}{i\Delta+\Gamma_{3}}\right)\\ -&\Lambda^{z}_{22}=-1-\left(\frac{1}{2}-b^{2}\right)\left(\frac{\Gamma_{1}\Lambda^{z}_{34}}{-i\Delta+\Gamma_{3}}+\frac{\Gamma_{1}\Lambda^{z}_{43}}{i\Delta+\Gamma_{3}}\right)\\ &\Lambda^{z}_{34}=1+\left(\frac{1}{2}-b^{2}\right)\left(\frac{\Gamma_{3}\Lambda^{z}_{11}}{\Gamma_{1}}+\frac{\Gamma_{3}\Lambda^{z}_{22}}{\Gamma_{1}}\right)\\ &\Lambda^{z}_{43}=1+\left(\frac{1}{2}-b^{2}\right)\left(\frac{\Gamma_{3}\Lambda^{z}_{11}}{\Gamma_{1}}+\frac{\Gamma_{3}\Lambda^{z}_{22}}{\Gamma_{1}}\right).\end{split}

(42)

Note that this equation treats the phonon-induced spin flips and the spin precessional dynamics self-consistently.
The first two lines in Eq. (42) give $\Lambda^{z}_{11}=\Lambda^{z}_{22}$ , while from the third and fourth lines we obtain $\Lambda^{z}_{34}=\Lambda^{z}_{43}$ . Therefore, the solution of the spin-phonon BSE is

\begin{split}&\Lambda^{z}_{11}=\Lambda^{z}_{22}=\frac{\Gamma_{3}^{2}+\Delta^{2}+\Gamma_{3}\Gamma_{1}(1-2b^{2})}{\Gamma_{3}^{2}4b^{2}(1-b^{2})+\Delta^{2}}\\ &\Lambda^{z}_{34}=\Lambda^{z}_{43}=\frac{(\Gamma_{3}^{2}+\Delta^{2})(\Gamma_{1}+\Gamma_{3}(1-2b^{2}))}{(\Gamma_{3}^{2}4b^{2}(1-b^{2})+\Delta^{2})\Gamma_{1}}.\end{split}

(43)

The state-dependent spin relaxation times $\tau^{z}_{11}$ , $\tau^{z}_{22}$ , $\tau^{z}_{34}$ , and $\tau^{z}_{43}$ are obtained using Eq. (20):

\begin{split}&\tau^{z}_{11}=\tau^{z}_{22}=\frac{\Gamma_{3}^{2}+\Delta^{2}+\Gamma_{3}\Gamma_{1}(1-2b^{2})}{(\Gamma_{3}^{2}4b^{2}(1-b^{2})+\Delta^{2})\Gamma_{1}}\\ &\tau^{z}_{34}=\tau^{z}_{43}=\frac{(\Gamma_{3}^{2}+\Delta^{2})(\Gamma_{1}+\Gamma_{3}(1-2b^{2}))}{(\Gamma_{3}^{2}4b^{2}(1-b^{2})+\Delta^{2})\Gamma_{1}\Gamma_{3}}.\end{split}

(44)

Figure 3(a) shows the vertex correction for the lowest state, $\Lambda_{11}^{z}$ in Eq. (43), and Fig. 3(b) shows the spin relaxation time ${\tau_{11}^{z}}$ from Eq. (44), both plotted as a function of the $e$ -ph collision time $\tau_{1}=1/\Gamma_{1}$ . These quantities show three distinct regimes with a qualitatively different dependence on the $e$ -ph collision time, which correspond to the EY, DP and strong-precession regimes. Our formalism encompasses these three regimes [29] because it can capture both spin-flip scattering and spin precession. The physics of these regimes is discussed below.

IV.3 Elliott-Yafet regime

We first focus on the EY regime, where spin relaxation occurs primarily through spin-flip transitions. In this regime, the spin-flip probability is small since $b\ll 1$ , and the spin precession rate is much smaller than the scattering rate, $\Delta\ll 2\Gamma_{3}b$ . Using these conditions in Eq. (43), we obtain the following vertex corrections in the EY regime:

\begin{split}\Lambda^{z}_{11}=\Lambda^{z}_{22}=\frac{\tau_{1}+\tau_{3}}{4b^{2}}\frac{1}{\tau_{1}}\\ \Lambda^{z}_{34}=\Lambda^{z}_{43}=\frac{\tau_{1}+\tau_{3}}{4b^{2}}\frac{1}{\tau_{3}},\end{split}

(45)

where $\tau_{1}=1/\Gamma_{1}$ and $\tau_{3}=1/\Gamma_{3}$ are state-dependent $e$ -ph collision times. These vertex corrections determine the spin relaxation times via Eq. (34):

\tau^{z}_{11}=\tau^{z}_{22}=\tau^{z}_{34}=\tau^{z}_{43}=\frac{\tau_{1}+\tau_{3}}{4b^{2}}.

(46)

These results, shown in Fig. 3, approximate well the solution of the spin-phonon BSE in the EY regime.
In the conventional theory of EY spin relaxation, the spin relaxation times are proportional to $1/b^{2}$ and to the $e$ -ph collision times (here, $\tau_{1}$ and $\tau_{3}$ ) [28, 44, 45]. Our results in Eq. (46) are consistent with that trend, although the spin relaxation times are proportional to the average of the $e$ -ph collision times of the two levels, $(\tau_{1}+\tau_{3})/2$ , and not to their individual values. This difference is a result of considering both forward- and back-scattering processes between the two electronic levels in the spin-phonon BSE [4], different from the simpler sRTA [28, 44, 45] which neglects electron back-scattering.

IV.4 Dyakonov-Perel regime

Next, we discuss the regime where spin precession is important and governs spin relaxation. Here, the spin-flip probability is still small ( $b\ll 1$ ) but the internal magnetic field due to SOC is significant, such that $\Delta\gg 2\Gamma_{3}b$ . Inserting these conditions in Eqs. (43)-(44), the vertex corrections become

\begin{split}&\Lambda^{z}_{11}=\Lambda^{z}_{22}=\frac{\Gamma_{3}^{2}+\Gamma_{3}\Gamma_{1}+\Delta^{2}}{\Delta^{2}},\\ &\Lambda^{z}_{34}=\Lambda^{z}_{43}=\frac{(\Gamma_{3}^{2}+\Delta^{2})(\Gamma_{1}+\Gamma_{3})}{\Delta^{2}\Gamma_{1}},\end{split}

(47)

and for the spin relaxation times we obtain

\begin{split}&\tau^{z}_{11}=\tau^{z}_{22}=\frac{\Gamma_{3}^{2}+\Gamma_{3}\Gamma_{1}+\Delta^{2}}{\Delta^{2}\Gamma_{1}},\\ &\tau^{z}_{34}=\tau^{z}_{43}=\frac{(\Gamma_{3}^{2}+\Delta^{2})(\Gamma_{1}+\Gamma_{3})}{\Delta^{2}\Gamma_{1}\Gamma_{3}}.\end{split}

(48)

These results describe spin relaxation governed by spin precession and renormalized by the $e$ -ph interactions.
In the DP regime, the $e$ -ph scattering rates are much greater than the bare spin precession rate $\Delta$ , and thus $\Gamma_{1,3}\gg\Delta$ . The vertex corrections in the DP regime become

\begin{split}&\Lambda^{z}_{11}=\Lambda^{z}_{22}=\frac{\Gamma_{3}^{2}+\Gamma_{3}\Gamma_{1}}{\Delta^{2}\Gamma_{1}}\Gamma_{1}\\ &\Lambda^{z}_{34}=\Lambda^{z}_{43}=\frac{\Gamma_{3}^{2}+\Gamma_{3}\Gamma_{1}}{\Delta^{2}\Gamma_{1}}\Gamma_{3},\end{split}

(49)

while for the spin relaxation times we obtain

\tau^{z}_{11}=\tau^{z}_{22}=\tau^{z}_{34}=\tau^{z}_{43}=\frac{\Gamma_{3}^{2}+\Gamma_{3}\Gamma_{1}}{\Delta^{2}\Gamma_{1}}.

(50)

These results, shown in Fig. 3, are an excellent approximation to the BSE solution in the DP regime.
A hallmark of DP relaxation is the inverse proportionality between the spin relaxation and $e$ -ph collision times [19], a trend captured by our treatment of the DP regime [see Eq. (50)]. Our formalism shows in Eq. (49) that this trend originates from the inverse-square scaling of the vertex corrections with the $e$ -ph collision times, $\Lambda\sim\Gamma_{1,3}^{2}\sim 1/\tau_{1,3}^{2}$ . Note that in the EY regime, rescaling the $e$ -ph collision times by a constant factor has no effect on the vertex corrections, which depend only on the ratio $\tau_{1}/\tau_{3}$ , and thus the EY spin relaxation times are proportional to the $e$ -ph collision times. The situation is different in the DP regime, where rescaling the $e$ -ph collision times changes the vertex corrections. These scaling trends can be employed in real materials to identify the dominant microscopic mechanisms for spin relaxation and decoherence [29].

IV.5 Strong-precession regime

A third, distinct regime is realized when the Zeeman splitting is much greater than the $e$ -ph scattering rates, $\Delta\gg\Gamma_{1,3}$ , such that the spins precess much faster than the rate of $e$ -ph collisions. In this strong-precession regime [43, 46, 47], the vertex corrections become

\begin{split}&\Lambda^{z}_{11}=\Lambda^{z}_{22}=1,\\ &\Lambda^{z}_{34}=\Lambda^{z}_{43}=1+\frac{\Gamma_{3}}{\Gamma_{1}}\end{split}

(51)

and can be approximated as $\Lambda\approx 1$ . Therefore, the vertex corrections are less important than in the EY or DP regimes. It follows that in the strong-precession regime the spin relaxation times are proportional to the $e$ -ph collision times:

\begin{split}&\tau^{z}_{11}=\tau^{z}_{22}=\frac{1}{\Gamma_{1}},\\ &\tau^{z}_{34}=\tau^{z}_{43}=\frac{1}{\Gamma_{1}}+\frac{1}{\Gamma_{3}}.\end{split}

(52)

These results, shown in Fig. 3, approximate well the BSE solution in the strong-precession regime. Their physical interpretation is interesting: in the strong-precession regime, the spins precess for many full cycles between phonon collisions, randomizing the spin direction. Consequently, the spin relaxation times become equal to the $e$ -ph collision times, and thus the vertex corrections $\Lambda\approx 1$ can be neglected.

IV.6 Uncovering the three regimes in GaAs

To identify the three spin relaxation mechanisms in a real material, we study spin relaxation in GaAs using our first-principles implementation of the spin-phonon BSE [29]. We focus on spin relaxation for conduction band electrons in GaAs at room temperature (300 K), and investigate how the spin relaxation times depend on the $e$ -ph collision times.
We compute the ground state and band structures of GaAs with density functional theory (DFT), using a plane-wave basis in the Quantum ESPRESSO code [48] and employing the HSE06 hybrid functional [49] to obtain an accurate band gap and electronic structure. We use fully relativistic norm-conserving pseudopotentials, generated in the local-density approximation (LDA) with Pseudo Dojo [50], together with a kinetic energy cutoff of 72 Ry and a relaxed lattice constant of 5.60 Å. The phonon energies and perturbation potentials are computed using density functional perturbation theory (DFPT) [7]. Using our perturbo code [31], we compute the $e$ -ph matrix elements on coarse Brillouin zone grids with $8\times 8\times 8$ $\bm{k}$ and $\bm{q}$ points, following which we generate spinor Wannier functions with the Wannier90 code [51] and use them in perturbo [31], with a method we developed in Ref. [28], to jointly interpolate the $e$ -ph matrix elements and spin matrices. The long-range quadrupole $e$ -ph interactions [10, 52, 11, 53] are included to fully account for $e$ -ph coupling with long-wavelength phonons. We interpolate the $e$ -ph matrix elements to fine Brillouin zone grids with up to $200\times 200\times 200$ $\bm{k}$ and $\bm{q}$ points, and a 10 meV Gaussian broadening for the energy-conserving delta functions in the $e$ -ph scattering rates [41]. The spin relaxation times in Eq. (30) are computed using the tetrahedron integration method [54], assuming a non-degenerate electron concentration of $10^{16}~{}\text{cm}^{-3}$ . The spin-phonon BSE in Eq. (32) is solved with an augmented iterative approach described in Ref. [29].

Using our first-principles spin-phonon BSE, we compute the spin relaxation time for electron spins in GaAs, and obtain a value 51 ps at 300 K in excellent agreement with the experimental value of 42 ps [55]. We also obtain an average $e$ -ph collision time of 410 fs, computing using

\langle\tau^{\text{e-ph}}\rangle=\frac{\sum_{n\bm{k}}\tau^{\text{e-ph}}_{nn\bm{k}}(\varepsilon_{n\bm{k}})(-\frac{df_{n\bm{k}}}{d\varepsilon})}{\sum_{n\bm{k}}\,\big{(}\!-\frac{df_{n\bm{k}}}{d\varepsilon}\big{)}}

(53)

and an average vertex correction $\langle\Lambda^{z}\rangle\!\approx\!171$ , computed from

\langle\Lambda^{z}\rangle=\frac{\sum_{n\bm{k}}\absolutevalue{s^{z}_{nn\bm{k}}}^{2}\Lambda^{z}_{nn\bm{k}}(\varepsilon_{n\bm{k}})(-\frac{df_{n\bm{k}}}{d\varepsilon})}{\sum_{n\bm{k}}\absolutevalue{s^{z}_{nn\bm{k}}}^{2}\,\big{(}\!-\frac{df_{n\bm{k}}}{d\varepsilon}\big{)}}.

(54)

These “real-material” values for GaAs are shown with a dot in Fig. 4(a) and (b), where we also show spin relaxation times and vertex corrections obtained by artificially varying the average $e$ -ph collision time (by rescaling the $e$ -ph matrix elements). The resulting trends show that the spin relaxation times are inversely proportional to the $e$ -ph collision time, placing GaAs in the DP spin-relaxation regime.
Upon decreasing the average $e$ -ph collision time $\langle\tau^{\text{e-ph}}\rangle$ in Fig. 4(a), the system evolves from the DP to the EY regime, and the vertex correction first increases and then saturates to a maximal value of $\sim$ 1900 in the short $e$ -ph collision time limit. As a consequence, the spin relaxation time in Fig. 4(b) peaks at $\langle\tau^{\text{e-ph}}\rangle\!\approx\!0.09$ ps and then decreases linearly at shorter $e$ -ph collision times. This crossover from the DP to the EY mechanism, obtained here by artificially tuning the $e$ -ph collision time in GaAs, is consistent with the results from our model system.
Conversely, when the $e$ -ph collision time is increased, the system transitions from the DP to the strong-precession regime: Fig. 4(a) shows that the vertex correction first decreases as $\langle\tau^{\text{e-ph}}\rangle^{-2}$ and then plateaus to a minimal value close to unity for long $e$ -ph collision times. Accordingly, the spin relaxation time in Fig. 4(b) reaches a minimum for $\langle\tau^{\text{e-ph}}\rangle\approx 3$ ps and then increases linearly at longer $e$ -ph collision times.
In summary, real materials exhibit the same spin relaxation mechanisms and vertex correction trends as our two-level system. These regimes for phonon-induced spin dynamics are summarized in Table 2. The presence of three distinct spin-relaxation regimes is general, and we expect it to be valid beyond the case of a simple semiconductor (GaAs) studied here.

Table 2: Summary of the Elliott-Yafet, Dyakonov-Perel, and strong-precession regimes.

	EY	DP	Strong-precession
	Spin-flip	Spin precession	Spin precession
$\Lambda$	$\sim{1}/{b^{2}}$	$\sim{1}/{(\tau^{\text{e-ph}})^{2}}$	$\sim 1$
$\tau^{s}$	$\sim{\tau^{\text{e-ph}}}/{b^{2}}$	$\sim{1}/{\tau^{\text{e-ph}}}$	$\sim\tau^{\text{e-ph}}$

V Discussion

Our formalism can capture all three of the EY, DP and strong-precession regimes in a unified framework. The reason can be inferred from the diagrammatic representation of the spin-phonon BSE in Fig. 1(c). In the diagrams, the EY relaxation is due to $e$ -ph scattering in the presence of spin mixing, which originates from the $e$ -ph interactions $\left[g_{n^{\prime}m^{\prime}\nu}(\bm{k},\bm{q})\right]^{*}$ and $g_{nm\nu}(\bm{k},\bm{q})$ and the wiggly line in the kernel of the BSE. The DP and strong-precession mechanisms are included by virtue of the electron propagators with two different band indices, $\mathcal{G}_{ml\bm{k}+\bm{q}}\mathcal{G}_{l^{\prime}m^{\prime}\bm{k}+\bm{q}}$ , placed between the wiggly line and the spin vertex. These propagators take into account spin precession (due to the SOC field) between $e$ -ph collisions. This elegant formalism captures a wide range of spin physics in a single diagram.
Although our discussion has focused on $T_{1}$ spin relaxation times, the spin decoherence times $T_{2}$ at finite magnetic fields can also be computed, as we plan to show in future work. Finally, the approach presented in this work for the phonon-dressed vertex is general and goes beyond spin relaxation. It can be employed to study the dynamics of any observable that couples with phonons, for which we also expect to find the three phonon-induced relaxation regimes discussed above for spin dynamics.

VI Conclusion

We have formulated a theory for the vertex corrections from $e$ -ph interactions to the susceptibility of a non-diagonal operator. The key result is a self-consistent BSE to calculate the phonon-dressed vertex, which encodes the dynamics of the operator coupling with phonons. When applied to spin, this approach enables quantitative calculations of spin relaxation and decoherence [29]. We have shown that our spin-phonon BSE captures both spin-flip transitions and spin precession, unifying the treatment of three spin decoherence mechanisms (EY, DP, and strong-precession) conventionally treated with separate heuristic models. By leveraging efficient workflows for first-principles $e$ -ph calculations [31], our method enables quantitative studies of spin relaxation and decoherence in a wide range of bulk and two-dimensional materials, as we show in the companion paper [29]. These advances open new avenues for understanding spin relaxation and decoherence in spintronics, magnetism, multiferroics, quantum materials and quantum technologies.

Acknowledgements.

This work was supported by the National Science Foundation under Grants No. DMR-1750613 and QII-TAQS 1936350, which provided for method development, and Grant No. OAC-2209262, which provided for code development. J.P. acknowledges support by the Korea Foundation for Advanced Studies.

Appendix A Ward identity

We derive a Ward identity for our BSE in Eq. (13). The Ward identity relates the vertex corrections with the electron self-energy [4, 1, 36], and guarantees that diagrams are taken into account consistently in the self-energy and in the BSE for the vertex [4, 1, 36].
For a system with Hamiltonian $H$ , the operator $\hat{A}$ is related to the the negative derivative of the Hamiltonian with respect to the external field $\mathcal{F}$ , $\hat{A}=-\gradient_{\mathcal{F}}H$ . We compute the change of the Fan-Migdal self-energy in Eq. (2) with respect to $\mathcal{F}$ , by taking the derivative

\begin{split}&{\gradient_{\mathcal{F}}\Sigma_{nn^{\prime}\bm{k}}(i\omega_{a})}\\ &~{}~{}=-\frac{1}{\beta N_{q}V_{\text{uc}}}\sum_{mm^{\prime}ll^{\prime}\bm{q}\nu,iq_{c}}\left[g_{n^{\prime}m^{\prime}\nu}(\bm{k},\bm{q})\right]^{*}g_{nm\nu}(\bm{k},\bm{q})\\ &~{}~{}~{}\times\mathcal{D}_{\nu\bm{q}}(iq_{c})\left(\gradient_{\mathcal{F}}\mathcal{G}_{mm^{\prime}\bm{k}}(i\omega_{a}+iq_{c})\right).\end{split}

(55)

Employing the matrix identity

\begin{split}{\gradient_{\mathcal{F}}\mathcal{G}_{mm^{\prime}\bm{k}}}&=-\sum_{ll^{\prime}}\mathcal{G}_{ml\bm{k}}\left[{\gradient_{\mathcal{F}}(\mathcal{G}^{-1}})\right]_{ll^{\prime}\bm{k}}\mathcal{G}_{l^{\prime}m^{\prime}\bm{k}}\\ &=-\sum_{ll^{\prime}}\mathcal{G}_{ml\bm{k}}\left(-\gradient_{\mathcal{F}}H_{ll^{\prime}\bm{k}}-\gradient_{\mathcal{F}}\Sigma_{ll^{\prime}\bm{k}}\right)\mathcal{G}_{l^{\prime}m^{\prime}\bm{k}},\end{split}

(56)

we obtain a self-consistent equation for the self-energy derivatives:

\begin{split}&{\gradient_{\mathcal{F}}\Sigma_{nn^{\prime}\bm{k}}(i\omega_{a})}\\ &~{}~{}=-\frac{1}{\beta N_{q}V_{\text{uc}}}\sum_{mm^{\prime}ll^{\prime}\bm{q}\nu,iq_{c}}\left[g_{n^{\prime}m^{\prime}\nu}(\bm{k},\bm{q})\right]^{*}g_{nm\nu}(\bm{k},\bm{q})\\ &~{}~{}~{}\times\mathcal{D}_{\nu\bm{q}}(iq_{c})\mathcal{G}_{ml\bm{k}+\bm{q}}(i\omega_{a}+iq_{c})\mathcal{G}_{l^{\prime}m^{\prime}\bm{k}+\bm{q}}(i\omega_{a}+iq_{c})\\ &~{}~{}~{}\times\left(\gradient_{\mathcal{F}}H_{ll^{\prime}\bm{k}}+\gradient_{\mathcal{F}}\Sigma_{ll^{\prime}\bm{k}}(i\omega_{a}+iq_{c})\right).\end{split}

(57)

Comparing Eq. (57) and Eq. (8) in the $i\nu_{b}\to 0$ limit, we discover the Ward identity expressed in terms of the Hamiltonian, vertex corrections, and self-energy:

(\gradient_{\mathcal{F}}H)\bm{\Lambda}(i\omega_{a},i\omega_{a})=\gradient_{\mathcal{F}}H+{\gradient_{\mathcal{F}}\Sigma(i\omega_{a})},

(58)

where $(\gradient_{\mathcal{F}}H)\bm{\Lambda}(i\omega_{a},i\omega_{a})=(\frac{\partial H_{nn^{\prime}\bm{k}}}{{\mathcal{\partial F^{\alpha}}}}){\Lambda}^{\alpha}_{nn^{\prime}\bm{k}}(i\omega_{a},i\omega_{a})$ . Equation (58) can be equivalently expressed in terms of the operator matrix elements $A^{\alpha}_{nn^{\prime}\bm{k}}$ ,

A^{\alpha}_{nn^{\prime}\bm{k}}\Lambda^{\alpha}_{nn^{\prime}\bm{k}}(i\omega_{a},i\omega_{a})=A^{\alpha}_{nn^{\prime}\bm{k}}-\frac{\partial\Sigma_{nn^{\prime}\bm{k}}(i\omega_{a})}{\partial\mathcal{F}^{\alpha}}.

(59)

Our expression for the Ward identity in Eqs. (58)-(59) is consistent with the Ward identity for the velocity operator derived in Ref. [4]:

v^{\alpha}_{nn\bm{k}}\Lambda^{\alpha}_{nn\bm{k}}(i\omega_{a},i\omega_{a})=v^{\alpha}_{nn\bm{k}}+\frac{\partial\Sigma_{nn\bm{k}}(i\omega_{a})}{\partial k^{\alpha}},

(60)

as the velocity operator is defined as $\hat{v}=\gradient_{\mathcal{F}}H$ with $\mathcal{F}=\bm{k}$ . This result further validates our BSE for the phonon-dressed vertex.

References

Mahan [2000] G. D. Mahan, Many-Particle Physics, 3rd ed., Physics of Solids and Liquids (2000).
Rohlfing and Louie [2000] M. Rohlfing and S. G. Louie, Phys. Rev. B 62, 4927 (2000).
Sekine et al. [2017] A. Sekine, D. Culcer, and A. H. MacDonald, Phys. Rev. B 96, 235134 (2017).
Kim et al. [2019] S. Kim, S. Woo, and H. Min, Phys. Rev. B 99, 165107 (2019).
Zhou and Bernardi [2019] J.-J. Zhou and M. Bernardi, Phys. Rev. Res. 1, 033138 (2019).
Desai et al. [2021] D. C. Desai, B. Zviazhynski, J.-J. Zhou, and M. Bernardi, Phys. Rev. B 103, L161103 (2021).
Baroni et al. [2001] S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Rev. Mod. Phys. 73, 515 (2001).
Giustino et al. [2007] F. Giustino, M. L. Cohen, and S. G. Louie, Phys. Rev. B 76, 165108 (2007).
Agapito and Bernardi [2018] L. A. Agapito and M. Bernardi, Phys. Rev. B 97, 235146 (2018).
Jhalani et al. [2020] V. A. Jhalani, J.-J. Zhou, J. Park, C. E. Dreyer, and M. Bernardi, Phys. Rev. Lett. 125, 136602 (2020).
Brunin et al. [2020a] G. Brunin, H. P. C. Miranda, M. Giantomassi, M. Royo, M. Stengel, M. J. Verstraete, X. Gonze, G.-M. Rignanese, and G. Hautier, Phys. Rev. Lett. 125, 136601 (2020a).
Lu et al. [2020] I.-T. Lu, J. Park, J.-J. Zhou, and M. Bernardi, npj Comput. Mater. 6, 1 (2020).
Lu et al. [2021] I.-T. Lu, J.-J. Zhou, J. Park, and M. Bernardi, Phys. Rev. Mater. 6, L010801 (2022).
Zhou et al. [2021a] J.-J. Zhou, J. Park, I. Timrov, A. Floris, M. Cococcioni, N. Marzari, and M. Bernardi, Phys. Rev. Lett. 127, 126404 (2021a).
Bernardi [2016] M. Bernardi, Eur. Phys. J. B 89, 239 (2016).
Pickett [1989] W. E. Pickett, Rev. Mod. Phys. 61, 433 (1989).
Balasubramanian et al. [2008] G. Balasubramanian, I. Y. Chan, R. Kolesov, M. Al-Hmoud, J. Tisler, C. Shin, C. Kim, A. Wojcik, P. R. Hemmer, A. Krueger, T. Hanke, A. Leitenstorfer, R. Bratschitsch, F. Jelezko, and J. Wrachtrup, Nature 455, 648 (2008).
Jansen [2012] R. Jansen, Nat. Mater. 11, 400 (2012).
Žutić et al. [2004] I. Žutić, J. Fabian, and S. Das Sarma, Rev. Mod. Phys. 76, 323 (2004).
Veldhorst et al. [2015] M. Veldhorst, C. H. Yang, J. C. C. Hwang, W. Huang, J. P. Dehollain, J. T. Muhonen, S. Simmons, A. Laucht, F. E. Hudson, K. M. Itoh, A. Morello, and A. S. Dzurak, Nature 526, 410 (2015).
Noiri et al. [2022] A. Noiri, K. Takeda, T. Nakajima, T. Kobayashi, A. Sammak, G. Scappucci, and S. Tarucha, Nature 601, 338 (2022).
Petit et al. [2020] L. Petit, H. G. J. Eenink, M. Russ, W. I. L. Lawrie, N. W. Hendrickx, S. G. J. Philips, J. S. Clarke, L. M. K. Vandersypen, and M. Veldhorst, Nature 580, 355 (2020).
Ziman [2001] J. M. Ziman, Electrons and Phonons: The Theory of Transport Phenomena in Solids (Oxford University Press, New York, 2001).
Lundstrom [2009] M. Lundstrom, Fundamentals of Carrier Transport (Cambridge University Press, New York, 2009).
Elliott [1954] R. J. Elliott, Phys. Rev. 96, 266 (1954).
Yafet [1963] Y. Yafet, in Solid State Physics, Vol. 14, edited by F. Seitz and D. Turnbull (1963) pp. 1–98.
Dyakonov and Perel [1972] M. Dyakonov and V. Perel, Sov. Phys. Solid State USSR 13, 3023 (1972).
Park et al. [2020a] J. Park, J.-J. Zhou, and M. Bernardi, Phys. Rev. B 101, 045202 (2020a).
Park et al. [2022] J. Park, J.-J. Zhou, and M. Bernardi arXiv:2203.06401 (2022).
Cardona [2001] M. Cardona, Phys. Status Solidi A 188, 1209 (2001).
Zhou et al. [2021b] J.-J. Zhou, J. Park, I.-T. Lu, I. Maliyov, X. Tong, and M. Bernardi, Comput. Phys. Commun. 264, 107970 (2021b).
Stefanucci and van Leeuwen [2013] G. Stefanucci and R. van Leeuwen, Nonequilibrium Many-Body Theory of Quantum Systems: A Modern Introduction (Cambridge University Press, 2013).
Lihm and Park [2020] J.-M. Lihm and C.-H. Park, Phys. Rev. B 101, 121102(R) (2020).
Allen and Heine [1976] P. B. Allen and V. Heine, J. Phys. C: Solid State Phys. 9, 2305 (1976).
Vollhardt and Wölfle [1980] D. Vollhardt and P. Wölfle, Phys. Rev. B 22, 4666 (1980).
Ward [1950] J. C. Ward, Phys. Rev. 78, 182 (1950).
Shen et al. [2014] K. Shen, G. Vignale, and R. Raimondi, Phys. Rev. Lett. 112, 096601 (2014).
Offidani et al. [2018] M. Offidani, R. Raimondi, and A. Ferreira, Condens. Matter 3, 18 (2018).
Coleman [2015] P. Coleman, Introduction to Many-Body Physics (Cambridge University Press, 2015).
Li [2015] W. Li, Phys. Rev. B 92, 075405 (2015).
Zhou and Bernardi [2016] J.-J. Zhou and M. Bernardi, Phys. Rev. B 94, 201201(R) (2016).
Culcer et al. [2017] D. Culcer, A. Sekine, and A. H. MacDonald, Phys. Rev. B 96, 035106 (2017).
Burkov and Balents [2004] A. A. Burkov and L. Balents, Phys. Rev. B 69, 245312 (2004).
Baral et al. [2016] A. Baral, S. Vollmar, S. Kaltenborn, and H. C. Schneider, New J. Phys. 18, 023012 (2016).
Vollmar et al. [2017] S. Vollmar, D. J. Hilton, and H. C. Schneider, Phys. Rev. B 96, 075203 (2017).
Szolnoki et al. [2017] L. Szolnoki, A. Kiss, B. Dóra, and F. Simon, Sci. Rep. 7, 9949 (2017).
Csősz et al. [2020] G. Csősz, L. Szolnoki, A. Kiss, B. Dóra, and F. Simon, Phys. Rev. Research 2, 033058 (2020).
Giannozzi et al. [2009] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, J. Phys.: Condens. Matter 21, 395502 (2009).
Heyd et al. [2003] J. Heyd, G. E. Scuseria, and M. Ernzerhof, J. Chem. Phys. 118, 8207 (2003).
van Setten et al. [2018] M. J. van Setten, M. Giantomassi, E. Bousquet, M. J. Verstraete, D. R. Hamann, X. Gonze, and G. M. Rignanese, Comput. Phys. Commun. 226, 39 (2018).
Pizzi et al. [2020] G. Pizzi, V. Vitale, R. Arita, S. Blügel, F. Freimuth, G. Géranton, M. Gibertini, D. Gresch, C. Johnson, T. Koretsune, J. Ibañez-Azpiroz, H. Lee, J.-M. Lihm, D. Marchand, A. Marrazzo, Y. Mokrousov, J. I. Mustafa, Y. Nohara, Y. Nomura, L. Paulatto, S. Poncé, T. Ponweiser, J. Qiao, F. Thöle, S. S. Tsirkin, M. Wierzbowska, N. Marzari, D. Vanderbilt, I. Souza, A. A. Mostofi, and J. R. Yates, J. Phys. Condens. Matter 32, 165902 (2020).
Park et al. [2020b] J. Park, J.-J. Zhou, V. A. Jhalani, C. E. Dreyer, and M. Bernardi, Phys. Rev. B 102, 125203 (2020b).
Brunin et al. [2020b] G. Brunin, H. P. C. Miranda, M. Giantomassi, M. Royo, M. Stengel, M. J. Verstraete, X. Gonze, G.-M. Rignanese, and G. Hautier, Phys. Rev. B 102, 094308 (2020b).
Blöchl et al. [1994] P. E. Blöchl, O. Jepsen, and O. K. Andersen, Phys. Rev. B 49, 16223 (1994).
Oertel et al. [2008] S. Oertel, J. Hübner, and M. Oestreich, Appl. Phys. Lett. 93, 132112 (2008).