Pseudospin-driven spin relaxation mechanism in graphene

The possibility of transporting spin information over long distances in graphene, owing to its small intrinsic spin-orbit coupling (SOC) and the absence of hyperfine interaction, has led to intense research into spintronic applications. However, measured spin relaxation times are orders of magnitude smaller than initially predicted, while the main physical process for spin dephasing and its charge-density and disorder dependences remain unconvincingly described by conventional mechanisms. Here, we unravel a spin relaxation mechanism for nonmagnetic samples that follows from an entanglement between spin and pseudospin driven by random SOC, which makes it unique to graphene. The mixing between spin and pseudospin-related Berry's phases results in fast spin dephasing even when approaching the ballistic limit, with increasing relaxation times away from the Dirac point, as observed experimentally. The SOC can be caused by adatoms, ripples or even the substrate, suggesting novel spin manipulation strategies based on the pseudospin degree of freedom.

The possibility of transporting spin information over long distances in graphene, owing to its small intrinsic spin-orbit coupling (SOC) and the absence of hyperfine interaction, has led to intense research into spintronic applications. However, measured spin relaxation times are orders of magnitude smaller than initially predicted, while the main physical process for spin dephasing and its charge-density and disorder dependences remain unconvincingly described by conventional mechanisms. Here, we unravel a spin relaxation mechanism for nonmagnetic samples that follows from an entanglement between spin and pseudospin driven by random SOC, which makes it unique to graphene. The mixing between spin and pseudospin-related Berry's phases results in fast spin dephasing even when approaching the ballistic limit, with increasing relaxation times away from the Dirac point, as observed experimentally. The SOC can be caused by adatoms, ripples or even the substrate, suggesting novel spin manipulation strategies based on the pseudospin degree of freedom.
The electronic properties of monolayer graphene strongly differ from those of two-dimensional metals and semiconductors in part because of inherent electron-hole band structure symmetry and a particular density of states which vanishes at the Dirac point [1]. Additionally, the sublattice degeneracy and honeycomb symmetry lead to eigenstates that hold an additional quantum (Berry's) phase, associated with the so-called pseudospin quantum degree of freedom. All of these electronic features are manifested through the Klein tunneling phenomenon [2], weak antilocalization [3] or the anomalous quantum Hall effect [4]. The possibility of using the pseudospin as a means to transport and store information has also been theoretically proposed [5,6]. There, the role of the pseudospin is equivalent to that of the spin in spintronics, such as in the pseudospin analogue of the giant magnetoresistance in bilayer graphene [6].
Even though pseudospin-related effects drive most of the unique transport signatures of graphene, the role of the pseudospin on the spin relaxation mechanism has not been explicitly addressed and quantified.
Pseudospin and spin dynamics are usually perceived as decoupled from each other, with pseudospin lifetimes being much shorter and pseudospin dynamics much faster than those for spins. However, this picture breaks down in the vicinity of the Dirac point, a region that is usually out of reach of perturbative approaches and that is particularly relevant for experiments, because Fermi energies can only be shifted by about 0.3 eV via electrostatic gating. Moreover, in the presence of SOC, spin couples to orbital motion, and therefore to pseudospin [7], so that spin and pseudospin dynamics should not be treated independently. Actually, any initially spin polarized state injected at low energy (for instance |Ψ ⊗ | ↑ ), should evolve under the time-evolution operator (including the spin-orbit coupling term) towards an introducing the pseudospin as a column vector), exhibiting a spin-pseudospin locking effect (see discussion associated to Figs. S1 and S2 in the Supplementary Information), and therefore entangled spin and pseudospin dynamics.
The reason for overlooking the role of the pseudospin on the spin dynamics is perhaps rooted in the fact that the spin transport properties appear remarkably similar to those found in common metals and semiconductors [8]. Indeed, spin precession measurements in nonlocal devices result in experimental signatures that would be indistinguishable from those obtained in a metal such as aluminium [9], or a semiconductor such as GaAs [10], with extracted spin relaxation times τ s that are also typically of the same order of magnitude (a few nanoseconds or lower). Spin relaxation in graphene has therefore been interpreted using the conventional experimental manifestations of either the Elliot-Yafet (EY) or Dyakonov-Perel (DP) mechanism [11][12][13][14][15][16]. In the EY scenario, the spin relaxation time is determined by the spin mixing of carriers and the SOC of the scattering potential, and thus it is usually assumed to be proportional to the momentum relaxation time as τ s ≈ α · τ p , with α ≫ 1 (in alkali metals α ∼ 10 4 − 10 6 ) [8]. In contrast, in the DP mechanism spin precesses about an effective magnetic field whose orientation is fixed by the momentum direction during free propagation of electrons. Such orientation changes at each scattering event, which results in a different scaling behavior as 1/τ DP s ∼ Ω 2 τ p [8] (with Ω the average magnitude of the intrinsic Larmor frequency over the momentum distribution).
Experimental estimates of τ s and τ p are generally obtained in a phenomenological way by fitting the experimental resistivity curves to the theoretical formulae obtained using semi-classical transport equations [12,17]. However, this phenomenological analysis is not well connected with the microscopic interpretation. First of all, the weak SOC in graphene would suggest τ s in the microsecond range [18][19][20], in clear disagreement with experimental data. In addition, the τ s estimated in high-mobility graphene with long mean free paths remains unsatisfactorily interpreted with a single relaxation mechanism, say EY or DP [13,21,22]. The suppression of τ s in clean graphene has been tentatively associated to an enhanced (intrinsic or extrinsic) spin-orbit coupling due to mechanical deformations such as ripples [23], or unavoidable adatoms incorporated during the device fabrication process [24,25], but the ultimate and microscopic nature of spin relaxation at play remains controversial and elusive [26].
Here, we explore this key fundamental issue by investigating the effect of weak perturbation induced by low densities of ad-atoms (down to 10 12 cm −2 ), which introduce a random Rashba field in real space but vanishingly small intervalley scattering, yielding long mean free paths. For typical electron densities within [10 10 , 10 12 ]cm −2 , the Fermi wavelength (λ F = 2 π/n, n the charge density) lies between 20 and 200 nm and thus exceeds the mean separation between adatoms (∼ 10nm) where spin-orbit scattering occurs, which questions the use of a standard semiclassical description. To study spin dynamics (and spin relaxation), we use a non-perturbative method by solving the full time-dependent evolution of initially spin polarized wavepackets, either through a direct diagonalization of a continuum model, or a real space algorithm and a tight-binding model for a microscopic disorder. We describe the system of a graphene monolayer functionalized with a random distribution of adatoms. The electronic structure of clean graphene is captured by the usual π-π* orthogonal tight-binding model (with a single p z -orbital per carbon site, zero onsite energies and nearest neighbors hopping γ 0 ). The presence of non-magnetic adatoms randomly adsorbed at the hollow positions on the graphene sheet introduces additional local spin-orbit coupling terms ( Fig. 1a,b), defined as [27].
The first term is the nearest neighbor hopping with γ 0 = 2.7 eV. The second term is a complex next nearest neighbor hopping term which represents the intrinsic SOC induced by adatoms, with d kj and d ik the unit vectors along the two bonds connecting second neighbors, s is a vector defined by the Pauli matrices (s x , s y , s z ), and V I the intrinsic SOC strength. The third term describes the Rashba SOC (V R ) which explicitly violates z → − z symmetry, with z being a unit vector normal to the graphene plane.
The last term denotes a potential shift µ associated with the carbon atoms in the random plaquettes R adjacent to adatoms (Fig. 1b). Such shift is due to weak electrostatic effects that arise from charge redistribution induced very locally around the adatom [27].
A Rashba splitting has been observed experimentally at the graphene/nickel and graphene/gold (Au) interfaces with spin splitting of up to 100 meV [28,29]. Gold and nickel as well as other materials like titanium, cobalt or chromium, are usually present during the fabrication of the nonlocal spin valves that are used to determine τ s and likely leave residues on the exposed graphene surface. Hereafter, we consider the case of Au adatoms whose influence on the transport properties of graphene has been studied experimentally [24]. The tight-binding parameters to describe both intrinsic and Rashba SOC induced by such adatoms are extracted from ab-initio calculations [29]. We then explore how the spin relaxation times scale as a function of the adatom density and adatom-induced local potential shift.
The spin dynamics is investigated by computing the time-dependence of the spin polarization defined by assuming that spins are initially injected out-of-plane (z direction), i.e.|Ψ(t = 0) =|ψ ↑ . The time evolution of the wavepackets |Ψ(t) is obtained by solving the time-dependent Schrödinger equation and are evaluated from the spreading of wavepackets by using real space propagation methods [30]. We focus on the expectation value of the spin z-component Figure 1 shows the typical behavior of S z (E, t) for two selected energies (at the Dirac point and at E = 150 meV) and two adatom densities ρ = 0.05% (about 10 12 adatoms per cm 2 ) (c) and ρ = 8%  gives τ s and T Ω extracted from the fits of S z (E, t) for varying adatom density. One first observes that the spin precession period is energy independent and is precisely equal to T Ω = πh/λ R (withλ R = 3ρV R an average SOC strength) even for the lowest coverage, which agrees with the estimate based on the continuum model [19]. In contrast, the spin relaxation time displays a significant energy dependence.
A V-shape is obtained for low energy, with τ s being minimal at the Dirac point with values ranging from 0.1 ps to 200 ps when tuning the adatom density from 8% to 0.05% (as given in Fig.3a, main frame).
Based on the observed scaling τ s ∼ 1/ρ (see Fig. 3b), one can further extrapolate the spin relaxation times for even smaller defect density, obtaining τ s ∼ 1 − 10ns for adsorbate densities decreasing from 10 11 cm −2 down to 10 10 cm −2 . The obtained V-shaped energy dependence and the absolute values of τ s are remarkably similar to those reported experimentally [11,12,17,24].
The faster relaxation at the Dirac point is actually evident in Figs. 1c and 1 d. The reason for this behaviour is the decay of the coupling between the pseudospin and momentum and the enhanced contribution of the SOC interaction, which leads to spin-pseudospin entanglement. The details of the entanglement are further described in Eq. (3) below and the Supplementary Information.
As discussed above, the usual approach to discriminate between conventional Elliot-Yafet and Dyakonov-Perel relaxation mechanisms in metals and semiconductors is to scrutinize the scaling of τ s versus τ p . Such procedure does not necessarily apply if the dominant processes that lead to momentum and scattering relaxation are not the same. For instance, in 2D membranes that respect mirror inversion symmetry, it was demonstrated that the carrier scattering by flexural phonons leads to fast spin flips but not to momentum scattering and, therefore, the spin transport is decoupled from the carrier mobility [23]. In the following discussion, we show that simple EY or DP scaling is also not suitable to describe our findings.
Within our microscopic calculations, we analyze the time-dependence of the diffusion coefficient for varying energies and ad-atom densities (Fig. 2c,d). For the lowest impurity density (0.05%, Fig. 2c), regardless of the considered energy, D(E, t) is seen to increase in time with no sign of saturation within our computational capability, indicating a ballistic-like regime for the considered timescales. Only for the largest ad-atom density (8%) does D(t) eventually saturate at high enough energies (above 100 meV, D(t) → D max ), allowing for the evaluation of the transport time using τ p (E) = D max (E)/2v 2 (E) (see dashed lines in Fig. 2b). A sharp increase of τ p is seen when approaching the Dirac point, where τ s reaches its minimum value, with τ s ≪ τ p . This energy dependence in τ p is not unique to gold adatoms but has also been observed for other types of disorder with a weak intervalley scattering contribution, such as epoxide defects or long range scatterers [30]. As seen in Fig. 3b, τ s ∼ 1/ρ, which does not allow us to discriminate between EY and DP processes. However, the absolute values of τ s and τ p (with τ s ≪ τ p ) are a clear manifestation of the breakdown of the typical scaling associated to both mechanisms. Even the unconventional DP regime described in Ref. [8] for the case of τ p /T Ω ≥ 1 where 1/τ s ∼ ∆Ω (with ∆Ω an effective width of the distribution of precession frequencies) cannot account for the observation that a weak variation in the local disorder affects the absolute values of τ s (while ρ is unchanged) as observed in Fig. 2. Here local disorder is monitored by the µ parameter. (Although µ belongs to the TB parameterization of the adatom, we use it temporarily to increase local disorder.) In fact, its value could slightly change when modifying the substrate screening or in presence of a more strongly bonded adsorbant than Au. As a consequence of the above findings, the spin relaxation mechanism at play is incompatible with both the EY and the DP mechanisms, a fact which could shed new light on the current debate on the microscopic nature of spin relaxation in clean graphene [13,21,22].
We now further study the origin of the τ s minimum at low energy and its unconventional scaling with τ p . Given that our simulations with the microscopic model give τ s ≪ τ p , we further explore the lowenergy spin dynamics with an effective continuum model, in which the spin-orbit scattering is treated as a homogeneous potential [19]. We solve the Dirac equation in the continuum model by using a 4× 4 effective Hamiltonian, taking into account the pseudospin degree of freedom While the hopping from three nearest neighbors h 0 ( k) =hv F (ζσ x k x + σ y k y ) ⊗ 1 s dominates at high energy and vanishes at the Dirac point (ζ = ±1 for K and K ′ valleys, σ are pseudospin Pauli matrices and 1 s is a 2 × 2 identity matrix), the intrinsic SOC h I ( k) =λ I ζ [σ z ⊗ s z ] and the Rashba interaction Within the continuum model spin relaxation is achieved by introducing an ad-hoc energy broadening.
We use an initially z-polarized state for injection and consider only the K valley. A certain density of Au impurities (inducing local spin-orbit coupling) is described by the effective spin-orbit couplinḡ λ R = 3ρV R andλ I = 3 √ 3ρV I . Note that no additional local (static) scattering potential is introduced here (µ = 0). By computing the spin dynamics of initially spin-polarized wavepackets, one also obtains a spin relaxation effect defined by the two timescales T Ω and τ s (see Supplementary Material).
It is instructive to contrast the results of the continuum model (Fig. 3a, inset) with those from the microscopic model (Fig. 3a, main frame). Although the spin precession period T Ω obtained by both models is identical (Fig. 3b) and the energy dependence of τ s is similar, the absolute values of τ s differ substantially, especially in the high energy regime, where τ s is clearly overestimated using the continuum model. Such difference also becomes increasingly large upon decreasing the adatom density because τ s presents a different scaling with defect coverage (Fig.3b). This clearly evidences the importance of disorder and illustrates the limits of a phenomenological approach using the continuum model for a quantitative comparison with experimental data. Notwithstanding, the qualitative agreement between both models (particularly for high coverage) and the weak momentum relaxation effects observed in the microscopic model (as seen in the long τ p ) suggest some generality in the unconventional spin relaxation observed near the Dirac point.
To further substantiate the origin of the spin relaxation, we scrutinize the spin and pseudospin dynamics of wavepackets using the continuum model. Pseudospin is intrinsically related to the graphene sublattice degeneracy and, as long as valley mixing is negligible, pseudospin is aligned in the direction of the momentum at high energy (h 0 ( k) dominates the Hamiltonian (3)). The Rashba spin-orbit term h R ( k) entangles spin s with the lattice pseudospin σ, overriding the locking rule between pseudospin and momentum since h 0 ( k) becomes vanishingly small in the vicinity of the Dirac point [7,15]. regularly showing an oscillatory pattern of S z (t) dominated by a single period T Ω = πh/λ R =0.19 ps (Fig. 4a). The spin precession occurs about an effective magnetic field B R dictated by the Rashba SOC and pointing tangentially to the Fermi circle (as seen from the precession from blue to pink in middle panels from t 1 to t 4 ). In contrast, the pseudospin σ(t) points approximately in the same direction of the momentum (evolving from green to orange). Its oscillatory pattern is driven by the Rashba period T Ω together with a superimposed and more rapid oscillation (see Supplementary Material).
The situation at low energy is markedly different (Fig. 4b,c). We observe a highly unconventional spin and pseudospin motion which is analyzed more closely for the spin and pseudospin z-components at the Dirac point and at E = −5 meV. Here, the amplitude of the pseudospin oscillation is strongly enhanced since pseudospin is no longer locked with momentum but starts to precess about an effective pseudo-magnetic field. The pseudo-magnetic field strongly depends on the spin orientation, thus yielding complex time-dependent dynamics of spin and pseudospin (see middle panels of Figure 4 corresponding to 4b,c). Such an effect derives from the increased pseudospin precession period T ps 0 = πh/E (about B ps 0 ), which decays significantly at low energy. Therefore σ i can no longer be replaced by its time average σ i , which in consequence also holds for the Rashba field B R . The time dependence of B R with variability on a timescale similar to the Rashba period leads then to strong non-linear dynamics of spin and pseudospin motion. As a result of such coupled dynamics, the spin precession cannot be described by a single period T Ω as becomes evident from the complex Fourier spectra of S z (t) in Fig.   4d. The time dependence of B R includes also changes of its direction, thus impacting the pseudospin and lifting the pseudospin-momentum locking. Both of these effects finally produce a joint spin/pseudospin motion prohibiting the de-coupling of driving forces (B ps 0 , B R ) that was possible at higher energies. In conclusion, our spin transport study in chemically modified graphene has revealed a hitherto unknown phenomenon related to the entangled dynamics of spin and pseudospin, induced by spin-orbit coupling and leading to fast spin relaxation in a quasi-ballistic transport regime. Entanglement between spin and orbital degrees of freedom has been discussed for ballistic semiconducting nanowires [31]. Here, faster spin relaxation develops when spin-pseudospin entanglement is maximized at the Dirac point, where the momentum scattering time becomes increasingly large because disorder preserves pseudospin symmetry. Such mechanism, occurring in clean graphene with long mean free paths, has no equivalent in condensed matter and cannot be described by EY or DP. It is here described for gold adsorbates, but should also be at play for other sources of local spin-orbit coupling (ripples, defects, etc.), thus contributing to the understanding of spin transport in graphene devices [11][12][13][14]17].
where the density matrix ρ(0) accounts for the initial spin polarization (out-of-plane) and s(t) = e iHt h se −iHt h is the spin operator in Heisenberg representation. As any trace, it can be replaced by the expectation value with random-phase states according to Tr [· · · ] → ϕ ′ RP |· · ·|ϕ ′ RP where the random-  |j is not an energy eigenstate at Fermi energy but samples the full spectrum. The trace in Eq. (4) (and equally the form with |ϕ ′ RP ) includes all states of the system at higher and lower energies, and not only those accessible in transport experiments (which are restricted to the Fermi energy). Accordingly, Eq. (4) is not appropriate when aiming at a comparison to experiment and another quantity needs to be computed as explained below. The quantum average of a given operator Q at a selected energy E can be generally written as average over all eigenstates at this energy through where |ψ i E are N degenerate eigenstates of H at energy E which are obtained by Hamiltonian diagonalization. It is next straightforward to show that where δ(E − H) is the continous projection operator, and the normalization with the number of states the trace with spin-polarized initial random phase states |ϕ RP = 1 where the time evolution of the wavepackets |Ψ(t) = e Real space implementation of the wavepacket quantum dynamics. The transport times are deduced from the numerical analysis of the spreading of wavepacket through [30,33]: A key quantity in the analysis of the transport properties is the diffusion coefficient: Assuming the system to be isotropic for the in-plane x and y directions, then The diffusion coefficient contains all information about the semiclassical effects of scattering leading to diffusive behavior, but also the quantum interference effects which lead to localization effects. D(t) increases ballistically at short times, then saturates due to elastic scattering events, and finally decays as a result of quantum interference effects (when significant). The elastic mean free path is derived from the maximum of the diffusion coefficient: with v(E) being the carrier velocity and D max the maximum value of D(t). The momentum relaxation can be extracted from elastic mean free path τ p (E) = ℓ e (E)/v(E).