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

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

Manuel Offidani , Roberto Raimondi  and Aires Ferreira 

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, allowing to capture spin relaxation processes and associated charge-to-spin inter- conversion phenomena in response to generic external perturbations, including spin density fluctuations and electric fields. By means of a self-consistent diagrammatic treatment of the generalized spin susceptibility, where the impurity potential is taken at all orders in perturbation theory, we recover the Dyakonov-Perel relation for the spin relaxation times (SRTs), which is shown to hold both for weak scatterers (Gaussian approximation) and in the unitary limit of resonant impurities. Finally, we show that the SRTs 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.

Department of Physics, University of York, York YO10 5DD, United Kingdom

Dipartimento di Matematica e Fisica, Università Roma Tre, 00146 Rome, Italy

1 Introduction

1.1 Spin relaxation in graphene

Graphene is considered a promising material for spintronics applications due to its negligible hyperfine interactions and low spin–orbit coupling (SOC) [1, 2]. Early theoretical estimates hinted at ultra-long spin lifetime (1–100 s) [3], whereas experiments found to be limited to a few nanoseconds [4]. The microscopic mechanisms responsible for the relatively fast spin relaxation in high-mobility graphene samples remain controversial [5], but recent findings indicate that spinful scatterers, such as magnetic adatoms, are the primary cause of spin relaxation [6, 7, 8, 9].

The spin dynamics in graphene is conventionally probed by means of non-local transport measurements [10, 11]. In this approach, a spin current is injected from ferromagnetic electrodes into graphene and allowed to diffuse under the effect of a perpendicular magnetic field. The Larmor precession of the electron’s spin about the external field modulates the average spin accumulation detected away from the injection point (Hanle curve), a bona fide spin signal from which can be deduced. Hanle precession measurements found a large spread in from tens of picoseconds up to a few nanoseconds [12, 13, 14, 15, 16, 17, 18, 19, 20], reflecting the different sample preparation and device fabrication methods. Theoretical studies have revealed a number of possible spin relaxation sources, including magnetic impurities, spin–orbit active adatoms, ripples and other substrate effects [21, 22, 23, 24, 25, 26, 27, 28]. Despite the relatively short , the high mobility of clean graphene samples allows spins to diffuse over extremely long distances up to 13 m at room temperature [29, 30, 31].

The paradigmatic model for studies of spin relaxation in graphene is the two-dimensional (2D) Hamiltonian of massless Dirac fermions supplemented with a (uniform or random) Bychkov-Rashba interaction [32]. This type of SOC has its origin in perturbations breaking the inversion symmetry, which include substrate-induced electric fields, adatoms, and ripple-induced gauge fields [3, 4]. The Bychkov-Rashba interaction in graphene (hereafter referred to as Rashba SOC for brevity) can be seen as a non-Abelian gauge field that couples to the intrinsic pseudospin of Dirac fermions, enabling spin relaxation upon impurity scattering e.g., via the familiar Dyakanov-Perel (DP) mechanism [33]. Previous theoretical descriptions of spin relaxation in such 2D Dirac-Rashba models were based on semiclassical approximations [34, 35]. On the other hand, a fully quantum-mechanical theory of spin–orbit-coupled transport in the static (DC) limit has been formulated recently by the authors [36, 37]. Analogously to the 2D electron gas (2DEG) case [38, 39, 40], it was shown that impurity scattering corrections exactly balance the intrinsic generation of spin Hall current for nonmagnetic disorder, , where and is an external DC electric field [36]. The vanishing of the spin Hall effect in this model is connected to the establishment of a robust nonequilibrium in-plane spin polarization [inverse spin-Galvanic effect (ISGE)] [37]. However, a time-dependent framework able to unveil how the steady state is reached within the 2D Dirac-Rashba model is yet to be developed. In this paper, we address this problem. We derive the coupled spin-charge drift-diffusion equations for nonmagnetic disorder and generic homogeneous perturbations by means of the diagrammatic technique for disordered electrons. A similar approach has been adopted very recently in the context of 2DEGs with both Bychkov-Rashba and Dresselhaus interactions [41], where it was shown perfect agreement between the Kubo diagrammatic formalism and the Keldysh SU(2) gauge theory [42]. In this work, we extend the standard quantum diagrammatic formalism to accommodate the enlarged (spin) 2 (pseudospin) Clifford structure of the 2D Dirac–Rashba model, for which the diffuson operator is 16-dimensional. We find that the typical DP relation connecting the spin relaxation time (SRT) and the momentum lifetime in the weak SOC regime, that is for , where is the SOC strength, holds at all orders in the scattering potential strength. The meaning and interpretation of our results for the SRTs can be also clarified by the SU(2) covariant conservation laws inherent to the diagrammatic (perturbative) structure, whose usage allows us derive the DP relation even in the zero-frequency limit. In particular, we provide the analytical expression of in the unitary limit of very strong scattering.

1.2 Dirac–Rashba model

The effective low-energy Hamiltonian describing the electronic properties of 2D Dirac fermions subject to a uniform Rashba interaction around the point reads as [43]


where is the bare velocity of massless Dirac fermions, is the 2D kinematic momentum operator, is the SOC strength and () are Pauli matrices associated with sublattice (pseudospin) and spin degrees of freedom, respectively. Here, is a disorder potential describing elastic scattering from nonmagnetic short-range impurities. For simplicity, in this work we neglect intervalley scattering processes and thus it suffices to consider the low-energy dynamics around the point.

The energy dispersion relation of the free Hamiltonian in Eq. (1.2) is


where labels the various subbands (Fig. 1(a)).

Fig. 1: (a) Energy dispersion around the point. The splitting of the Dirac bands leads to a spin gap or pseudogap. (b) Tangential winding of the spin texture in regimes I and II.

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 (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 . 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 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 in the remainder of the work.

1.3 Disorder effects

The random potential in Eq. (1.2) affects the spin dynamics by inducing elastic transitions between electronic states associated with different effective Larmor fields, 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 () is


where is the time-ordering symbol and is the Heaviside step function. Changing to frequency representation (), one obtains


where 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, , where denotes the average over all impurity configurations (Fig. 2(a)). Its momentum representation is


where is the Fourier transform of (see Appendix A) and


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 [48, 49]. The self-energy induced by short-range impurities is -independent, , and hence we drop this index in what follows.

Fig. 2: (a) Dyson equation for the disordered averaged Green’s function. (b-c) Approximation schemes for evaluation of the self energy: Gaussian (b) and -matrix approximation (TMA) (c). Box shows Feynman rules for the disorder potential insertions (dashed lines) and impurity density insertion (red crosses).

To account for the characteristic resonant (unitary) scattering regime of graphene with relaxation time [50, 51], we adopt a -matrix approach by evaluating the self energy at all orders in . We obtain


where parameterizes the scattering strength of the spin-transparent (scalar) impurities, is the impurity areal density and is the single-impurity -matrix. Note that multiple impurity scattering diagrams 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. (A)], where is the identity matrix, , and


In the above, selects the energy regime and


with denoting the ultraviolet cutoff of the low-energy theory [50].

The self energy simplifies in two important limiting cases: (i) weak Gaussian disorder () and (ii) unitary disorder (). In the weak scattering regime, it suffices to only take into account the ’rainbow’ diagram in the Dyson expansion ; 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


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 [36]. We thus find


where the functions , proportional to the imaginary parts of Eqs. (1.3)-(1.3), have different forms depending on the Fermi level position. In this work, we will restrict the analysis to diffusive systems with weak SOC and . It is thus convenient to express the various quantities in in terms of the quasiparticle broadening in regime II, i.e.,


Explicitly, we have


Within the -matrix approach, while can be again safely neglected, the nondiagonal part of acquires a finite value. However, in the unitary limit, we have 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].

2 Microscopic linear response theory for spin relaxation

2.1 General formalism

We consider the long-wavelength spin dynamics generated by a generic external perturbation


where () is the current density operator () or density operator () and is a generalized vector potential [36]. We will consider in detail two important cases: (i) an electric field perturbation e.g., and (ii) a spin density fluctuation . 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 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 field or a ’spin injection field’ [58]. Expressing the above equation in terms of the Fourier transform in the long-wavelength limit we have


where () 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 and thus are neglected.

Fig. 3: Diagrammatic technique for evaluation of generalized spin susceptibilities. (a) Two-particle ladder diagram. (b) BS equation for the vertex renormalization. (c) Skeleton expansion of the ladder diagram in terms of an infinite series of two-particle, noncrossing diagrams. Full (open) square denotes a () matrix insertion.

The disorder average in Eq. (2.1) 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. (1.3). In Sec. 2.3, we provide the connection with the full -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 (for the -matrix extension see Eq. (2.3) and text therein). Projecting onto the elements of the Clifford algebra


we recast the BS equation into the form




Introducing the 16-dimensional vectors and a more compact matrix form for Eq. (2.1) is given in terms of the diffuson operator 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 [59].

The formal result Eq. (2.1) deserves a few comments. Firstly, 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, 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 and spin current 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 diffusion operator Eq. (2.1) 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 . 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. (2.1) and the renormalized vertex:


where and we have used Eq. (2.1). The above equation states the spin response is related solely to the associated component of the renormalized vertex . A similar relation holds for other response functions. For example the AC longitudinal (Drude) conductivity is written as


Therefore Eq. (2.1) 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 by exploring symmetry. The model of Eq. (1.2) is invariant under the group , which is an emergent symmetry of the continuum (long-wavelength) theory. As rotations in the continuum do not describe the sublattice symmetry of the graphene system, a representation for the relevant set of discrete operations has to be considered. Relevant to us are , the rotation of around the -axis exchanging sublattice (and valleys), and , the reflection over the -axis. We have


with and are Pauli matrices acting on the valley degree of freedom. We also make use of isospin (valley) rotations [60, 61]


For scalar disorder it suffices to examine the form of the clean-system susceptibility at


For any of the symmetries listed above, we have , and inserting resolutions of the identity in the form into Eq. (2.1) we find


where is the parity of under . 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 . The allowed couplings and parities under are shown in the Tab. 1. As anticipated above, the in-plane components are coupled to orthogonal charge currents , as well as spin Hall currents and staggered magnetizations [36, 37]. The out-of-plane component is instead coupled to a mass term and in-plane spin currents .

Polarization Couplings
-1 -1 +1
-1 +1 +1
+1 -1 +1
Tab. 1: Table summarizing the couplings to the spin polarizations along the three spatial directions in the 2D Dirac–Rashba model with nonmagnetic scalar disorder.

2.2 Diffusive equations and SRTs

In the following, we choose to consider the in-plane spin response to an AC electric field , . 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 (see Tab. 1).

2.2.1 In-plane spin dynamics

Without loss of generality, let us consider the dynamics of the polarization. According to Tab. 1, is coupled to three operators: and . However, leading terms in the expansion are only contained in the sub-block. Hence, to capture the SRTs it suffices to restrict to this algebra. As anticipated above, we consider here the response to an AC electric field , associated with the vertex . (Details of calculation and full form of the diffuson operator is given in Appendices C and D.) To capture purely diffusive processes, we expand in the low-frequency and small SOC limits, and , respectively. In this regime, Eq. (2.1) is written then as


where . In the light of previous discussions (cf. Eqs. (2.1) and (2.1)), and are connected by a linear transformation to the steady-state charge current and the spin polarization (Appendix D).

Off-diagonal elements of Eq. (2.2.1) 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 () term of the spin-spin part of the diffuson


Inverting Eq. (2.2.1), we find


from which, by using Eqs. (2.1) and (2.1) is it possible, upon Fourier transform, to derive the diffusive equation of motion for coupled charge-spin dynamics as


where and . Note charge current relaxation is regulated by the transport time , indicating the absence of backscattering [36, 49, 50].

2.2.2 Out-of-plane spin dynamics

For the out-of-plane spin dynamics we consider the renormalized vertex . The off diagonal components of the associated diffuson block contains sub-leading terms in the expansion (Appendix C), such that the out-of-plane SRTs can be calculated similarly to Eq. (2.2.1) as


The generalization of the equations of motion Eqs. (2.2.1),(2.2.1) in this case is written as


where is the effect of the external perturbation (spin-injection field). The in-plane and out-of-plane 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 .

2.3 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 . This remarkable result is rooted in the conservation laws associated to the disordered Dirac–Rashba Hamiltonian Eq. (1.2) [36]. The first step is to write the Heisenberg equation of motion for the spin polarizations


where are the second and third rank Levi-Civita tensors and is the -component of the pure spin current (with polarization ). As before, we consider an electric field applied along the direction. We find


where 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 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. (2.1) for , together with Eq. (2.1) we have


In the latter we have neglected the terms and which, as said above, provide higher order corrections in the expansion. Multiplying both sides of Eq. (2.3) by the electric field , and using together with Eq. (2.1), we find


Analogously to the 2DEG case, in the steady state 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.(2.3), imposes the out-of-equilibrium value


Evaluating the above quantities explicitly and we recover the ISGE obtained in Ref. [37]. Using Eq. (2.3) we finally arrive at


and therefore


where we have identified the spin relaxation time


in perfect accordance with the result obtained above, Eq. (2.2.1). The bubble 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 -matrix approximation. The Bethe Salpeter equation Eq. (2.1) now reads


where is the single-impurity -matrix in the R/A sectors introduced in Eq. (1.3). Projecting onto the Clifford algebra, similarly to Eq. (2.1), we have


where we have defined


Recasting Eq. (2.3) in vector notation, in the same spirit of Eq. (2.1), we have


and consequently


The latter equation allows again to find a connection with the observables. For example, the generalization of Eq. (2.1) 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 -matrix limit in principle all components of would contribute, each of them with weight given by . In the limiting case of unitary limit , where we find a simplification as


This implies that for Eq.(2.3) a relation similar to the Gaussian case is obtained


where we have restricted ourselves again to the dominant subspace . After standard algebra, we arrive at


and the SRT defined as


where we have used the definition of in the unitary limit, Eq. (1.3). We conclude that the the formal expression connecting 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. (1.3) and Eq. (1.3)—one has


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, (2.3), allows us to conclude


in the Gaussian limit, and a similar relation for the unitary limit.

3 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 previously-developed 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 . We demonstrated how the Dyakonov-Perel relation between momentum and spin lifetime 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 is not possible to study the dynamics in the region of Fermi energies comparable to the Rashba pseudogap region , 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 and towards the Elliot-Yafet type . 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




and 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. (2.1),(2.1). The retarded function is displaced in energy by the amount . Similarly to the bare Green’s function decomposition Eqs. (A) and (A), we write the renormalized (disorder averaged) propagators as