Reexamination of the Elliott-Yafet spin-relaxation mechanism

We analyze spin-dependent carrier dynamics due to incoherent electron-phonon scattering, which is commonly referred to as Elliott-Yafet (EY) spin-relaxation mechanism. For this mechanism one usually distinguishes two contributions: (1) from the electrostatic interaction together with spin-mixing in the wave functions, which is often called the Elliott contribution, and (2) the phonon-modulated spin-orbit interaction, which is often called the Yafet or Overhauser contribution. By computing the reduced electronic density matrix, we improve Yafet's original calculation, which is not valid for pronounced spin mixing as it equates the pseudo-spin polarization with the spin polarization. The important novel quantity in our calculation is a torque operator that determines the spin dynamics. The contribution (1) to this torque vanishes exactly. From this general result, we derive a modified expression for the Elliott-Yafet spin relaxation time.


I. INTRODUCTION
Spin relaxation plays a role in spin-dependent dynamics both on long and short timescales.
While its contribution to (precessional) magnetization damping is usually ascribed to spinlattice relaxation and treated phenomenologically, the analysis of ultrafast demagnetization dynamics has often been based on the microscopic concept of spin-flip processes due to electron-phonon interactions as developed for semiconductors in the 1950s. 1 Overhauser 2 was the first to identify the modulation of the spin-orbit interaction by lattice vibrations as the microscopic mechanism for spin relaxation due to incoherent electron-phonon scattering processes. Elliott 3 argued shortly thereafter that there is an additional contribution to the spin relaxation due to the momentum-dependent spin mixing in the wave functions and that, consequently, even spin-diagonal incoherent scattering processes due to spin-diagonal interactions can contribute to spin relaxation. The subject was taken up again by Gerasimenko and Andreev, 4 and Yafet. 5 The latter calculated the spin-flip matrix element due to electron-phonon interaction (as a function of electronic momentum transfer q) including both the Overhauser and the Elliott contributions and showed that the first few orders in q vanish due to a cancellation of Overhauser and Elliott contributions. Nowadays, this combination of Overhauser and Elliott contributions is usually called Elliott-Yafet mechanism because Yafet derived a relatively simple result for the close-to-equilibrium spin relaxation time, which is suitable for evaluation from ab-initio input and can be used to justify an approximate relation between the spin relaxation time and the momentum relaxation time.
At present, both simplified and ab-initio based expressions due to the EY mechanism are widely used, in particular for magnetization dynamics in metals. [6][7][8][9][10] In this paper, we present a new analysis of the spin-relaxation problem due to incoherent electron-phonon scattering, as it was originally considered by Overhauser, Elliott, and Yafet.
We derive the dynamical equation for the change of the reduced electronic spin density matrix by expressing the spin dynamics in terms of phonon assisted density matrices. This approach achieves a correct description of the dynamics of the spin vector, as opposed to the Yafet treatment, which gets spins and pseudospins mixed up, and thus cannot correctly account for the "amount of spin-flip" in each scattering transition. If one correctly describes the spin vector, it becomes obvious that the important quantity for spin dynamics is a torque matrix element, which is not present in the conventional derivation. 5 Remarkably, there is no contribution to this torque matrix element from spin-diagonal scattering mechanisms. Put differently, spin-diagonal electron-phonon scattering and spin-orbit coupling alone, which is usually referred to as the Elliott spin-relaxation mechanism, yields no spin dynamics. Based on this observation, we derive a modified result for the EY spin relaxation time.
This paper is organized as follows. In Sec. II we briefly review the conventional Elliott-Yafet treatment, some more recent contributions, and how the Elliott-Yafet spin relaxation mechanism has been applied in theoretical models for the demagnetization of ferromagnets.
In Sec. III we set up the electron-phonon interaction hamiltonian, discuss long and shortrange contributions, and derive the equations of motion for the spin density matrix and the spin expectation value. Section IV is devoted to the derivation of a spin relaxation time for the special case of Kramers degenerate bands. An important ingredient for this derivation is the form of the quasi-equilibrium spin-density matrix in the presence of spin-orbit coupling, which is discussed in some detail. The conclusions are presented in Sec. V, and appendix A contains a short demonstration concerning the form of the spinor wavefunctions on Kramers degenerate bands.

II. ELLIOTT-YAFET APPROACH
The Elliott-Yafet approach has been reviewed often. 1,11,12 Only for comparison with our calculations, we repeat here some of the main points of Yafet's derivation in Ref. 5 using his notation. The objective of Yafet is to calculate the rate of a spin-flip transition for two Kramers degenerate bands including spin-orbit coupling. The Kramers degeneracy implies that for a band index b we have two wave functions with the same energy E bk , whereK is the time-reversal operator. Focusing on one band and dropping the corresponding b index he then calculates the Golden-Rule transition probability W ⇑k,⇓k ′ for a spin-flip transition ψ ⇑k → ψ ⇓k due to the electron-phonon interaction. More precisely, Yafet calculates the dynamics of the spin polarization, which is defined as where n λk denote the carrier distributions with momentum k and pseudospins λ =⇑, ⇓. It is clear that D spin cannot be a good approximation to the spin polarization for pronounced spin mixing. Then (including a factor of 2 as in Yafet's derivation) one obtains Thus the change of spin polarization is determined essentially by the number of transitions (per unit time), which is obtained by adding the in-scattering and out-scattering Golden-Rule probabilities where and q = |k ′ − k|. The other symbols have an obvious meaning and are defined below. To obtain a relaxation time valid for a small spin polarization, the distribution functions n λk are assumed to be of the quasi-equilibrium form f (ǫ µk −µ λ ) with µ λ the (pseudo)spin dependent chemical potentials, and the distributions are expanded for small µ λ 's. This treatment has recently been extended to ferromagnets. 13 Finally, Yafet shows from a symmetry argument that there is a cancellation between the two contributions to the electron-phonon (e-pn) coupling matrix element in the long-wavelength limit. In this paper, we refer to these two contributions as the Elliott and Overhauser contributions, which we explain in detail below.
In our view this method has three problems.
1) The spin polarization D spin , defined in (3), is computed as the difference of pseudospin occupation numbers as if the electrons were in pure spin states for all k values. Here, the Yafet derivation ignores spin mixing, which leads to a k dependent expectation value ofŝ z . Importantly, the modulus of this spin expectation value may be significantly smaller than /2.
2) The calculation of spin dynamics is based on transition probabilities between pseudospin states or, equivalently, pseudo-spin occupation numbers n ⇑k and n ⇓k . Such a treatment neglects coherences between the pseudospin states, which are the offdiagonal elements of the spin density matrix 3) Yafet proves the cancellation between the Elliott and Overhauser contributions to the spin-flip matrix element M ⇑k,⇓k ′ for the short range part of the electron-phonon coupling matrix element. As shown by Grimaldi and Fulde 14 there is also a long-range contribution, for which this cancellation does not hold, and which is larger than the short-range contribution in the long-wave limit.
The Yafet method is so widely accepted that it may be worthwhile to mention that ours is not the first paper to point to these problems. For instance, Ref. 13 states Obviously, m λ k [the magnetic moment of state (λ, k)] is different for different k-vectors, but it is about ±µ B as long as the spin mixing is small. Yafet neglects the k-dependence of m λ k and assumes a constant magnetic moment m or −m for all dominant spin-up or all dominant spin-down states, respectively, which is certainly problematic in systems with spin "hot spots," i.e., with regions in the Brillouin zone where the spin mixing is very large.
Also, for semiconductors, Yu et al. 15 have used spin projection operators to avoid the first problem in their calculation of the EY spin-relaxation time, but missed the Overhauser contribution, as they only include spin-diagonal scattering processes.
Going beyond spin-relaxation times, there are microscopic approaches to spin-dependent carrier dynamics that compute the reduced density matrix and thus avoid using the pseudospin polarization, [16][17][18][19] see Ref. 20 for a review. Numerically solving for the full dynamical spin density matrix also yields Dyakonov-Perel and Bir-Aronov-Pikus contributions spin relaxation, but is very CPU-time intensive and difficult because the microscopic carrier dynamics has to be computed to obtain spin relaxation times that are often orders of magnitude longer than typical scattering times, which have to be resolved in the numerics.
Numerically calculating the spin density matrix does not yield explicit expressions for spin relaxation times, which is one of the goals of our paper. Further, our calculation gives a more transparent description of the ensemble spin dynamics in terms of the torque matrix element than is possible by "brute-force" calculation of the dynamical spin-density matrix.
The Elliott-Yafet mechanism has also been applied to the demagnetization dynamics in ferromagnets. [6][7][8][9][10]21 None of these approaches avoids all of the three problems above.

III. SPIN-DEPENDENT ELECTRON-PHONON SCATTERING DYNAMICS
In this section we present a derivation of spin-dependent carrier dynamics due to electronphonon interactions that avoids the three problems listed above. We first derive the interaction hamiltonian between electrons and phonons including the long-range contribution of the Coulomb potential, and specialize to the long-wavelength limit. With this interaction we derive the equation of motion for the spin density matrix including phonon-assisted correlation functions. We identify the torque matrix element that determines the incoherent dynamics of the average spin. Finally, we derive the scattering limit of the dynamical equation.

A. Electron-phonon interaction hamiltonian
We start our derivation by writing down the interaction potential of a single electron in a lattice in terms of lattice site coordinates R n within the rigid-ion approximation aŝ Here, v eff is the effective electrostatic Coulomb interaction between the electron and the ionic core and ξ = /(4m 2 c 2 ). We denote single-particle operators acting on the electronic space and spin variables by small letters and a hat; for instance,p = −i ∇. 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 operatorx ×p of the electrons, and therefore does not conserve the electronic orbital angular momentum. The rigid-ion approximation places no fundamental restriction on the following development, but the equations of motion become more complicated without it.
We follow the usual treatment of the electron-phonon coupling by considering small de- with respect to R n . The resulting electron-phonon hamiltonian in second quantization is where Ψ ≡ (Ψ ↑ , Ψ ↓ ) T denotes a spinor field operator. Note that ∂v e-ion /∂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 We follow the standard approach 22 and expand the field operators according to where are Bloch spinor wavefunctions with band and crystal momentum labels (µk). c µk is the corresponding destruction operator, and V is the normalization volume of the crystal. We use the convention that a Greek index includes band and pseudospin index, i.e., µ = (b, λ).
Instead of the Bloch wavefunctions, which are orthogonal on the whole crystal volume V, we will mainly deal with the lattice periodic u µ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â that occur in the present paper have the meaning By virtue of Eq. (8), the electron-phonon interaction hamiltonian (9) can be split into two terms the first one of which we call Elliott and the second Overhauser hamiltonian. The Elliott hamiltonian itself is spin diagonal whereas the Overhauser hamiltonian contains the electronic spin operatorŝ We next expand the phonon displacement operators according to 22 q,λ annihilates (creates) phonons with momentum q in the first Brillouin zone and mode index λ. We use the abbreviation of x which has the unit of length, where ̺ is the density. The phonon dispersion is denoted by ω q,λ . Further, ε q,λ is the polarization vector of the phonon mode with the property (ε q,λ ) * = −ε −q,λ . We introduce the Fourier transformation of the interaction potential in which has a long-range part (G = 0) and short-range part (G = 0).
Now we split the integral over the crystal volume V according to use the periodicity u µk (x + R (0) n ) = u µk (x) and the relation n e iR (0) n ·(q+k−k ′ ) = Nδ q,k−k ′ , which neglects Umklapp processes. We thus obtain for the e-pn coupling hamiltonian (9) with the matrix element in terms of the e-pn interaction operator v (λ) The "1" in the curly braces in the above expression comes from the Elliott term (15) and the rest from the Overhauser term (16). The q and k vectors are restricted to the first Brillouin zone, and we have already neglected Umklapp processes, which contribute if k + q lies outside the first BZ. The long-range part of this expression results from the G = 0 contribution to the sum, whereas the sum over the G = 0 defines the short-range part of the matrix elements. Grimaldi and Fulde 14 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. 20 Grimaldi and Fulde 14 also demonstrate how the longrange contribution is screened so that v eff (q) has a well defined q → 0 (long-wavelength) limit. In this limit, the long-range interaction operator iŝ v (λ) Note that for small q the electrons couple exclusively to longitudinal phonons, where the polarization vector ε q,λ points in the same direction as the wave-vector q.

B. Equation of motion for spin density-matrix
In this section we derive the equation of motion for the reduced spin density matrix 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 (10) in diagonal form Here, we use k · p theory, 23,24 i.e., the u µk 's and energy dispersions ǫ µk are determined as the solution of the eigenvalue problem of the self-adjoint k · p operatorĥ eff (k) for lattice-periodic u µk 's on Ω. This guarantees that, for each k, the eigenvalue spectrum ǫ µk is discrete and the u µk form a complete set, i.e., µ |u µk u µk | = 1.
Further, the u µk 's for each k are orthogonal with regard to the scalar product (13). In semiconductors, the k · p hamiltonian operator can be approximated on a subset of bands by an effective hamiltonian, such as the Luttinger hamiltonian. 24 The number of bands included in (25) determines the dimension of the reduced density matrix ρ k .
The equation of motion for the electronic spin-density matrix ρ µµ ′ k due to (20) and (25) is Equation (28) is the full equation of motion that contains the coherent contribution as well as the e-pn interaction contribution (29), in which the phonon-assisted electronic density matrix b q,λ c † νk+q c µ ′ k appears. 24,25 The Elliott-Yafet mechanism arises from the e-pn interaction contribution in (28). To study the EY mechanism in isolation, we will in the following neglect the coherent contributions to the equation of motion of the spin density-matrix.
When we specialize the results to the case of a pair of Kramers degenerate bands in Sec. IV below, the coherent contributions vanish exactly.

C. Dynamics of the average spin due to incoherent e-pn interaction
Before we make approximations to the phonon-assisted density matrix, we compute the dynamics of the average spin By combining (29) and (30) we obtain We can employ the completeness property (27) of the u µk 's to show that this has the form with the torque vector operator due to the e-pn interaction We will denote the matrix element of this operator by This matrix element completely determines the spin change that occurs in a transition (µk) → (µ ′ k + q). This information is not explicit the matrix element (21). Using (22), one The expression (32) together with (35) is already an important result, as it shows thatonly the Overhauser contribution to the e-pn interaction gives rise to a non-vanishing torque (35).
The Elliott contribution is spin diagonal and therefore its contribution to the torque vanishes.
We stress that this result is general enough to include materials with Kramers degeneracy (and "spin hot spots") as well as non-spin degenerate systems with avoided level crossings due to spin-orbit coupling. In all cases, the Elliott contribution to the torque matrix element does not change the average spin. This conclusion does not depend on the single-particle basis u µk in the sense that one could apply a k-dependent unitary transformation to the u µk 's and the c µk 's in (32) and (35). This argument shows that it is not the spin mixing in the wave functions that determines the Elliott torque, and-in this restricted sense-there is no Elliott contribution to spin dynamics in any basis.
The torque operatort k+q,k , which is due to the Overhauser contribution only, has the long-wavelength limit Finally, we note also that (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.

D. Dynamics of the average spin in the scattering limit
To compare the results of the present approach with the conventional EY analysis, we specialize the dynamical equation (32) for S α for the case of incoherent scattering. To this end, we need to evaluate the phonon-assisted density matrix b q,λ c † µk+q c µ ′ k . We do this by truncating its equation of motion at the scattering level to obtain where N q,λ = b( ω q,λ ) are phonon occupation numbers, and performing a Markov approximation, see, e.g., Refs. 24 and 25.
The equation of motion for the phonon-assisted density matrix and, consequently, for the reduced density-matrix ρ µµ ′ k is determined by the full matrix element including both the Elliott and Overhauser contributions, and one finds for the ensemble averaged spin We will assume that the phonon system is not changed appreciably by the electronic dynamics so that the N q,λ = b( ω q,λ ) are the equilibrium phonon occupations given by the Bose function b ( ω q,λ ). At this level of generality, we still have to deal with the reduced density matrix on the right-hand side.

IV. RELAXATION TIME FOR KRAMERS DEGENERATE BANDS
We now specialize to the case of a pair of Kramers degenerate bands with a quasiequilibrium spin-expectation value to obtain a relation for the spin relaxation, which generalizes Yafet's treatment.

A. Eigenstates for Kramers degenerate bands
Kramers degenerate states including spin-orbit coupling can be assumed to have the form 26-28 Part of this assumption is that these states are labeled by pseudo spin indices and that, Because of the degeneracy, there is an ambiguity to the definition of these states, as any superposition will also be an eigenstate. In accordance with Refs. 26-28 we choose them to fulfill This is equivalent to the statement that the chosen states u ⇑k , u ⇓k diagonalizeŝ z in the degenerate subspace. The choice is important, if one wants to attach some importance to the magnitude of 1 Ω´Ω |b k (x)| 2 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 (41) is fulfilled.
The condition (41) does not imply that the electrons are in pure spin-up or spin-down states characterized by the same two-dimensional spinors at each k. In fact, we have as well as | u λk |ŝ z u λk | ≤ /2 for λ =⇑, ⇓. For more details, see appendix A.

B. Quasi-equilibrium with spin polarization
In order to derive a characteristic rate for the relaxation of a small excess spin polarization δ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 δS z . Thus we must determine the quasi-equilibrium distribution for the eigenenergies of the generalized grand-canonical k · p hamiltonian Here, ζ z , the spin accumulation, is the Lagrange parameter needed to obtain the finite spin expectation value. Sinceĥ eff (k) includes spin-orbit coupling it does not commute withŝ z , so that the eigenstates and eigenenergies of are, in general, not identical to the u µk 's that result as eigenfunctions ofĥ eff (k), cf. Eq. (26).
We assume that ζ z is small and determine E µk perturbatively in the degenerate subspace of the u µk 's. Since we chose the u µk 's that diagonalizeŝ z in the degenerate subspace, cf. (41), this gives If we also assume that other bands are sufficiently far away, the eigenvectors are w µk ≃ u µ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 µ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 µk , i.e., where f (ǫ) = [exp(βǫ) + 1] −1 is the Fermi function. Since ζ z is small, we keep only the first order result in ζ z and with (45) we find for the quasi-equilibrium distribution where we have defined ∆ k ≡ − df dǫ | ǫ k −µ . For low temperatures, this function approaches a δ-function concentrated at ǫ k − µ. From the reduced density matrix (47) we find the relation ζ z = δS z /N between δS z and ζ z with the normalization

C. EY relaxation time
To obtain the close-to-equilibrium spin relaxation rate defined by we insert (47) in (38) and linearize. Then, as in Yafet's original treatment and its extension by Fähnle and coworkers, 13 ζ z drops out and we find µk,µ ′ k+q u µk+q |ŝ z u µk+q ∆ k+q This expression for the spin relaxation, or T 1 , 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 accounting for pseudo-spin flips. Instead of |g µk,µ ′ k ′ | 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 electronphonon interaction and the spin appear. As shown in Sec. III C the torque matrix element only has an Overhauser contribution, whereas in the g matrix element both Elliott and Overhauser terms contribute. Regarding the g matrix element in the long wavelength limit, Yafet 5 and numerical evaluations 28 of his spin-relaxation result find a cancellation between Elliott and Overhauser contributions in the short range contribution, whereas Grimaldi and Fulde 14 find for the long-wavelength limit a dominating long-range Elliott contribution to g.
Our result (50) is not limited to small spin mixing and is qualitatively different from the original Yafet formula, as there is no pure Elliott contribution. 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.

V. CONCLUSION
In conclusion, we have presented an investigation of the carrier-spin dynamics due to incoherent electron-phonon scattering in the presence of spin-orbit coupling, as originally considered by Overhauser, Elliott, and Yafet. We examined the dynamical equation for the reduced density matrix in terms of the phonon assisted density matrix including the contributions of the spin-diagonal (electrostatic) electron-phonon interaction as well as the phonon-modulated spin-orbit coupling. We showed that the central quantity to account for the spin change is a torque matrix element, which gives a physically more appealing picture of spin dynamics due to the EY mechanism. Our results show that the electronphonon "Elliott contribution," i.e., the spin mixing in the Bloch states, has no impact on the torque matrix element due to the electron-phonon interaction. This important property of If c k = 0, diagonalization leads to the eigenstates u λk with λ =⇑, ⇓ where L normalizes the u λk 's. The corresponding spin eigenvalues are In the presence of inversion symmetry, as in aluminium and silicon, these expressions simplify further. Then d k = d −k and the eigenvalues are ±( /2) d 2 k + |c k | 2 for λ =⇑, ⇓, respectively. This result shows that pseudo spin labels are justified as these states generally have an expectation value ofŝ z of modulus smaller than (or almost equal to) ± /2. At the same time, the pseudo-spin states are not just rotated pure spin states, which would lead to a spin projection of exactly ± /2 along the quantization axis. Further, these states may have nonvanishingŝ x andŝ y expectation values. * Graduate School of Excellence Materials Science in Mainz, Gottlieb-Daimler-Strasse 47, 67663