# Microscopic Theory of Spin Relaxation Anisotropy in Graphene with Proximity-Induced Spin–Orbit Coupling

###### Abstract

Inducing sizable spin–orbit interactions in graphene by proximity effect is establishing as a successful route to harnessing two-dimensional Dirac fermions for spintronics. Semiconducting transition metal dichalcogenides (TMDs) are an ideal complement to graphene because of their strong intrinsic spin–orbit coupling (SOC) and spin/valley-selective light absorption, which allows all-optical spin injection into graphene. In this study, we present a microscopic theory of spin dynamics in weakly disordered graphene samples subject to uniform proximity-induced SOC as realized in graphene/TMD bilayers. A time-dependent perturbative treatment is employed to derive spin Bloch equations governing the spin dynamics at high electronic density. Various scenarios are predicted, depending on a delicate competition between interface-induced Bychkov-Rashba and spin–valley (Zeeman-type) interactions and the ratio of intra- to inter-valley scattering rates. For weak SOC compared to the disorder-induced quasiparticle broadening, the anisotropy ratio of out-of-plane to in-plane spin lifetimes agrees qualitatively with a toy model of spins in a weak fluctuating SOC field recently proposed by Cummings and co-workers [PRL 119, 206601 (2017)]. In the opposite regime of well-resolved SOC, qualitatively different formulae are obtained, which can be tested in ultra-clean heterostructures characterized by uniform proximity-induced SOC in the graphene layer.

## I Introduction

The tailored control of electronic properties in van der Waals heterostructures built from the assembly of two-dimensional (2D) crystals has provided a unique route to explore interface-induced phenomena (Belleste_Nanoscale2011, ; Butler_ACSNano2013, ; Zhang_JPDAppl2017, ). Heterostructures combining graphene and semiconducting group-VI dichalcogenides [MX (e.g., MMo, W; XS, Se)] could enable low-power spin-logic devices harnessing the unique interplay between quantum (spin and valley) degrees of freedom in honeycomb layers (Zibouche14, ; Soumyanarayanan_16, ; Garcia18, ). This thrust has been fueled by the prospect of enhancing spin–orbital effects in graphene (Geim_Nature2013, ; Soumyanarayanan_Nature2016, ), while preserving the quintessential Dirac character of its 2D quasiparticles. The much sought after interface-induced SOC has been recently demonstrated in graphene/TMD bilayer heterostructures (Avsar_NatComm2014, ; Wang_NatComm2015, ; Wang_PRX2016, ; Yang_2DMat2016, ; Volk_PRB2017, ; Zihlmann_PRB2018, ), where sharp weak antilocalization features in the magnetoconductance data (Wang_PRX2016, ; Yang_2DMat2016, ; Volk_PRB2017, ; Zihlmann_PRB2018, ) and dramatic reduction of spin lifetimes (Raes_NatComm2016, ; Benitez_Nature2017, ; Ghiasi_NanoLett2017, ) hint at a massive enhancement of spin–orbit interactions in the 2D carbon layer (up to 10 meV), consistent with the predictions of model calculations and first-principles studies (Wang_NatComm2015, ; Gmitra_PRB2016, ; Alsharari_PRBR2016, ).

The modification of electronic states in graphene-based van der Waals heterostructures due to proximity-induced SOC can be understood within a weak interlayer coupling picture, where Dirac states located in the band gap of a 2D semiconductor are perturbed in two fundamental ways. Firstly, the interfacial breaking of mirror inversion symmetry leads to the familiar Bychkov-Rashba effect (BychkovRashba84, ). The spin rotational invariance is lifted (point group symmetry reduction ), which causes the spin splitting of the Dirac states. Secondly, the proximity to different atoms (metal or chalcogen elements) located beneath the graphene flake () effectively “transfers” the sublattice-resolved SOC of the TMD substrate onto graphene (and hence spin–valley interactions). The relative magnitude of the spin–orbit effects experienced by -electrons in graphene depend on type and number of TMD layers, degree of vertical strain, and possible presence of resonant spin–orbit scatterers (Gmitra16, ; Singh_GTMDs2018, ; Pachoud15, ; Huang16, ; Wakamura18, ). The proximity spin–orbital effects couple all internal degrees of freedom of graphene (i.e. spin, sublattice and valley), enabling interesting spin-dependent non-equilibrium phenomena, including highly anisotropic spin dynamics (Cummings_PRL2017, ), spin-Galvanic and spin-Hall effects (Offidani_PRL2017, ; Milletari_PRL2017, ; Garcia_NanoLett2017, ).

In this work, we investigate how spin relaxation times in weakly disordered monolayer graphene are affected by proximity-induced SOC. The spin–orbit (SO) interaction enters the long-wavelength continuum Hamiltonian as an additional uniform term , that is (we choose natural units with )

(1) |

where is the Fermi velocity of massless Dirac fermions and is a disorder potential describing scattering from nonmagnetic impurities. The Hamiltonian is expressed in the basis and we have introduced () with as Pauli matrices in the valley (sublattice) space, respectively (here, and denote identity matrices). While knowing exactly the SO interaction is generally not possible, first-principles calculations and transport data provide a mean to estimate the various SO terms allowed by symmetry (Bir_Book1974, ; Winkler_Book2003, ; Kochan_PRB2017, ; Cysne_PRB2018, ). It is straightforward to show that there are only three such terms compatible with symmetry, , respectively, intrinsic-like SOC (McClure_1962, ; Kane_2005, ), Bychkov-Rashba SOC (Bychov_1984, ) and spin–valley interaction (Gmitra_PRB2016, ; Kochan_PRB2017, ). We note in passing that, beyond SOC, charge carriers in graphene can also experience an orbital sublattice-staggered potential (Gmitra_PRB2016, ). This effect is believed to be very weak in graphene/TMD bilayers (in contrast to rotationally aligned graphene on h-BN (Jung_NanoLett2012, )) and will be neglected in the following discussion (Foot_1, ).

The intrinsic-type SOC reads

(2) |

where is the spin–orbit energy. This term is invariant under all symmetry operations of the group, and thus it is already present in pristine graphene. As shown in a seminal work by Kane and Mele (Kane_2005, ), a large would drive graphene into a nontrivial topological insulating phase. However this term is very weak in graphene on typical substrates (Min_PRB2006, ; Hernando_PRB2006, ; Gmitra_PRB2009, ). Furthermore, in 2D heterostructures, the interfacial breaking of mirror inversion symmetry favours the appearance of an in-plane pseudo-magnetic field, that is, the familiar Bychkov-Rashba effect. This term (invariant under the point group) directly couples to the velocity of electrons, thus acting as a Lorentz pseudomagnetic field (Milletari_PRL2017, ):

(3) |

Finally, in honeycomb layers with interpenetrating triangular lattices made up of chemically distinct species, another spin-conserving SOC is allowed (Zhu_PRB2011, ; Xu_NatPhys2014, ). The sublattice inversion asymmetry can be captured by introducing sublattice-resolved next-nearest neighbours hoppings reducing the point group symmetry to (Kochan_PRB2017, ). This leads to a Zeeman-type spin-valley coupling

(4) |

The scenario faithfully describes graphene on TMDs, where the small lattice mismatch produces different SO energy on carbon sublattices (Wang_NatComm2015, ; Kochan_PRB2017, ; Gmitra_PRB2016, ). Uniform proximity-induced SO terms are block diagonal in valley space due to absence of interlayer hoppings connecting inequivalent valleys in graphene (Phong_2Dmat2017, ; Cysne_PRB2018, ).

Interface-induced Bychkov-Rashba and spin-valley interactions in graphene/TMD bilayers can in principle be large as tens of meV. With such a sizable imprinted in-plane (Lorentz-type) and out-of-plane (Zeeman-type) SO fields, the spin relaxation times for in-plane () and out-of-the-plane () polarization channels can be dramatically different. A recently-introduced figure of merit for the competition of the SOC along orthogonal spatial directions is the spin relaxation time anisotropy (SRTA): , which in graphene on TMDs has been estimated to be of order (Raes_NatComm2016, ; Benitez_Nature2017, ; Ghiasi_NanoLett2017, ). A simple treatment to obtain SRTA ratios has been put forward in Ref. (Fabian_Acta2007, ), which assumes that the electronic motion of bare quasiparticles (without SOC) is affected by a perturbing spin–orbit field with its precession axis randomly changing due to impurity scattering. The model applied to graphene/TMDs systems yields analytic formulas relating to the ratio and , where and are, respectively, the intra- and inter-valley momentum lifetimes (Cummings_PRL2017, ). However, the formalism presented there is limited to weak SOC, that is, with . This can be a strong constraint when trying to model ultra-clean samples with high charge carrier mobility, in which can be as large as unity (Wang_PRX2016, ). Also, a microscopic approach able to provide more physical insight to spin relaxation would be desirable, calling for a detailed study of how the spin dynamics is affected by the interplay of uniform proximity-induced SOC and impurity scattering. Here, we address theoretically this problem by means of the single-particle density matrix formalism. We obtain a set of coupled spin Bloch equations governing the spin dynamics for high electronic density — being the Fermi energy—assuming Gaussian-type (white-noise) disorder leading to intra- and inter-valley scattering processes. A variety of scenarios is shown to emerge, from simple- or multi-exponentially decaying spin dynamics to purely damped oscillating modes, depending on the relative magnitude of the three main energy scales: , and . We provide analytic expressions for the SRTA in the asymptotic limits of weak SOC (compatible with the findings of Ref. (Cummings_PRL2017, )) and strong SOC, which should be used to fit experimental data when .

The paper is organized as follows. Section II derives the general spin Bloch equations starting from the quantum Liouville equation. In Sec. III, we provide analytic solutions in the presence or absence of intervalley scattering and in the limiting cases of weak and strong SOC. Section IV discusses the obtained SRTA, putting it in relation with recent theoretical and experimental results and Sec. IV presents our conclusions.

## Ii Formalism: Spin bloch equations

The starting point of our approach is the quantum Liouville equation for the single-particle density matrix operator (Tarasenko06, ; Culcer_PRB2007_2008, ; Dugaev09, ; Wu_2010, )

(5) |

We consider a scattering potential generated by dilute short-range impurities at random locations ,

(6) |

where () are reals parameterizing the amplitude of intravalley (intervalley) scattering processes and characterize the spatial profile of the scattering potential. In order to derive the spin Bloch equations for high electronic density, we follow closely the treatment by Culcer and Winkler (Culcer_PRB2007_2008, ). The first step is to project Eq. (5) onto plane-wave eigenstates of the unperturbed graphene Hamiltonian, namely

(7) |

where is the wavevector around a Dirac point ( is the wavevector angle) and are quantum indices for sublattice, valley and spin, respectively. The free eigenvalues read as , where . is then a matrix of dimension , whose matrix elements are written as and is short-hand for the set of quantum indices (we use a similar notation for , and ). The proximity-induced SOC term has non-zero matrix elements between conduction and valence states leading to interband transitions. However, we focus here on the large Fermi energy regime , where interband coherence effects are strongly suppressed. Hence, we take . For simplicity of notation, we consider positive energies , henceforth considering electrons in the conduction band and dropping the sublattice index from all expressions. To simplify the treatment we also neglect valley coherence (Foot_2, ). The two inequivalent Dirac points can only be connected then by scattering events, according to Eq. (6).

Following Ref. (Culcer_PRB2007_2008, ), we split the density matrix into diagonal and off-diagonal elements: where for it is assumed . We have

(8) | ||||

(9) |

To simplify the analytical treatment, we neglect the term in the commutator on the left-hand side of Eq. (9). The approximation is valid in the limit of high Fermi energy, that is, . Also, only contains off-diagonal elements in , such that the commutator on the right-hand side of Eq. (8) only contains . We are ultimately interested in the diagonal part , as the spin observables are defined as

(10) |

We hence solve Eq. (9) and substitute the solution into the right-hand side of Eq. (8), which gives the collision integral. As customary, we treat Eq. (9) perturbatively for weak disorder with Gaussian (white-noise) statistics

(11) | ||||

(12) |

where is the impurity areal density. After a somewhat lengthy but straightforward calculation, where Eqs. (8)-(9) are expressed in the interaction picture and the evolution operator is expanded in powers of , one arrives at the following equation for the spin components

(13) |

with a Larmor precession term

(14) |

A few comments are in order. Central to the derivation of the quantum kinetic equation for the reduced spin density matrix [Eq. (13)] is the assumption of Gaussian disorder. The latter is equivalent to the first Born approximation (Offidani_MDPI2018, ) and thus it neglects any effects from skew scattering (allowed in the model (Milletari_PRL2017, )) and modifications to the energy dependence of the collision integral due to scattering resonances. Nevertheless, the relation between spin lifetime and momentum scattering time is expected to be preserved at all orders in perturbation theory, as shown explicitly in the minimal Dirac–Rashba model () with (Offidani_MDPI2018, ). This means that inclusion of higher-order scattering processes beyond the first Born approximation should not affect the SRTA ratios in the regime of validity of the quantum kinetic treatment (), consistently with the findings from exact numerical simulations (Cummings_PRL2017, ).

Next, we use the quantum kinetic equation Eq. (13) to obtain the spin Bloch equations governing the spin dynamics. Firstly, we separate the collision integral into intra and inter-valley parts, , with the corresponding matrix elements of the scattering potential

(15) | ||||

(16) |

where we have assumed that the impurity potential has a common matrix structure i.e., and (the generalization of our results to an arbitrary number of uncorrelated disorders can be easily accomplished using the standard Mathiessen’s rule). We can then write

(17) | ||||

(18) |

where . To solve the coupled system of 6 equations (3 spin 2 valley) Eq. (13), we expand in cylindric harmonics

(19) |

We note that the Dirac-delta function in Eq. (13) imposes energy conservation i.e., , such that the components of also depend on . Substituting Eq. (19) into Eq. (13), and retaining only the lowest-order harmonics we finally obtain (see Appendix for details)

(20) | ||||

(21) | ||||

(22) |

and

(23) | ||||

(24) | ||||

(25) |

where

(26) |

with . We have introduced the ratio of inter- to intra-valley energy scales defined as , as well as the intravalley momentum scattering time

(27) |

The spin Bloch equations [Eqs. (20)-(25)] together with the corresponding expressions for the barred component at —obtained by the formal replacement and —are the central result of this section.

## Iii Results

We are mostly interested in the zeroth harmonics of the various spin components, which according to Eq. (10) completely determine the spin density observables (foot_3, ). In most cases it is not possible to derive a simple closed expressions for arbitrary . Therefore in the following we solve the equations in the two limiting cases and , which is also helpful to get physical insight.

### iii.1 Intravalley scattering only:

The calculations are carried out explicitly for the out-of-plane component . The spin Bloch equations are recast in the following form

(28) |

where we introduced the following admixtures of in-plane spin harmonics

(29) | ||||

(30) |

The eigenfunctions can be written as

(31) |

where

are the the solution of the algebraic equation

(32) |

and are the corresponding eigenvectors. The coefficients are determined by imposing the Cauchy boundary conditions . The analytical solution to Eq. (32) is rather cumbersome. It is more transparent instead to find a solution perturbatively by expanding

(33) |

where and () representing the case of dominant Bychkov-Rashba (spin-valley) spin–orbit interaction. We find for

(34) |

where and

(35) |

For the minimal Dirac–Rashba model with , we recover the familiar Dyakonov-Perel relation (Zhang_NJP2012, ), resulting in an exponentially decaying solution with spin relaxation time

(36) |

In the latter regime, the spin polarization is lost due to motional narrowing, yielding its characteristic dependence on the momentum scattering time (see e.g., Refs. (Fabian_Acta2007, ; Wu_2010, )). In the opposite limit of resolved spin-splitting , electrons complete full Larmor coherent precession cycles between scattering events, which induce spin-memory loss (see Fig. 1 and discussion below). In this limit, the spin lifetime is of the order of the momentum scattering time, similarly to two-dimensional electron gases with large spin splitting (Gridnev01, ; Schwab_PRB2006, ; Liu_PRB2011, ). Combining the two limiting cases, we have

(37) |

For dominant spin–valley SOC (), we find instead

(38) |

which provides the asymptotic behaviour

(39) |

For the in-plane component a similar procedure leads to

(40) |

and

(41) |

Interestingly, the two weak SOC limits and display the same spin dynamics. The spin–valley term only provides a small correction to the Dyakonov-Perel spin-relaxation time. From these results, the SRTA ratio for pure intravalley disorder is readily obtained

(42) |

that is

the conventional SRTA of the minimal () model with weak Rashba SOC, i.e. (Fabian_Acta2007, ) is observed irrespective of the spin-valley coupling. On the contrary, in the case of strong Bychkov-Rashba SOC, the quantum kinetic treatment predicts . This result is related to the role of the Bychkov-Rashba field in the two opposite limits and , cf. Fig. 1. Note that because of the totally-in plane Bychkov-Rashba SOC, simple commutator algebra for the precession term gives that the period of out-of-plane spins is half of that of in-plane ones: . In the Dyakonov-Perel limit , where electrons’ spin only precesses a small angle before being scattered, the spin dynamics can be understood as the result of a random walk, with unit step . Spin relaxation is achieved after collisions when the accumulated phase is of the orde of unit, that is, . The faster precession of reflects in a different unit step , which immediately implies , i.e. spins along reach the critical value in half of the time compared to initially in-plane spins. On the contrary, when , spin relaxation is achieved on the time scale of a single impurity scattering event, as spins can coherently complete many precession cycles on a time scale . The anisotropic spin precession reflects instead in this case in the oscillating term rather than the spin decay, as found in Eqs. (37) and (40).

### iii.2 Intervalley scattering case:

Short-range scatterers and atomically-sharp defects responsible for a finite intervalley scattering time are invariably present in realistic conditions (Wang_PRX2016, ). Thus, the inclusion of intervalley processes in the collision integral is crucial to understand the spin dynamics in graphene-based heterostructures. Let us start again from the out-of-plane component. We find in this case

(43) |

with

(44) |

where we have identified the intervalley momentum lifetime

(45) |

Proceeding as shown above, we obtain after standard algebraic manipulations

(46) |

(47) |

Considering the in-plane components, we were able to reduce the initial set of 8 coupled equations to two equations coupling and (see Appendix for details), reading as

(48) |

where we have set

(49) | ||||

(50) |

Solving Eq. (48) with the same boundary conditions as above, i.e. and all the other functions being zero at the initial time, we find

(51) |

and

(52) |

with . In the large limit, second line of the latter equation, the solution includes a second term with , giving overall a multi-exponential solution. This term is subleading in the cases we are interested in, hence we neglected it in Eq. (52). In Figs. 2 and (3), we show representative examples of the spin polarization dynamics in the strong and weak SOC limits, respectively, according to our results.

## Iv Spin relaxation Time Anisotropy

We discuss now in more detail how the spin dynamics evolves from weak proximity-induced SOC () to well-resolved SOC (). The explicit form of the SRTA ratio is

(53) |

Together with the microscopic derivation of the spin Bloch equation for this model, Eqs. (20)-(25) and their solution—showing a a crossover between a purely damped to oscillating damped spin dynamics—these are the most important results of this paper. The first observation concerns the strong Bychkov-Rashba case with , which can in principle be achieved in clean graphene-based heterostructures, where also the lattice mismatch is sizable enough to produce . Contrary to the other two presented cases (first and third lines of Eq. (53)), in this limit a direct estimation of or is not possible. Hence, whenever is measured, the extraction of other parameters from spin precession measurements alone should be considered unfeasible.

We focus in the following on the two more interesting cases and For the weak SOC case—first line of Eq. (53)—we report a visualization of the obtained result in terms of contour lines for fixed , see Fig. (4)(a). Our results agrees very well with the toy model supporting the numerical findings in Ref. (Cummings_PRL2017, ), i.e.

(54) | ||||

(55) |

Note the the different pre-factors in front of the second term with respect to the results obtained from the microscopic Hamiltonian, Eq. (53)), first line. The inset of Fig. 4(a) shows a detailed comparison for the case . Following the analysis performed in Ref. (Ghiasi_NanoLett2017, ) (), assuming for graphene/MoSe (Gmitra_PRB2016, ), a is obtained, which taking ps gives ps (against ps following Ref. 54). These estimtes (obtained from modeling of spin precession data for ) agree qualitatively well with typical relaxation times obtained from weak localization data (Wang_PRX2016, ; Tikhonenko_PRL2008, ).

However different scenarios are possible. For instance, in Ref. (Yang_2DMat2016, ), the authors estimate , with for graphene/WS heterostructures. In this case the weak SOC approximation might fail. In fact, assuming as above, using Eq. (54) from Ref. (Cummings_PRL2017, ) one would get an unphysical result , where the intervalley scattering time is shorter than the (intravalley) momentum scattering time. The usage of Eq. (53) in the limit of strong spin-valley (third line) then is needed. Using this relation, we estimate , pointing to dominant intravalley processes.

## V Conclusions

In this work, we investigated theoretically the spin dynamics in graphene with proximity-induced SOC. Starting from the quantum Liouville equation, we derived the effective spin Bloch equations governing the spin dynamics of 2D Dirac fermions subject to in-plane (Bychkov-Rashba) and out-of-plane (spin-valley) interactions. We discussed in detail the irreversible loss of spin information with origin in intra- and inter-valley scattering processes within the standard Gaussian approximation for the disorder potential, obtaining the time dependence of the spin polarization vector and associated spin-relaxation times. We finally discussed the interesting results for the spin relaxation-time anisotropy . The result reported Ref. (Cummings_PRL2017, ) for weak SOC is qualitatively reproduced by our microscopic theory. Crucially, we have shown that the weak SOC approximation to the spin relaxation anisotropy ratios fails in when the proximity-induced SOC is of the same order or larger than the disorder-induced quasiparticle broadening. Our results for well-resolved SOC then should be used to fit spin precession measurements.

We remark that the adopted formalism is only valid in the highly-doped regime of large Fermi energy, where it is assumed that SOC only induces Larmor precession. This is a strong assumption that might break down at low electronic density in samples with large interface-induced SOC of order meV. In that case the spin texture of the bands is well established; momentum is then strongly correlated with the direction of the psuedomagnetic field, which can favour or inhibit certain matrix elements of the scattering potential . For instance, intervalley scattering has been suggested detrimental for the out-of-plane spin component, as producing transitions between states with opposite Zeeman pseudomagnetic field (Cummings_PRL2017, ). However our treatment, where the kinetic equations are projected onto eigenstates of bare graphene, has not the capability of capturing such an effect (see Figs. 2(b), 3(b), where is virtually unaffected by the value of ). A possibility to incorporate the SOC self-consistently is by adopting the quantum diagrammatic formalism for Dirac fermions, according to the procedure outlined in Ref. (Offidani_MDPI2018, ). We will address this problem in a future publication.

## Appendix A: Details on the Derivation of the Spin Bloch Equations

In this appendix we report more details about the derivation of the spin Bloch equations [Eqs. (20)-(25)] starting from the collision integral in Eq. (13)

(56) | ||||

(57) |

Note the collision integral in Eq. (56) is diagonal in valley space, i.e. which was justified in the main text. Intervalley processes are still taken into account internally to the collision integral, i.e. by considering transitions of the type where electrons initially at are scattered at and then scattered back at .

For point-like impurities, the different matrix elements of the scattering potential are written as

(58) | ||||

(59) |

which plugged into Eq. (56) and after having taken after disorder average as prescribed in Eqs. (11),(12) gives Eqs. (17)-(18) of the main text. Using the notation in the main text and the relation

(60) |

we have explicitly

(61) | ||||

(62) | ||||

(63) |

with

(64) |

and for the intervalley part

(65) | ||||

(66) | ||||

(67) |

and

(68) |

The SOC couples different harmonics .. Let us neglect that for a while. We find the following system of equations

(69) | ||||