Microscopic linear response theory of spin relaxation and relativistic transport phenomena in graphene

We present a unified theoretical framework for the study of spin dynamics and relativistic transport phenomena in disordered two-dimensional Dirac systems with pseudospin-spin coupling. The formalism is applied to the paradigmatic case of graphene with uniform Bychkov-Rashba interaction and shown to capture spin relaxation processes and associated charge-to-spin interconversion phenomena in response to generic external perturbations, including spin density fluctuations and electric fields. A controlled diagrammatic evaluation of the generalized spin susceptibility in the diffusive regime of weak spin-orbit interaction allows us to show that the spin and momentum lifetimes satisfy the standard Dyakonov-Perel relation for both weak (Gaussian) and resonant (unitary) nonmagnetic disorder. Finally, we demonstrate that the spin relaxation rate can be derived in the zero-frequency limit by exploiting the SU(2) covariant conservation laws for the spin observables. Our results set the stage for a fully quantum-mechanical description of spin relaxation in both pristine graphene samples with weak spin-orbit fields and in graphene heterostructures with enhanced spin-orbital effects currently attracting much attention.


Introduction
impurities. For simplicity, in this work we neglect intervalley scattering processes and thus it suffices to consider the low-energy dynamics around the K point.
The energy dispersion relation of the free Hamiltonian H 0 = H − V in Eq. (1) is where µ, ν = ±1 labels the various subbands ( Fig. 1(a)). The Rashba interaction aligns the electron spin at right angles to the wavevector, the so-called spin-momentum locking configuration ( Fig. 1(b)) [44,45]. For | | > 2|λ| (region II), the split Fermi surface displays counter-rotating spin textures reminiscent of (non-chiral) 2DEGs with Rashba interaction [32]. A regime (pseudogap, region I) where the chemical potential intersects a single subband, with electronic states having well-defined spin helicity, extends for energies | | < 2|λ|. In the conventional 2DEG this circumstance only happens at a single point i.e., the intersection between the parabolic bands [46].
The spin texture of energy bands in the 2D Dirac-Rashba model is modulated by the band velocity i.e., where σ µνk = (1/v)∇ k µν (k) is the pseudospin polarization vector. The entanglement between pseudospin and spin degrees of freedom in the model is responsible for a rich energy dependence of transport coefficients [36,37]. For brevity of notation, we assume , λ > 0 in the remainder of the work.

Disorder effects
The random potential in Eq. (1) affects the spin dynamics by inducing elastic transitions between electronic states (µνk) → (µ ν k ) associated with different effective Larmor fields, Ω µνk = λ s µνk ≈ −µνλk ×ŝ for λ. This random change in the spin precession axis is responsible for the irreversible loss of spin information. To describe the effects of disorder, we employ standard many-body perturbation theory methods. We work within the zero-temperature Green's function formalism.
The retarded(R)/advanced(A) single-particle Green's function (a = A, R ≡ −, +) is where T is the time-ordering symbol and θ(.) is the Heaviside step function. Changing to frequency representation (ω ≡ ), one obtains where G a 0 ( ) = ( + ıvσ · ∇ + λ (σ × s) ·ẑ ± ı0 + ) −1 is the Green's function of free 2D Dirac-Rashba fermions and is the self-energy operator. The central quantity in our approach is the disorder averaged Green's function, G a (x − x , ) = G a (x, x ; ) , where ... denotes the average over all impurity configurations ( Fig. 2(a)). Its momentum representation is where is the disordered averaged self energy evaluated in the noncrossing approximation. The latter neglects corrections arising from quantum coherent multiple impurity scattering, which is justified in the diffusive regime with τ 1 [48,49]. The self-energy induced by short-range impurities is k-independent, Σ a k ( ) ≡ Σ a ( ), and hence we drop this index in what follows.
To account for the characteristic resonant (unitary) scattering regime of graphene with relaxation time τ ∝ [50, 51], we adopt a T -matrix approach by evaluating the self energy Σ a ( ) at all orders in V . We obtain where u 0 parameterizes the scattering strength of the spin-transparent (scalar) impurities, n i is the impurity areal density and T a ( ) is the single-impurity T -matrix. Note that multiple impurity scattering diagrams ∝ O(n 2 i ) can be neglected in the diffusive regime of interest. We have also introduced as the momentum integrated Green's function of the clean system [Eq.
The self energy simplifies in two important limiting cases: (i) weak Gaussian disorder (|u 0 | |g a 0 | −1 ) and (ii) unitary disorder (u 0 → ±∞). In the weak scattering regime, it suffices to only take into account the 'rainbow' diagram in the Dyson expansion Σ a ( ) = V +G a 0 ( )V +... ; see Fig. 2 (b). Note that this approximation is equivalent to assuming that the randomness of the disorder potential satisfies white-noise statistics In this case we have Σ a ( )| Gauss. = n i u 2 0 g a 0 ( ) .
The real part of the self-energy provides a parametrically small renormalization of the band structure, which can be safely neglected in the diffusive regime with 1/( τ ) 1 [36]. We thus find where the functions η 0 , η r , η KM , proportional to the imaginary parts of Eqs. (11)- (13), have different forms depending on the Fermi level position. In this work, we will restrict the analysis to diffusive systems with weak SOC λτ 1 and λ. It is thus convenient to express the various quantities in Σ a ( ) in terms of the quasiparticle broadening in regime II, i.e., 1 2τ Explicitly, we have Within the T -matrix approach, while Σ a ( ) can be again safely neglected, the nondiagonal part of Σ a ( ) acquires a finite value. However, in the unitary limit, we have Σ a ( ) = −n i /g a 0,0 ( ) and we recover a scalar self-energy, with The unitary result captures the typical energy dependence τ ∝ observed in high-mobility graphene samples [50], where the charge carrier mobility is likely limited by short-range scatterers, including adsorbates, short-range ripples and vacancies [52,53,54,55].

General formalism
We consider the long-wavelength spin dynamics generated by a generic external perturbation where J αβ ∝ σ α s β (α, β = 0, i) is the current density operator (α = x, y) or density operator (α = 0, z) and A αβ is a generalized vector potential [36]. We will consider in detail two important cases: (i) an electric field perturbation e.g., The induced spin polarization density is evaluated within the framework of linear response theory. This approach has been applied to derive charge-spin diffusion equations describing spin dynamics and magnetoelectric effects in 2DEGs [41,56,57]. As shown below, a suitable extension of this approach to accommodate the enlarged (spin ⊗ pseudospin) Clifford algebra γ αβ = σ α s β will allow us to obtain a rigorous microscopic theory of diffusive transport and spin relaxation for 2D Dirac systems. The linear response of theî-component of the spin polarization vector at zero temperature reads as where is the generalized spin susceptibility associated to the external perturbation i.e., an electric [58]. Expressing the above equation in terms of the Fourier transform χ i,αβ (q, ω) in the long-wavelength limit q → 0 we have where κ = v (κ = 1/2) for a electric (spin injection) field and Tr is the trace over all degrees of freedom. Terms involving products of Green's functions on the same sector (RR and AA) are smaller by a factor of ( τ ) −1 and thus are neglected. The disorder average in Eq. (25) is evaluated by means of the diagrammatic technique (Fig. 3). For brevity of notation, we first present the formalism within the Gaussian approximation for the self-energy, Eq. (17). In Sec. 2.3, we provide the connection with the full T -matrix result.
A summation of non-crossing two-particle (ladder) diagrams leads to where tr is the trace over internal degrees of freedom (spin and sublattice). The dressed vertexγ αβ satisfies the Bethe-Salpeter (BS) equationγ where N 0 ≡ /πv 2 (for the T -matrix extension see Eq. (61) and text therein). Projecting onto the elements of the Clifford algebraγ we recast the BS equation into the form where Introducing the 16-dimensional vectorsγ αβ (ω) = (γ αβ00 (ω), ...,γ αβzz (ω)) t and γ αβ = (0, 0, ...., γ αβαβ , ..., 0) t a more compact matrix form for Eq. (27) is given in terms of the diffuson operator D −1 as The spin relaxation rates are simply identified as the poles of the generalized susceptibility in the complex ω-plane. The determination of the SRTs is thus reduced to the analysis of the behavior of D −1 = D −1 (ω) [59]. The formal result Eq. (31) deserves a few comments. Firstly, D −1 spans in principle the entire Clifford algebra, which physically encodes the coupled dynamics of spin and other observables associated with the elements γ αβ . However, by exploiting symmetries, D −1 can be reduced into block diagonal form, such that only some observables are coupled to the spin polarizations along the three spatial directions. Secondly, a distinct feature of Dirac systems is that spin densities are coupled to charge currents even in the case (considered here) of spatially uniform external perturbations q=0. The linear Dirac dispersion of graphene is reflected in the form of the charge current J i = vσ i and spin current J a i = vσ i s a /2 vertices, which do not depend explicitly on momentum; by virtue of that they can be directly identified (apart from constants) as elements of the Clifford algebra. Therefore all the relevant information about coupling between currents and densities is built-in on the 16×16 diffusion operator Eq. (31) in our formalism. This will allows us to obtain a unified description of spin relaxation processes and relativistic transport phenomena (e.g., charge-to-spin conversion) within our q=0 formalism. We analyze the implications below.
The coupling of the electrons' spin to currents or other observables in the long wavelength limit also suggests two equivalent scenarios to study spin relaxation. The first natural choice is to consider spin injection and investigate the relaxation of the spin density profile (density-density response); alternatively, one can probe the spin response indirectly by exciting an observable coupled to the spin density through D −1 . For instance, as we will see in the following, one can drive a charge current via application of an electric field to obtain a in-plane spin polarization of carriers (ISGE). In that case, the information about the in-plane SRTs is readily accessible by examining how the steady state (Edelstein) polarization is achieved (density-current response).
Before moving on, let us stress that within the Gaussian approximation, a useful relation can be derived connecting the generalized susceptibility Eq. (26) and the renormalized vertex: where α ≡ (2πτ N 0 ) −1 and we have used Eq. (29). The above equation states the spin response is related solely to the associated component of the renormalized vertexγ αβ0i . A similar relation holds for other response functions. For example the AC longitudinal (Drude) conductivity is written as Therefore Eq. (32) and similar relations allow to identify the components of a renormalized vertex with the associated observables, and will turn useful in the following. Let us now determine the allowed couplings to S x,y,z by exploring symmetry. The model of Eq. (1) is invariant under the group C ∞v , which is an emergent symmetry of the continuum (long-wavelength) theory. As rotations in the continuum do not describe the sublattice symmetry A ↔ B of the graphene system, a representation U for the relevant set of discrete operations has to be considered. Relevant to us are C 2 , the rotation of π around theẑ-axis exchanging sublattice (and valleys), and R x , the reflection over thex-axis. We have with r x : (x, y) → (x, −y) and τ i=x,y,z are Pauli matrices acting on the valley degree of freedom. We also make use of isospin (valley) rotations Λ x,y,z [60,61] For scalar disorder it suffices to examine the form of the clean-system susceptibility at ω = 0 For any of the symmetries S listed above, we have , and inserting resolutions of the identity in the form S † S into Eq. (38) we find where p αβ (p 0i ) = ±1 is the parity of γ αβ (γ 0i ) under S. From this result, we see that a non-zero response requires the operator γ αβ to have the same parity of the spin vertex under the action of any of S. The allowed couplings and parities under S are shown in the Tab. 1. As anticipated above, the in-plane components S x(y) are coupled to orthogonal charge currents σ y(x) , as well as spin Hall currents γ xz (γ yz ) and staggered magnetizations γ zy (γ zx ) [36,37]. The out-of-plane component S z is instead coupled to a mass term σ z and in-plane spin currents γ xx , γ yy .

Diffusive equations and SRTs
In the following, we choose to consider the in-plane spin response to an AC electric field y. This choice, as discussed above, is equivalent to consider in-plane spin injection [62], but has the advantage to allow for a unified description of spin dynamics and charge-spin interconversion, e.g. to capture spin-Galvanic (Edelstein) effect. For the out-of-plane spin dynamics, we take a spin-density perturbation H ext ⊥ = 1 2 s z B z (ω) (see Tab. 1).

In-plane spin dynamics
Without loss of generality, let us consider the dynamics of theŷ polarization. According to Tab. 1, s y is coupled to three operators: σ x , σ y s z and σ z s x . However, leading terms in the ( τ ) −1 expansion are only contained in the s y /σ x sub-block. Hence, to capture the SRTs it suffices to restrict to this 2 × 2 algebra. As anticipated above, we consider here the response to an AC electric field E x (ω), associated with the vertex κγ x0 = vσ x ≡ v x . (Details of calculation and full form of the 4 × 4 diffuson operator is given in Appendices C and D.) To capture purely diffusive processes, we expand D −1 (ω) in the low-frequency and small SOC limits, ωτ 1 and λτ 1, respectively. In this regime, Eq. (31) is written then as where Γ s = τ /τ s = 2λ 2 τ . In the light of previous discussions (cf. Eqs. (32) and (33)),ṽ x0 andṽ 0y are connected by a linear transformation to the steady-state charge current and the spin polarization (Appendix D).
Off-diagonal elements of Eq. (40) carry in relation to diagonal ones an extra order of smallness λ/ , suggesting spin and charge to be weakly coupled in this limit. Their inclusion however encodes charge-to-spin interconversion and it is essential to get a correct physical description. The eigenvalues −ıω ± are found as and can be associated with charge current and spin relaxation times, respectively. We see then the SRT can be identified as the mass (ω = 0) term of the spin-spin part of the diffuson Inverting Eq. (40), we findγ from which, by using Eqs. (32) and (33) is it possible, upon Fourier transform, to derive the diffusive equation of motion for coupled charge-spin dynamics as where J 0 x (t) ≡ 2v 2 E x (t)/α and S 0 y (t) ≡ −λE x (t)/ α. Note charge current relaxation is regulated by the transport time τ tr ≡ 2τ , indicating the absence of backscattering [36,49,50].

Out-of-plane spin dynamics
For the out-of-plane spin dynamics we consider the renormalized vertex κγ 0z = 1 2s z . The off diagonal components of the associated 4 × 4 diffuson block contains sub-leading terms in the ( τ ) −1 expansion (Appendix C), such that the out-of-plane SRTs can be calculated similarly to Eq. (43) as The generalization of the equations of motion Eqs. (46), (47) in this case is written as where S 0 z (t) =Ḃ z (t)/4α is the effect of the external perturbation (spin-injection field). The in-plane and out-ofplane SRTs are in the following relation which is nothing but the well-known DP ratio for 2DEGs [56]. The above result has also been obtained for graphene within the time-dependent perturbation theory for the density matrix [35]. The agreement between graphene and the Rashba 2DEG results is expected at high electronic density λ.

SRT from the conservation laws in the DC limit
In this section, we demonstrate how the SRTs we have obtained above can be equivalently extracted in the static limit ω = 0. This remarkable result is rooted in the conservation laws associated to the disordered Dirac-Rashba Hamiltonian Eq. (1) [36]. The first step is to write the Heisenberg equation of motion for the spin polarizations where lj , c li are the second and third rank Levi-Civita tensors and J c j = J c j is theĵ-component of the pure spin current (with polarization "c"). As before, we consider an electric field applied along thex direction. We find where J z y is identified as the spin Hall current according to the chosen geometry. The spin Hall current is written in response to the electric field where σ z yx is the DC spin Hall conductivity. As for now no assumption has been made for the self-energy approximation associated to the scalar impurities field. Let us start from the more transparent Gaussian case. Using the corresponding version of Eq. (32) for σ z yx , together with Eq. (29) we have In the latter we have neglected the terms M yzyz and M zxyz which, as said above, provide higher order corrections in the ( τ ) −1 expansion. Multiplying both sides of Eq. (54) by the electric field E x , and using S y = χ y,x0 E x together with Eq. (32), we find Analogously to the 2DEG case, in the steady state ω = 0 the latter is identically zero, despite the Dirac character of fermions [36]. The absence of spin Hall current as imposed by the steady-state version of the continuity equation Eq.(52), imposes the out-of-equilibrium value Evaluating the above quantities explicitlyγ x0x0 = 2, M x0yz /M 0yyz = λ/ and we recover the ISGE obtained in Ref. [37]. Using Eq. (52) we finally arrive at and therefore where we have identified the spin relaxation time in perfect accordance with the result obtained above, Eq. (43). The bubble M 0yyz is therefore what completely determines the in-plane spin relaxation. We now ask how the above result is modified when treating the self-energy in the T -matrix approximation. The Bethe Salpeter equation Eq. (27) now reads where T R/A ( ) is the single-impurity T -matrix in the R/A sectors introduced in Eq. (9). Projecting onto the Clifford algebra, similarly to Eq. (29), we havẽ where we have defined Recasting Eq. (62) in vector notation, in the same spirit of Eq. (31), we havẽ and consequently The latter equation allows again to find a connection with the observables. For example, the generalization of Eq. (32) is written as The spin Hall conductivity instead is found as Differently to the Gaussian case, where we were able to relate the response of an observable uniquely to the associated component of the renormalized vertex, in the T -matrix limit in principle all components ofγ x0 would contribute, each of them with weight given by Y −1 . In the limiting case of unitary limit u 0 → ∞, where lim u0→∞ , we find a simplification as This implies that for Eq.(68) a relation similar to the Gaussian case is obtained where we have restricted ourselves again to the dominant subspace σ x /s y . After standard algebra, we arrive at and the SRT defined as where we have used the definition of τ in the unitary limit, Eq. (21). We conclude that the the formal expression connecting τ s and τ (the DP relation) is the same as found in the Gaussian limit for the self-energy. However, given the different dependence of τ on the Fermi level in the two approximations-cf. Eq. (19) and Eq. (21)-one has Unitary. (73) The SRT associated to the out-of-plane component can be derived along the same lines. The relevant Heisenberg equation now reads and a similar reasoning that lead to Eq, (60), allows us to conclude in the Gaussian limit, and a similar relation for the unitary limit.

Conclusions
In the present work, we laid the foundations of a general microscopic theory of diffusive transport and spin relaxation in 2D Dirac systems subject to spin-orbit interactions. Our work represents the logical extension of the previouslydeveloped diagrammatic treatments [56,57] to all orders in the scattering potential, for disordered electron systems with an enlarged pseudospin ⊗ spin Clifford algebra [36,37]. We applied the formalism to the paradigmatic case of 2D Dirac fermions with Rashba spin-orbit coupling considering the purely diffusive regime λτ 1 τ . We demonstrated how the Dyakonov-Perel relation between momentum and spin lifetime τ ∝ τ −1 s holds in both the Gaussian (weak short-range scatterers) and the unitary (strong short-range scatterers) regimes, despite the drastic different dependence momentum scattering times τ = τ ( ) in the two regimes. We derived the same result both by direct diagrammatic resummation (in the non-crossing approximation) and by exploiting the conservation laws of the theory in the zero-frequency limit. Under the diffusive regime λτ 1 τ is not possible to study the dynamics in the region of Fermi energies comparable to the Rashba pseudogap region ∼ 2λ, which was recently predicted to display interesting out-of-equilibrium phenomena [37]. The strong spin-momentum locking approaching this regime lets us infer a modification of the relation between τ s and τ towards the Elliot-Yafet type τ s ∝ τ . Our theory sets the stage to the study the spin dynamics in that regime. This topic has become of renewed great interest due to recent progresses in graphene-based heterostructures, where the spin relaxation anisotropy has been recognized as a viable tool to estimate the induced large spin-orbital effects [65,66,67,68].

A Single-Particle Green's Function: clean system
The explicit form of the clean single-particle Green's function is where L and φ k is the wavevector.

B Integrals and expansion
The current work makes extensive use of momentum integrations involving products of two renormalized Green's functions with analyticity in opposite halves of the complex plane, see e.g. Eqs. (25), (27). The retarded function is displaced in energy by the amount ω. Similarly to the bare Green's function decomposition Eqs. (76) and (78), we write the renormalized (disorder averaged) propagators as where M a ik ( ) = M are matrix coefficients and the kernels L a µ = L a ik are obtained in the Gaussian limit from the functions L 0µ of Eq. (77) by analytical continuation → + a ı 2τ . In the T matrix approach, the analytical continuation has to be performed as to include the other matrix structure of the self-energy ∝ γ r , γ KM [37]. We can generically write where z a i ( ) are complex quantities. Given the decomposition in Eq. (79), the integrals we need to solve are reduced to product of two kernels in different combinations, accompanied or not by a factor v 2 k 2 . Terms proportional to v 4 k 4 can be shown to vanish upon angular integration dφ k . We write below an exact solution and then expand at linear order in ω. For simplicity we show the results for the Gaussian approximation. The first type of integrals is where the principal branch of the log function has been chosen. Note z a 1 (λ) → z a 2 (−λ). Thus, Γ 11 (λ) = Γ 22 (−λ) and Γ 12 (λ) = Γ 21 (−λ). At linear order in ω we find where we have retained leading order terms in ( τ ) −1 . The other class of integrals we need to solve is where the ultraviolet cutoff Λ/v k F has been introduced to regularize the integrals. A careful evaluation yields Γ and the expressions for 1 ↔ 2 again obtainable with the replacement λ → −λ.

C Full form of the diffuson
Here we report the full form of the two relevant blocks of the diffuson, involving S y,z . The expressions are provided at leading order in the expansions for ωτ λτ 1 τ .
• Subspace σ x , s y , σ y s z , σ z s x • Subspace s z , σ x s x , σ y s y , σ z

D Equation for operators instead of vertices
In the main text, we have written equations of motion for the renormalized vertices, rather than for the observables themselves. As an example, here we report the diffusive matrix D −1 for the observables, in the relevant sub-block s y /σ x for the in-plane spin dynamics. Also here we consider the response to an external electric field E x . To this aim we recall in the Gaussian approximation (cf. Eq. (32) and (33)) Manipulating Eq. (31) we have where we have defined the matrix Consequently by subtracting to both sides v 2 We conclude the diffusive matrix for the observables is As anticipated D −1 obs has the same eigenvalues as D −1 , justifying the approach in the main text.