Re-examination of the Elliott–Yafet spin-relaxation mechanism

We start our derivation by writing down the interaction potential of a single electron in a lattice in terms of lattice site coordinates ${{\bf{R}}}_{n}$ as ${\hat{v}}_{{\rm{e}}-{\rm{L}}}({\bf{x}},\{{{\bf{R}}}_{n}\})={\displaystyle \sum }_{n}{\hat{v}}_{{\rm{e}}-\mathrm{ion}}({\bf{x}}-{{\bf{R}}}_{n})$, where

$$\begin{eqnarray}&&{\hat{v}}_{{\rm{e}}-\mathrm{ion}}({\bf{x}})={v}_{\mathrm{eff}}({\bf{x}})+\xi [{\rm{\nabla }}{v}_{\mathrm{eff}}({\bf{x}})\times \hat{{\bf{p}}}]\cdot \hat{{\bf{s}}}.\end{eqnarray} \tag{ 8 }$$

Here, ${v}_{\mathrm{eff}}$ is the effective electrostatic Coulomb interaction between the electron and the ionic core and $\xi ={\hslash }/(4\;{m}^{2}{c}^{2})$. We have made here the rigid-ion approximation, that is, we assume that the core electrons and the nuclei form ions, which set up a (pseudo-)potential for the valence electrons ${\hat{v}}_{{\rm{e}}-\mathrm{ion}}({\bf{x}})$, which follows rigidly the motion of the ions. This is, in general, a good approximation for semiconductors [29] and simple metals [30]. We stress that this approximation places no fundamental restriction on the following development, but the equations of motion become more complicated without it. We denote single-particle operators acting on the electronic space and spin variables by small letters and a hat; for instance, $\hat{{\bf{p}}}=-{\rm{i}}{\hslash }{\rm{\nabla }}$. The second term in (8) is the spin–orbit interaction. We remind the reader that, due to lack of rotational symmetry, this potential does not commute with the electronic angular momentum operator $\hat{{\bf{x}}}\times \hat{{\bf{p}}}$ of the electrons, and therefore does not conserve the electronic orbital angular momentum.

For ferromagnets, the electron–ion interaction is sometimes written as

$$\begin{eqnarray}&&{\hat{v}}_{{\rm{e}}-\mathrm{ion}}({\bf{x}})=\hat{v}({\bf{x}})+\xi [{\rm{\nabla }}\hat{v}({\bf{x}})\times \hat{{\bf{p}}}]\cdot \hat{{\bf{s}}},\end{eqnarray} \tag{ 9 }$$

with

$$\begin{eqnarray}&&\hat{v}={v}_{\mathrm{eff}}({\bf{x}})\hat{1}+\delta v({\bf{x}}){\hat{\sigma }}_{z},\end{eqnarray} \tag{ 10 }$$

where ${v}_{\mathrm{eff}}=({v}_{\uparrow }+{v}_{\downarrow })/2$, ${\rm{\Delta }}v=({v}_{\uparrow }-{v}_{\downarrow })/2$, ${\hat{\sigma }}_{z}$ is a Pauli matrix, and $\hat{1}$ the identity in spin space. In both cases, the electron–ion interaction has a spin–diagonal and an explicitly spin-dependent contribution. We will use the electron–ion interaction hamiltonian in the form (8) as it applies to semiconductors, but distinguish in the following only between spin–diagonal and explicitly spin-dependent contributions, which also applies to (9). In semiconductors, these two contributions could be called, respectively, Elliott and Overhauser (or Yafet) contributions. However if it is not clearly stated to which one of the interaction hamiltonians (8) and (9) one refers, this terminology can be ambiguous.

We follow the usual treatment of the electron–phonon coupling [30] by considering small deviations ${{\bf{R}}}_{n}={{\bf{R}}}_{n}^{(0)}+{{\bf{Q}}}_{n}$ around the equilibrium configuration ${{\bf{R}}}_{n}^{(0)}$ and expand ${v}_{\mathrm{eff}}({\bf{r}}-{{\bf{R}}}_{n})$ with respect to ${{\bf{R}}}_{n}$. The resulting electron–phonon hamiltonian in second quantization is

$$\begin{eqnarray}&&\left.{H}_{{\rm{e}}-\mathrm{pn}}=\displaystyle \sum _{n}\int {{\rm{\Psi }}}^{\dagger }({\bf{r}}){\left[{{\bf{Q}}}_{n}\cdot \displaystyle \frac{\partial {\hat{v}}_{{\rm{e}}-\mathrm{ion}}({\bf{r}}-{{\bf{R}}}_{n})}{\partial {{\bf{R}}}_{n}}\right|}_{{{\bf{R}}}_{n}^{(0)}}\right]{\rm{\Psi }}({\bf{r}})\ {{\rm{d}}}^{3}r,\end{eqnarray} \tag{ 11 }$$

where ${\rm{\Psi }}\equiv {({{\rm{\Psi }}}_{\uparrow },{{\rm{\Psi }}}_{\downarrow })}^{T}$ denotes a spinor field operator. Note that $\partial {\hat{v}}_{{\rm{e}}-\mathrm{ion}}/\partial {{\bf{R}}}_{n}$ is an operator in spin space.

While the properties of the electron–phonon interaction can be discussed without specifying a model for the single-particle band structure, we also need the electronic single-particle contribution to the hamiltonian for the derivation of the equations of motion. We thus consider the model hamiltonian $H={H}_{{\rm{e}}}+{H}_{{\rm{e}}-\mathrm{pn}}$ with

$$\begin{eqnarray}&&{H}_{{\rm{e}}}={\displaystyle \int }_{{ \mathcal V }}{{\rm{\Psi }}}^{\dagger }({\bf{r}})\left[\displaystyle \frac{{\hat{{\bf{p}}}}^{2}}{2{m}_{0}}+\displaystyle \sum _{n}{\hat{v}}_{{\rm{e}}-\mathrm{ion}}({\bf{r}}-{{\bf{R}}}_{n}^{(0)})\right]{\rm{\Psi }}({\bf{r}}){{\rm{d}}}^{3}r.\end{eqnarray} \tag{ 12 }$$

We follow the standard approach [31] and expand the field operators according to

$$\begin{eqnarray}&&{\rm{\Psi }}({\bf{r}})=\displaystyle \frac{1}{\sqrt{{ \mathcal V }}}\displaystyle \sum _{\mu }\displaystyle \sum _{{\bf{k}}}{\varphi }_{\mu {\bf{k}}}({\bf{r}}){c}_{\mu {\bf{k}}},\end{eqnarray} \tag{ 13 }$$

where ${\varphi }_{\mu {\bf{k}}}({\bf{r}})={{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {\bf{r}}}{u}_{\mu {\bf{k}}}({\bf{r}})$ are Bloch spinor wavefunctions with band and crystal momentum labels $(\mu {\bf{k}})$. The destruction operator corresponding to ${\varphi }_{\mu {\bf{k}}}$ is ${c}_{\mu {\bf{k}}}$, and ${ \mathcal V }$ is the normalization volume of the crystal. A Greek index labels the bands including spin–orbit coupling. In the case of Kramers degenerate bands, to which we specialize below, it includes a band index b and an additional label λ to distinguish the degenerate states, which is sometimes called a pseudospin index, i.e., $\mu =(b,\lambda )$.

Instead of the Bloch wavefunctions, which are orthogonal on the whole crystal volume ${ \mathcal V }$, we will mainly deal with the lattice periodic ${u}_{\mu {\bf{k}}}$'s, which are defined and orthogonal on the unit cell Ω, because we intend to derive matrix elements on the unit cell Ω. Matrix elements of any single-electron operator $\hat{a}$ that occur in the present paper have the meaning

$$\begin{eqnarray}&&{\langle {u}_{\mu {\bf{k}}}| \hat{a}{u}_{{\mu }^{\prime }{{\bf{k}}}^{\prime }}\rangle }_{{\rm{\Omega }}}:= \displaystyle \sum _{\sigma }\displaystyle \frac{1}{{\rm{\Omega }}}{\displaystyle \int }_{{\rm{\Omega }}}{[{u}_{\mu {\bf{k}}}^{*}({\bf{r}})]}_{\sigma }{[\hat{a}{u}_{{\mu }^{\prime }{{\bf{k}}}^{\prime }}({\bf{r}})]}_{\sigma }{{\rm{d}}}^{3}r.\end{eqnarray} \tag{ 14 }$$

By virtue of equation (8), the electron–phonon interaction hamiltonian (11) can be split into two parts ${H}_{{\rm{e}}-\mathrm{pn}}={H}_{{\rm{e}}-\mathrm{pn}}^{(1)}+{H}_{{\rm{e}}-\mathrm{pn}}^{(2)}$, one of which

$$\begin{eqnarray}&&{H}_{{\rm{e}}-\mathrm{pn}}^{(1)}=-\displaystyle \sum _{n}\displaystyle \sum _{{\bf{k}}\mu ,{{\bf{k}}}^{\prime }{\mu }^{\prime }}{c}_{\mu {\bf{k}}}^{\dagger }{c}_{{\mu }^{\prime }{{\bf{k}}}^{\prime }}\int \displaystyle \frac{{{\rm{d}}}^{3}r}{{ \mathcal V }}{\varphi }_{\mu {\bf{k}}}^{*}({\bf{r}}){{\bf{Q}}}_{n}\cdot {\boldsymbol{\nabla }}{v}_{\mathrm{eff}}({\bf{r}}-{{\bf{R}}}_{n}^{(0)}){\varphi }_{{\mu }^{\prime }{{\bf{k}}}^{\prime }}({\bf{r}})\end{eqnarray} \tag{ 15 }$$

is spin–diagonal, whereas the other one explicitly depends on the electronic spin operator $\hat{{\bf{s}}}$

$$\begin{eqnarray}{H}_{{\rm{e}}-\mathrm{pn}}^{(2)} & = & -\xi \displaystyle \sum _{n}\displaystyle \sum _{{\bf{k}},{{\bf{k}}}^{\prime }}\displaystyle \sum _{\mu {\mu }^{\prime }}{c}_{\mu {\bf{k}}}^{\dagger }{c}_{{\mu }^{\prime }{{\bf{k}}}^{\prime }}\\ & & \displaystyle \int \displaystyle \frac{{{\rm{d}}}^{3}r}{{ \mathcal V }}{\varphi }_{\mu {\bf{k}}}^{*}({\bf{r}})[{\boldsymbol{\nabla }}({{\bf{Q}}}_{n}\cdot {\boldsymbol{\nabla }}{v}_{\mathrm{eff}}({\bf{r}}-{{\bf{R}}}_{n}^{(0)}))\times \hat{{\bf{p}}}]\cdot \hat{{\bf{s}}}\ {\varphi }_{{\mu }^{\prime }{{\bf{k}}}^{\prime }}({\bf{r}}).\end{eqnarray} \tag{ 16 }$$

We next expand the phonon displacement operators according to [31]

$$\begin{eqnarray}&&{{\bf{Q}}}_{n}(t)={\rm{i}}\displaystyle \sum _{{\bf{q}},\lambda }{{\ell }}_{{\bf{q}},\lambda }^{(0)}{B}_{{\bf{q}},\lambda }{{\rm{e}}}^{{\rm{i}}{\bf{q}}\cdot {{\bf{R}}}_{n}^{(0)}}\ {{\boldsymbol{\varepsilon }}}_{{\bf{q}},\lambda },\end{eqnarray} \tag{ 17 }$$

with ${B}_{{\bf{q}},\lambda }\equiv {b}_{-{\bf{q}},\lambda }^{\dagger }+{b}_{{\bf{q}},\lambda }$, where ${b}_{{\bf{q}},\lambda }^{(\dagger )}$ annihilates (creates) phonons with momentum ${\bf{q}}$ in the first Brillouin zone and mode index λ. We use the abbreviation of ${{\ell }}_{{\bf{q}},\lambda }^{(0)}=\sqrt{{\hslash }/(2\varrho { \mathcal V }{\omega }_{{\bf{q}},\lambda })}$, which has the unit of length, where ϱ is the density. The phonon dispersion is denoted by ${\omega }_{{\bf{q}},\lambda }$. Further, ${{\boldsymbol{\varepsilon }}}_{{\bf{q}},\lambda }$ is the polarization vector of the phonon mode with the property ${({{\boldsymbol{\varepsilon }}}_{{\bf{q}},\lambda })}^{*}=-{{\boldsymbol{\varepsilon }}}_{-{\bf{q}},\lambda }$. We introduce the Fourier transformation of the interaction potential in the form [31]

$$\begin{eqnarray}&&{v}_{\mathrm{eff}}({\bf{r}})=\displaystyle \frac{1}{N}\displaystyle \sum _{{\bf{q}}\in 1\mathrm{BZ}}\displaystyle \sum _{{\bf{G}}}{v}_{\mathrm{eff}}({\bf{q}}+{\bf{G}}){{\rm{e}}}^{{\rm{i}}({\bf{q}}+{\bf{G}})\cdot {\bf{r}}},\end{eqnarray} \tag{ 18 }$$

which has a long-range part (${\bf{G}}=0$) and short-range part (${\bf{G}}\ne 0$).

Now we split the integral over the crystal volume ${ \mathcal V }$ according to

$$\begin{eqnarray}&&\displaystyle \frac{1}{{ \mathcal V }}{\displaystyle \int }_{{ \mathcal V }}f({\bf{r}}){{\rm{d}}}^{3}r=\displaystyle \frac{1}{N{\rm{\Omega }}}\displaystyle \sum _{n}{\displaystyle \int }_{{\rm{\Omega }}}f({\bf{r}}+{{\bf{R}}}_{n}^{(0)}){{\rm{d}}}^{3}r,\end{eqnarray} \tag{ 19 }$$

use the periodicity ${u}_{\mu {\bf{k}}}({\bf{r}}+{{\bf{R}}}_{n}^{(0)})={u}_{\mu {\bf{k}}}({\bf{r}})$ and the relation ${\displaystyle \sum }_{n}{{\rm{e}}}^{{\rm{i}}{{\bf{R}}}_{n}^{(0)}\cdot ({\bf{q}}+{{\bf{k}}}^{\prime }-{\bf{k}})}=N{\delta }_{{\bf{q}},{\bf{k}}-{{\bf{k}}}^{\prime }}$. We thus obtain for the e–pn coupling hamiltonian (11)

$$\begin{eqnarray}&&{H}_{{\rm{e}}-\mathrm{pn}}=\displaystyle \sum _{{\bf{kq}}}\displaystyle \sum _{\mu {\mu }^{\prime }}\displaystyle \sum _{\lambda }{g}_{\mu {\bf{k}}+{\bf{q}},{\mu }^{\prime }{\bf{k}}}^{(\lambda )}{B}_{{\bf{q}},\lambda }{c}_{\mu {\bf{k}}+{\bf{q}}}^{\dagger }{c}_{{\mu }^{\prime }{\bf{k}}},\end{eqnarray} \tag{ 20 }$$

with the matrix element

$$\begin{eqnarray}&&{g}_{\mu {\bf{k}}+{\bf{q}},{\mu }^{\prime }{\bf{k}}}^{(\lambda )}={\langle {u}_{\mu {\bf{k}}+{\bf{q}}}| {\hat{v}}_{{\bf{k}}+{\bf{q}},{\bf{k}}}^{(\lambda )}{u}_{{\mu }^{\prime }{\bf{k}}}\rangle }_{{\rm{\Omega }}}\end{eqnarray} \tag{ 21 }$$

in terms of the e–pn interaction operator

$$\begin{eqnarray}{\hat{v}}_{{\bf{k}}+{\bf{q}},{\bf{k}}}^{(\lambda )} & = & {{\ell }}_{{\bf{q}},\lambda }^{(0)}\displaystyle \sum _{{\bf{G}}}{{\rm{e}}}^{{\rm{i}}{\bf{G}}\cdot {\bf{r}}}{v}_{\mathrm{eff}}({\bf{q}}+{\bf{G}})\\ & & [({\bf{q}}+{\bf{G}})\cdot {{\boldsymbol{\varepsilon }}}_{{\bf{q}},\lambda }]\{\underset{(1)}{\underbrace{1}}+\underset{(2)}{\underbrace{{\rm{i}}\xi [({\bf{q}}+{\bf{G}})\times (\hat{{\bf{p}}}+{\hslash }{\bf{k}})]\cdot \hat{{\bf{s}}}}}\}.\end{eqnarray} \tag{ 22 }$$

The first term in the curly braces in the above expression comes from the spin diagonal term (15) and the rest from the explicitly spin-dependent term (16). The ${\bf{q}}$ and ${\bf{k}}$ vectors are restricted to the first Brillouin zone, and we have already neglected Umklapp processes, which contribute if ${\bf{k}}+{\bf{q}}$ lies outside the first BZ. This is a common approximation for optically excited electron dynamics when one is not interested in transport processes. The long-range part of this expression results from the ${\bf{G}}={\bf{0}}$ contribution to the sum, whereas the sum over the ${\bf{G}}\ne {\bf{0}}$ defines the short-range part of the matrix elements. Grimaldi and Fulde [23] demonstrate that the long-range part of the matrix element is most important in the long-wavelength limit, so that we explicitly isolate the long-range part in the following.

We will use the e–pn interaction matrix element in the long-wavelength limit as is customary in semiconductor spintronics [28]. Grimaldi and Fulde [23] also demonstrate how the long-range contribution is screened so that ${v}_{\mathrm{eff}}({\bf{q}})$ has a well defined ${\bf{q}}\to 0$ (long-wavelength) limit. In this limit, the long-range interaction operator is

$$\begin{eqnarray}&&{\hat{v}}_{{\bf{k}}+{\bf{q}},{\bf{k}}}^{(\lambda )}={{\ell }}_{{\bf{q}},\lambda }^{(0)}{v}_{\mathrm{eff}}({\bf{q}}\to 0)({\bf{q}}\cdot {{\boldsymbol{\varepsilon }}}_{{\bf{q}},\lambda })\{\underset{(1)}{\underbrace{1}}-\underset{(2)}{\underbrace{{\rm{i}}\xi [{\bf{q}}\times (\hat{{\bf{p}}}+{\hslash }{\bf{k}})]\cdot \hat{{\bf{s}}}}}\}.\end{eqnarray} \tag{ 23 }$$

Note that for small q the electrons couple exclusively to longitudinal phonons, where the polarization vector ${{\boldsymbol{\varepsilon }}}_{{\bf{q}},\lambda }$ points in the same direction as the wave-vector ${\bf{q}}$.

In this section we derive the equation of motion for the reduced spin–density matrix

$$\begin{eqnarray}&&{\rho }_{{\bf{k}}}^{\mu {\mu }^{\prime }}\equiv \langle {\hat{c}}_{\mu {\bf{k}}}^{\dagger }{\hat{c}}_{{\mu }^{\prime }{\bf{k}}}\rangle \end{eqnarray} \tag{ 24 }$$

which determines all single-particle properties of the electronic ensemble.

We will for definiteness also include coherent terms and eventually take the scattering limit. Therefore we need the single-particle hamiltonian (12) in diagonal form

$$\begin{eqnarray}&&{H}_{{\rm{e}}}=\displaystyle \sum _{\mu }{\epsilon }_{\mu {\bf{k}}}{c}_{\mu {\bf{k}}}^{\dagger }{c}_{\mu {\bf{k}}}.\end{eqnarray} \tag{ 25 }$$

Here, we ${\bf{k}}\cdot {\bf{p}}$ theory [32, 33], i.e., the ${u}_{\mu {\bf{k}}}$'s and energy dispersions ${\epsilon }_{\mu {\bf{k}}}$ are determined as the solution of the eigenvalue problem of the self-adjoint ${\bf{k}}\cdot {\bf{p}}$ operator ${\hat{h}}_{\mathrm{eff}}({\bf{k}})$

$$\begin{eqnarray}&&{\hat{h}}_{\mathrm{eff}}({\bf{k}}){u}_{\mu {\bf{k}}}={\epsilon }_{\mu {\bf{k}}}{u}_{\mu {\bf{k}}}\end{eqnarray} \tag{ 26 }$$

for lattice-periodic ${u}_{\mu {\bf{k}}}$'s on Ω. This guarantees that, for each ${\bf{k}}$, the eigenvalue spectrum ${\epsilon }_{\mu {\bf{k}}}$ is discrete and the ${u}_{\mu {\bf{k}}}$ form a complete set, i.e.

$$\begin{eqnarray}&&\displaystyle \sum _{\mu }| {u}_{\mu {\bf{k}}}\rangle \langle {u}_{\mu {\bf{k}}}| ={\mathbb{1}}.\end{eqnarray} \tag{ 27 }$$

Further, the ${u}_{\mu {\bf{k}}}$'s for each ${\bf{k}}$ are orthogonal with regard to the scalar product (14). In semiconductors, the ${\bf{k}}\cdot {\bf{p}}$ hamiltonian operator can be approximated on a subset of bands by an effective hamiltonian, such as the Luttinger hamiltonian [33]. The number of bands included in (25) determines the dimension of the reduced density matrix ${\rho }_{{\bf{k}}}$.

The equation of motion for the electronic spin–density matrix ${\rho }_{{\bf{k}}}^{\mu {\mu }^{\prime }}$ due to (20) and  (25) is

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{\rho }_{{\bf{k}}}^{\mu {\mu }^{\prime }}=\displaystyle \frac{{\rm{i}}}{{\hslash }}({\epsilon }_{\mu {\bf{k}}}-{\epsilon }_{{\mu }^{\prime }{\bf{k}}}){\rho }_{{\bf{k}}}^{\mu {\mu }^{\prime }}+{\left.\displaystyle \frac{\partial }{\partial t}{\rho }_{{\bf{k}}}^{\mu {\mu }^{\prime }}\right|}_{{\rm{e}}-\mathrm{pn}},\end{eqnarray} \tag{ 28 }$$

with

$$\begin{eqnarray}{\left.\displaystyle \frac{\partial }{\partial t}{\rho }_{{\bf{k}}}^{\mu {\mu }^{\prime }}\right|}_{{\rm{e}}-\mathrm{pn}} & = & \displaystyle \frac{{\rm{i}}}{{\hslash }}\displaystyle \sum _{\lambda }\displaystyle \sum _{\nu {\bf{q}}}\{{\langle {u}_{\nu {\bf{k}}+{\bf{q}}}| {\hat{v}}_{{\bf{k}}+{\bf{q}},{\bf{k}}}^{(\lambda )}{u}_{\mu {\bf{k}}}\rangle }_{{\rm{\Omega }}}\langle {B}_{{\bf{q}},\lambda }{c}_{\nu {\bf{k}}+{\bf{q}}}^{\dagger }{c}_{{\mu }^{\prime }{\bf{k}}}\rangle \\ & & -{\langle {u}_{{\mu }^{\prime }{\bf{k}}}| {\hat{v}}_{{\bf{k}},{\bf{k}}+{\bf{q}}}^{(\lambda )}{u}_{\nu {\bf{k}}+{\bf{q}}}\rangle }_{{\rm{\Omega }}}\langle {B}_{-{\bf{q}},\lambda }{c}_{\mu {\bf{k}}}^{\dagger }{c}_{\nu {\bf{k}}+{\bf{q}}}\rangle \}.\end{eqnarray} \tag{ 29 }$$

Equation (28) is the full equation of motion that contains the coherent contribution, which includes a rotation of the off-diagonal elements of the spin–density matrix with energy difference ${\epsilon }_{\mu {\bf{k}}}-{\epsilon }_{{\mu }^{\prime }{\bf{k}}}$, and the e–pn interaction contribution (29), in which the phonon-assisted electronic density matrix $\langle {b}_{{\bf{q}},\lambda }{c}_{\nu {\bf{k}}+{\bf{q}}}^{\dagger }{c}_{{\mu }^{\prime }{\bf{k}}}\rangle$ appears [33, 34]. Equation (28) also describes spin relaxation due to the coherent rotation of ρ because the single-particle hamiltonian (26) includes spin–orbit coupling and thus may lead to spin–split bands. For the case of hole bands in GaAs, the influence of coherent rotation times has been studied, e.g., [27, 35, 36]. In the context of ferromagnets, the importance of the off-diagonal elements of the spin–density matrix has been also investigated in the framework of a quantum kinetic calculation recently [37]. In accordance with [27] we obtain the pure EY mechanism by neglecting the spin splitting, i.e., coherent contributions in (28), in the basis in which the single-particle hamiltonian (26) is diagonal. For the case of a pair of Kramers degenerate bands, to which we specialize in section 4 below, the coherent contributions vanish exactly.

Before we make approximations to the phonon-assisted density matrix, we compute the dynamics of the vector of the average spin

$$\begin{eqnarray}&&{\bf{S}}=\left\langle {\displaystyle \int }_{{ \mathcal V }}{{\rm{\Psi }}}^{\dagger }({\bf{r}})\hat{{\bf{s}}}\;{\rm{\Psi }}({\bf{r}}){{\rm{d}}}^{3}r\right\rangle =\displaystyle \sum _{{\bf{k}}}\displaystyle \sum _{\mu {\mu }^{\prime }}{\langle {u}_{\mu {\bf{k}}}| \hat{{\bf{s}}}{u}_{{\mu }^{\prime }{\bf{k}}}\rangle }_{{\rm{\Omega }}}\ {\rho }_{{\bf{k}}}^{\mu {\mu }^{\prime }}.\end{eqnarray} \tag{ 30 }$$

By combining (29) and (30) we obtain for the $\alpha =x$, y, z components

$$\begin{eqnarray}{\left.\displaystyle \frac{\partial }{\partial t}{S}_{\alpha }\right|}_{{\rm{e}}-\mathrm{pn}} & = & \displaystyle \frac{{\rm{i}}}{{\hslash }}\displaystyle \sum _{\lambda }\displaystyle \sum _{\nu \mu {\mu }^{\prime }}\displaystyle \sum _{{\bf{k}},{\bf{q}}}{\langle {u}_{\mu {\bf{k}}}| {\hat{s}}_{\alpha }{u}_{{\mu }^{\prime }{\bf{k}}}\rangle }_{{\rm{\Omega }}}\{{\langle {u}_{\nu {\bf{k}}+{\bf{q}}}| {\hat{v}}_{{\bf{k}}+{\bf{q}},{\bf{k}}}^{(\lambda )}{u}_{\mu {\bf{k}}}\rangle }_{{\rm{\Omega }}}\\ & & \times \langle {B}_{{\bf{q}},\lambda }{c}_{\nu {\bf{k}}+{\bf{q}}}^{\dagger }{c}_{{\mu }^{\prime }{\bf{k}}}\rangle -{\langle {u}_{{\mu }^{\prime }{\bf{k}}}| {\hat{v}}_{{\bf{k}},{\bf{k}}+{\bf{q}}}^{(\lambda )}{u}_{\nu {\bf{k}}+{\bf{q}}}\rangle }_{{\rm{\Omega }}}\langle {B}_{-{\bf{q}},\lambda }{c}_{\mu {\bf{k}}}^{\dagger }{c}_{\nu {\bf{k}}+{\bf{q}}}\rangle \}.\end{eqnarray} \tag{ 31 }$$

We can employ the completeness property (27) of the ${u}_{\mu {\bf{k}}}$'s to show that this has the form

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\partial }{\partial t}{S}_{\alpha }\right|}_{{\rm{e}}-\mathrm{pn}}=\displaystyle \sum _{{\bf{kq}}}\displaystyle \sum _{\mu {\mu }^{\prime },\lambda }\langle {B}_{{\bf{q}},\lambda }{c}_{{\mu }^{\prime }{\bf{k}}+{\bf{q}}}^{\dagger }{c}_{\mu {\bf{k}}}\rangle {\langle {u}_{{\mu }^{\prime }{\bf{k}}+{\bf{q}}}| {({\hat{{\bf{t}}}}_{{\bf{k}}+{\bf{q}},{\bf{k}}}^{(\lambda )})}_{\alpha }{u}_{\mu {\bf{k}}}\rangle }_{{\rm{\Omega }}},\end{eqnarray} \tag{ 32 }$$

with the torque vector operator due to the e–pn interaction

$$\begin{eqnarray}&&{\hat{{\bf{t}}}}_{{\bf{k}}+{\bf{q}},{\bf{k}}}^{(\lambda )}:= \displaystyle \frac{1}{{\rm{i}}{\hslash }}[\hat{{\bf{s}}},{\hat{v}}_{{\bf{k}}+{\bf{q}},{\bf{k}}}^{(\lambda )}].\end{eqnarray} \tag{ 33 }$$

We will denote the matrix element of this operator by

$$\begin{eqnarray}&&{{\bf{t}}}_{{\mu }^{\prime }{\bf{k}}+{\bf{q}},\mu {\bf{k}}}^{(\lambda )}:= \langle {u}_{{\mu }^{\prime }{\bf{k}}+{\bf{q}}}| {\hat{{\bf{t}}}}_{{\bf{k}}+{\bf{q}},{\bf{k}}}^{(\lambda )}{u}_{\mu {\bf{k}}}\rangle .\end{eqnarray} \tag{ 34 }$$

This matrix element completely determines the spin change that occurs in a transition $(\mu {\bf{k}})\to ({\mu }^{\prime }{\bf{k}}+{\bf{q}})$. This information is not explicit in the matrix element (21). Using (22), one finds

$$\begin{eqnarray}&&{\hat{{\bf{t}}}}_{{\bf{k}}+{\bf{q}},{\bf{k}}}^{(\lambda )}=\displaystyle \sum _{{\bf{G}}}{{\ell }}_{{\bf{q}},\lambda }^{(0)}{{\rm{e}}}^{{\rm{i}}{\bf{G}}\cdot {\bf{r}}}{v}_{\mathrm{eff}}({\bf{q}}+{\bf{G}})\\ &&\quad [({\bf{q}}+{\bf{G}})\cdot {{\boldsymbol{\varepsilon }}}_{{\bf{q}},\lambda }]\{\underset{(1)}{\underbrace{0}}+\underset{(2)}{\underbrace{(-{\rm{i}}\xi )[({\bf{q}}+{\bf{G}})\times ({\hslash }{\bf{k}}+\hat{{\bf{p}}})]\times \hat{{\bf{s}}}}}\},\end{eqnarray} \tag{ 35 }$$

where, again, the two terms result from spin–diagonal and explicitly spin-dependent contributions (15) and (16), respectively. The expression (32) together with (35) is an important result, as it shows that only the explicitly spin dependent contribution to the e–pn interaction gives rise to a non-vanishing torque (35). This conclusion does not depend on the single-particle basis ${u}_{\mu {\bf{k}}}$ in the sense that one could apply a ${\bf{k}}$-dependent unitary transformation to the ${u}_{\mu {\bf{k}}}$'s and the ${c}_{\mu {\bf{k}}}$'s in (32) and (35). Note that no approximation whatsoever has been made for the dynamics of the phonon assisted quantity. Further, equation (32) and the definition of the torque matrix element is a rather general result for the spin change due to any incoherent electron–boson scattering, for instance, with photons or magnons.

We want to illustrate this result in the language of semiconductor spintronics. Here, the Elliott spin relaxation mechanism refers to the spin relaxation due to a spin–diagonal interaction combined with a spin-mixing in the wavefunction. Our result shows that it is not the spin mixing in the wave functions that determines the Elliott torque matrix element: in this restricted sense, there is no Elliott contribution to spin dynamics in any basis, and the spin dynamics originates exclusively from the spin-dependent Overhauser (or Yafet) contribution. When we evaluate the phonon-assisted density matrix in the scattering limit, an electron–phonon matrix element appears, which contains both spin–diagonal and spin-dependent contributions. For this regular matrix element, the spin diagonal ('Elliott') contribution is likely much bigger than the spin-dependent ('Overhauser' or 'Yafet') contribution.

The torque operator ${\hat{{\bf{t}}}}_{{\bf{k}}+{\bf{q}},{\bf{k}}}^{(\lambda )}$, which is due to the spin-dependent contribution only, has the long-wavelength limit

$$\begin{eqnarray}&&{\hat{{\bf{t}}}}_{{\bf{k}}+{\bf{q}},{\bf{k}}}^{(\lambda )}=-{\rm{i}}\xi {{\ell }}_{{\bf{q}},\lambda }^{(0)}{v}_{\mathrm{eff}}({\bf{q}}\to 0)({\bf{q}}\cdot {{\boldsymbol{\varepsilon }}}_{{\bf{q}},\lambda })[{\bf{q}}\times ({\hslash }{\bf{k}}+\hat{{\bf{p}}})]\times \hat{{\bf{s}}}.\end{eqnarray} \tag{ 36 }$$

To compare the results of the present approach with the conventional EY analysis, we specialize the dynamical equation (32) for ${S}_{\alpha }$ for the case of incoherent scattering. To this end, we need to evaluate the phonon-assisted density matrix $\langle {b}_{{\bf{q}},\lambda }{c}_{\mu {\bf{k}}+{\bf{q}}}^{\dagger }{c}_{{\mu }^{\prime }{\bf{k}}}\rangle$. Its equation of motion couples to higher order correlation functions of the type $\langle {b}^{\dagger }{{bc}}^{\dagger }c\rangle$ and $\langle {c}^{\dagger }{c}^{\dagger }{cc}\rangle$. This equation-of-motion hierarchy is truncated at the scattering level according to the dynamical coupled-cluster method, and by performing a Markov approximation, see [34, 38]. Details of this procedure are given in appendix B. With these approximations we find for the ensemble averaged spin

$$\begin{eqnarray}\displaystyle \frac{\partial }{\partial t}{\bf{S}} & = & 2\mathrm{Re}\displaystyle \sum _{{\bf{kq}}}\displaystyle \sum _{\mu {\mu }^{\prime }}\displaystyle \sum _{\lambda }{{\bf{t}}}_{{\mu }^{\prime }{\bf{k}}+{\bf{q}},\mu {\bf{k}}}^{(\lambda )}{g}_{{\tau }^{\prime }{\bf{k}},\tau {\bf{k}}+{\bf{q}}}^{(\lambda )}\displaystyle \frac{1}{({\epsilon }_{{\mu }^{\prime }{\bf{k}}+{\bf{q}}}-{\epsilon }_{\mu {\bf{k}}}-{\hslash }{\omega }_{{\bf{q}},\lambda })+{\rm{i}}{\hslash }\gamma }\\ & & \quad \times \displaystyle \sum _{\tau {\tau }^{\prime }}\{[{\delta }_{\mu {\tau }^{\prime }}-{\rho }_{{\bf{k}}}^{{\tau }^{\prime }\mu }]{\rho }_{{\bf{k}}+{\bf{q}}}^{{\mu }^{\prime }\tau }(1+{N}_{{\bf{q}},\lambda })-[{\delta }_{\tau {\mu }^{\prime }}-{\rho }_{{\bf{k}}+{\bf{q}}}^{{\mu }^{\prime }\tau }]{\rho }_{{\bf{k}}}^{{\tau }^{\prime }\mu }{N}_{{\bf{q}},\lambda }\}.\end{eqnarray} \tag{ 37 }$$

where ${N}_{{\bf{q}},\lambda }=b({\hslash }{\omega }_{{\bf{q}},\lambda })$ are phonon occupation numbers. We will assume that the phonon system is not changed appreciably by the electronic dynamics so that the ${N}_{{\bf{q}},\lambda }=b({\hslash }{\omega }_{{\bf{q}},\lambda }$) are the equilibrium phonon occupations given by the Bose function $b({\hslash }{\omega }_{{\bf{q}},\lambda })$. In the result (37) the full reduced density matrix appears on the right-hand side, which needs to be determined to obtain the spin dynamics.

Kramers degenerate states including spin–orbit coupling can be assumed to have the form [39–41]

$$\begin{eqnarray}&&{u}_{{\bf{k}}\Uparrow }({\bf{x}})={a}_{{\bf{k}}}({\bf{x}})| \uparrow \rangle +{b}_{{\bf{k}}}({\bf{x}})| \downarrow \rangle ,\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}&&{u}_{{\bf{k}}\Downarrow }({\bf{x}})={a}_{-{\bf{k}}}^{*}({\bf{x}})| \downarrow \rangle -{b}_{-{\bf{k}}}^{*}({\bf{x}})| \uparrow \rangle .\end{eqnarray} \tag{ 39 }$$

In semiconductors and metals with relatively weak spin–orbit coupling, these spin-mixed states can be considered to be mostly spin-up and spin-down for the vast majority of the Brillouin zone, which is indicated by the labels ⇑, ⇓. More precisely, $\displaystyle \frac{1}{{\rm{\Omega }}}{\displaystyle \int }_{{\rm{\Omega }}}{| {b}_{{\bf{k}}}({\bf{x}})| }^{2}\;{{\rm{d}}}^{3}x$ is much smaller than $\displaystyle \frac{1}{{\rm{\Omega }}}{\displaystyle \int }_{{\rm{\Omega }}}{| {a}_{{\bf{k}}}({\bf{x}})| }^{2}\;{{\rm{d}}}^{3}x$. Because of the degeneracy, there is an ambiguity in the definition of these states, as any superposition will also be an eigenstate. In accordance with [39–41], we choose them to fulfill

$$\begin{eqnarray}&&\langle {u}_{\Uparrow ,{\bf{k}}}| {\hat{s}}_{z}{u}_{\Downarrow ,{\bf{k}}}\rangle =0.\end{eqnarray} \tag{ 40 }$$

This is equivalent to the statement that the chosen states ${u}_{\Uparrow {\bf{k}}}$, ${u}_{\Downarrow {\bf{k}}}$ diagonalize ${\hat{s}}_{z}$ in the degenerate subspace. The choice is important, if one wants to attach some importance to the magnitude of $\displaystyle \frac{1}{{\rm{\Omega }}}{\displaystyle \int }_{{\rm{\Omega }}}{| {b}_{{\bf{k}}}({\bf{x}})| }^{2}\;{{\rm{d}}}^{3}x$ as spin mixing parameter. In ab initio calculations one often uses a small external magnetic field to enforce a quantization direction so that (40) is fulfilled.

The condition (40) does not imply that the electrons are in pure spin-up or spin-down states characterized by the same two-dimensional spinors at each ${\bf{k}}.$ In fact, we have

$$\begin{eqnarray}&&{u}_{\Uparrow ,{\bf{k}}}\ne | \uparrow \rangle \qquad \mathrm{and}\qquad {u}_{\Downarrow ,{\bf{k}}}\ne | \downarrow \rangle ,\end{eqnarray} \tag{ 41 }$$

as well as $| \langle {u}_{\lambda {\bf{k}}}| {\hat{s}}_{z}{u}_{\lambda {\bf{k}}}\rangle | \leqslant {\hslash }/2$ for $\lambda =\Uparrow$, ⇓. For more details, see appendix A.

In order to derive a characteristic rate for the relaxation of a small excess spin polarization $\delta {S}_{z}$ at a temperature T, we first need to determine the reduced density matrix for this case. We assume the system to be in a quasi-equilibrium with a given spin-expectation value $\delta {S}_{z}$. Thus we must determine the quasi-equilibrium distribution for the eigenenergies of the generalized grand-canonical ${\bf{k}}\cdot {\bf{p}}$ hamiltonian

$$\begin{eqnarray}&&\hat{{\mathfrak{k}}}({\bf{k}})\equiv {\hat{h}}_{\mathrm{eff}}({\bf{k}})-\mu -{\zeta }_{z}{\hat{s}}_{z}.\end{eqnarray} \tag{ 42 }$$

Here, ${\zeta }_{z}$, the spin accumulation, is the Lagrange parameter needed to obtain the finite spin expectation value. Since ${\hat{h}}_{\mathrm{eff}}({\bf{k}})$ includes spin–orbit coupling it does not commute with ${\hat{s}}_{z}$, so that the eigenstates and eigenenergies of

$$\begin{eqnarray}&&\hat{{\mathfrak{k}}}({\bf{k}}){w}_{\mu {\bf{k}}}={E}_{\mu {\bf{k}}}{w}_{\mu {\bf{k}}}\end{eqnarray} \tag{ 43 }$$

are, in general, not identical to the ${u}_{\mu {\bf{k}}}$'s that result as eigenfunctions of ${\hat{h}}_{\mathrm{eff}}({\bf{k}})$, see equation (26).

We assume that ${\zeta }_{z}$ is small and determine ${E}_{\mu {\bf{k}}}$ perturbatively in the degenerate subspace of the ${u}_{\mu {\bf{k}}}$'s. Since we chose the ${u}_{\mu {\bf{k}}}$'s that diagonalize ${\hat{s}}_{z}$ in the degenerate subspace, see (40), this gives

$$\begin{eqnarray}&&{E}_{\mu {\bf{k}}}\simeq {\epsilon }_{{\bf{k}}}-\mu -{\zeta }_{z}\langle {u}_{\mu {\bf{k}}}| {\hat{s}}_{z}{u}_{\mu {\bf{k}}}\rangle .\end{eqnarray} \tag{ 44 }$$

If we also assume that other bands are sufficiently far away, the eigenvectors are ${w}_{\mu {\bf{k}}}\simeq {u}_{\mu {\bf{k}}}$. If there are other bands close to the bands for which the spin relaxation is computed, one needs to include an eigenvector correction for ${w}_{\mu {\bf{k}}}$.

The quasi-equilibrium reduced density matrix in the eigenbasis of the grand canonical hamiltonian is diagonal and depends only on the grand-canonical energies ${E}_{\mu {\bf{k}}}$, i.e.,

$$\begin{eqnarray}&&{\rho }_{{\bf{k}}}^{\mu {\mu }^{\prime }}={\delta }_{\mu ,{\mu }^{\prime }}f({E}_{\mu {\bf{k}}}),\end{eqnarray} \tag{ 45 }$$

where $f(\epsilon )={[\mathrm{exp}(\beta \epsilon )+1]}^{-1}$ is the Fermi function. Since ${\zeta }_{z}$ is small, we keep only the first order result in ${\zeta }_{z}$ and with (44) we find for the quasi-equilibrium distribution

$$\begin{eqnarray}&&{\rho }_{{\bf{k}}}^{\mu {\mu }^{\prime }}={\delta }_{\mu ,{\mu }^{\prime }}[f({\epsilon }_{{\bf{k}}}-\mu )+{{\rm{\Delta }}}_{{\bf{k}}}{\zeta }_{z}\langle {u}_{\mu {\bf{k}}}| {\hat{s}}_{z}{u}_{\mu {\bf{k}}}\rangle ],\end{eqnarray} \tag{ 46 }$$

where we have defined ${{\rm{\Delta }}}_{{\bf{k}}}\equiv -{\displaystyle \frac{{\rm{d}}f}{{\rm{d}}\epsilon }| }_{{\epsilon }_{{\bf{k}}}-\mu }$. For low temperatures, this function approaches a δ-function concentrated at ${\epsilon }_{{\bf{k}}}-\mu$. From the reduced density matrix (46) we find the relation ${\zeta }_{z}=\delta {S}_{z}/{ \mathcal N }$ between $\delta {S}_{z}$ and ${\zeta }_{z}$ with the normalization

$$\begin{eqnarray}&&{ \mathcal N }=\displaystyle \sum _{\mu {\bf{k}}}{\langle {u}_{\mu {\bf{k}}}| {\hat{s}}_{z}{u}_{\mu {\bf{k}}}\rangle }^{2}{{\rm{\Delta }}}_{{\bf{k}}}.\end{eqnarray} \tag{ 47 }$$

To obtain the close-to-equilibrium spin relaxation rate defined by

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}\delta {S}_{z}=-\displaystyle \frac{\delta {S}_{z}}{{\tau }_{\mathrm{SR}}},\end{eqnarray} \tag{ 48 }$$

we insert (46) in (37) and linearize. Then, as in Yafet's original treatment and its extension by Fähnle and coworkers, [20] ${\zeta }_{z}$ drops out and we find

$$\begin{eqnarray}\displaystyle \frac{1}{{\tau }_{\mathrm{SR}}} & = & -\displaystyle \frac{2}{{ \mathcal N }}\mathrm{Re}\displaystyle \sum _{{\bf{kq}}}\displaystyle \sum _{\mu {\mu }^{\prime }}\displaystyle \sum _{\lambda }{({t}_{{\mu }^{\prime }{\bf{k}}+{\bf{q}},\mu {\bf{k}}}^{(\lambda )})}_{z}{g}_{\mu {\bf{k}},{\mu }^{\prime }{\bf{k}}+{\bf{q}}}^{(\lambda )}\langle {u}_{\mu {\bf{k}}+{\bf{q}}}| {\hat{s}}_{z}{u}_{\mu {\bf{k}}+{\bf{q}}}\rangle {{\rm{\Delta }}}_{{\bf{k}}+{\bf{q}}}\\ & & \times \left\{\displaystyle \frac{{N}_{{\bf{q}},\lambda }+f({\epsilon }_{\mu {\bf{k}}}-\mu )}{{\epsilon }_{{\mu }^{\prime }{\bf{k}}+{\bf{q}}}-{\epsilon }_{\mu {\bf{k}}}+{\hslash }{\omega }_{{\bf{q}},\lambda }+{\rm{i}}{\hslash }\gamma }+\displaystyle \frac{1+{N}_{{\bf{q}},\lambda }-f({\epsilon }_{\mu {\bf{k}}}-\mu )}{{\epsilon }_{{\mu }^{\prime }{\bf{k}}+{\bf{q}}}-{\epsilon }_{\mu {\bf{k}}}-{\hslash }{\omega }_{{\bf{q}},\lambda }+{\rm{i}}{\hslash }\gamma }\right\}.\end{eqnarray} \tag{ 49 }$$

This expression for the spin relaxation, or T₁, time is a more physically transparent result compared to Yafet's because it relates the spin dynamics directly to a torque matrix element. It improves on the original Yafet result by including the correct bookkeeping for spin instead of counting the contribution of spin mixed states as if they were pure spin states. Instead of ${| {g}_{\mu {\bf{k}},{\mu }^{\prime }{{\bf{k}}}^{\prime }}^{(\lambda )}| }^{2}$, i.e., the squared modulus of the e–pn matrix element, as it occurs in Yafet's result, here matrix elements of the torque, the electron–phonon interaction and the spin appear.

In Yafet's original result, only one matrix element, namely our g appears, and the spin diagonal and explicitly spin-dependent contributions to this matrix element are conventionally called, respectively, the Elliott and Overhauser (or Yafet) contributions. In our derivation, two different matrix elements occur. If we stick to the established language, the torque matrix element only has an Overhauser contribution, whereas in the g matrix element both Elliott and Overhauser terms contribute. Thus we do not have a 'pure' Elliott contribution, but the leading term is a mixed one with an Elliott contribution in g, but Overhauser in t. Regarding the g matrix element in the long wavelength limit, Yafet [5] and numerical evaluations [41] of his spin-relaxation result find a cancellation between Elliott and Overhauser contributions in the short range contribution, whereas Grimaldi and Fulde [23] find for the long-wavelength limit a dominating long-range Elliott contribution to g.

Our result (49) is not limited to small spin mixing and is qualitatively different from the original Yafet formula. It will make it easier to compute accurate EY spin relaxation times numerically and to compare different contributions to spin dynamics by computing the relevant combination of torque and g matrix elements.