Spin Hall and inverse spin galvanic effects in graphene with strong interfacial spin-orbit coupling: a quasi-classical Green's function approach

van der Waals heterostructures assembled from atomically thin crystals are ideal model systems to study spin-orbital coupled transport because they exhibit a strong interplay between spin, lattice and valley degrees of freedom that can be manipulated by strain, electric bias and proximity effects. The recently predicted spin-helical regime in graphene on transition metal dichalcogenides, in which spin and pseudospin degrees of freedom are locked together [M. Offidani et al. Phys. Rev. Lett. 119, 196801 (2017)], suggests their potential application in spintronics. Here, by deriving an Eilenberger equation for the quasiclassical Green's function of two-dimensional Dirac fermions in the presence of} spin-orbit coupling\textcolor{black}{{} (SOC) and scalar disorder, we obtain analytical expressions for the dc spin galvanic susceptibility and spin Hall conductivity in the spin-helical regime. Our results disclose a sign change in the spin Hall angle (SHA) when the Fermi energy relative to the Dirac point matches the Bychkov-Rashba energy scale, irrespective of the magnitude of the spin-valley interaction imprinted on the graphene layer. The behavior of the SHA is connected to a reversal of the total internal angular momentum of Bloch electrons that reflects the spin-pseudospin entanglement induced by SOC. We also show that the charge-spin conversion reaches a maximum when the Fermi level lies at the edge of the spin-minority band in agreement with previous findings. Both features are fingerprints of spin-helical Dirac fermions and suggest a direct way to estimate the strength of proximity-induced SOC from transport data. The relevance of these findings for interpreting recent spin-charge conversion measurements in nonlocal spin-valve geometry is also discussed.


I. INTRODUCTION
van der Waals heterostructures [1] have become one of the most promising spintronic platforms, where both fundamental and applied aspects of spin transport can be addressed with exquisite electrical control in the atomically thin limit [2][3][4].
Soon after graphene became well established as a highperformance spin channel supporting spin transport over long distances at room temperature [5][6][7][8][9], the primary focus has shifted towards the study of emergent spincharge coupling effects in van der Waals heterostructures. An intriguing possibility consists of exploiting relativistic SOC phenomena to generate and manipulate spin-polarization flow in atomically thin planes. One of the first proposed schemes made use of proximityinduced SOC in graphene flakes with a dilute coverage of heavy adatoms [10][11][12], which are believed to be efficient extrinsic sources of spin Hall currents and nonequilibrium spin polarization [13,14]. An alternative approach, which recently meets a lot of attention, consists of enhancing the SOC by placing the graphene sheet on top of a layered semiconductor [15][16][17][18][19][20][21]. It is now understood that the breaking of inversion symmetry in van der Waals heterostructures results in dramatically enhanced intrinsic-and Bychkov-Rashba (BR)-type SOC [22][23][24], endowing spin-split Dirac cones with a robust skyrmionlike spin texture in k space [25]. The interface-induced SOC, whose precise spatial profile reflects the interlayer atomic registry and disorder landscape, can be further manipulated by applying strain and electric fields [26][27][28][29][30][31]. Low-temperature magnetotransport data for graphene on transition metal dichalcogenides (TMDs) is consistent with interface-induced SOC in the range 1-10 meV [32][33][34][35][36], up to two orders of magnitude higher than graphene's weak intrinsic SOC (λ I 42.2 µeV [37]), in good agreement with theoretical predictions based on density functional theory and semi-empirical methods [16,27,29]. Concurrently, room-temperature Hanle-type spin precession measurements revealed another fingerprint of proximitized graphene, that is, a giant ratio of out-of-plane to in-plane spin lifetimes (τ ⊥ /τ ≈ 10 [38][39][40]) driven by the competition of symmetry distinct spin-orbit interactions and intervalley scattering [41,42]. Meanwhile, quantum Hall effect measurements performed on ultraclean bilayer graphene/few-layer WSe 2 have shown that interfacial SOC can be made as large as 15 meV by removing contaminants from the device areas [20].
A high interfacial SOC with magnitude comparable to the quasiparticle lifetime broadening is a desirable feature because it allows efficient spin-charge conversion via spin-helical Dirac fermions [18] as demonstrated in a series of elegant spin-valve experiments on graphene/TMD bilayers [43][44][45][46][47]. The delicate interplay between intrinsic and extrinsic (impurity-driven) spin-orbit-coupled transport mechanisms in graphene-based heterostructures has been recently studied by means of a linear response formalism, supported by conservation laws [48][49][50][51]. Unlike BR-coupled two-dimensional (2D) electron gases [52][53][54], the spin Hall conductivity in an infinite system with non-magnetic defects was found to be finite due to the generally noncoplanar nature of the equilibrium spin texture at the Fermi energy [50]. The critical role played by impurity scattering in the context of SOC-driven spincharge conversion has also been investigated by means of the non-equilibrium Green's function technique, which is particularly suited to derive coupled spin-charge driftdiffusion equations [55][56][57]. In particular, the Keldysh technique in the so-called quasiclassical approximation, pioneered by Eilenberger [58,59] for dirty type II superconductors, has been applied to describe the locked spin-charge dynamics of topological insulators (TI) [60].
The aim of the present paper is to extend the quasiclassical approach developed in Refs. [55][56][57] to a system of dirty 2D Dirac fermions subject to strong proximityinduced SOC. Our focus will be on the low-density regime highlighted in Ref. [18], in which the Fermi energy crosses a single spin-split band, and thus the 2D Dirac fermions acquire a well-defined spin helicity akin to surface states of a three-dimensional topological insulator. Our work is organized as follows: In Section II we introduce the effective Hamiltonian of graphene/TMD heterostructures and discuss how proximity-induced SOC modifies the electronic structure at low energies. In Section III, we present the Keldysh technique and in Section IV we derive the Eilenberger equation for the quasiclassical Green's function in the spin-helical regime. In particular we discuss the T−matrix expansion for the disorder potential and derive the general expression for the collision integral. In Section V, by confining to the Born approximation for the self-energy, i.e. the second order in the T−matrix expansion, we solve the Eilenberger equation and find the expressions for the longitudinal electrical conductivity and the spin galvanic susceptibility, while there is no spin Hall effect due to the absence, at this order, of skew scattering. The latter is explicitly considered in Section VI by carrying out the T-matrix expansion up to the third order in the scattering potential. (Some technical details are relegated to the Appendix for clarity of exposition.) The explicit solution of the Eilenberger equation provides the expressions of the electrical conductivity, spin galvanic susceptibility and spin Hall conductivity in terms of the dimensionless disorder coupling strength and of the energy eigenstate at the Fermi level. Recent spin precession measurements of inverse spin Hall and spin galvanic effects in a graphene/WS 2 heterostructure [45] are put into context. Finally Section VII presents our conclusions.

II. THE MODEL
The low-energy excitations in graphene/TMD bilayers are governed by the following generalized Dirac-Rashba model (we use natural units = 1 throughout) where ξ = ±1 signs refer to inequivalent valleys K(+) and K (−), (Ψ ξ , Ψ † ξ ) ≡ (Ψ ξ (x), Ψ † ξ (x)) are 4-component spinor fields defined on the internal spaces of pseudospin and spin and v is the velocity of massless Dirac fermions (here, σ is the vector of Pauli matrices acting on the pseudospin space). The term V ξ SO describes the spatially uniform proximity-induced SOC on the graphene layer and comprises several contributions that reflect the C 3v point group symmetry of the heterostructure [16,29,61]. For instance, the breaking of mirror reflection symmetry about the graphene plane allows z → −z asymmetric SOC. This so-called BR effect [22] is generically present in graphene on a substrate and is responsible for the entanglement between pseudospin and spin degrees of freedom [23], with clear signatures in the spin dynamics [62] and current-induced spin-orbit torque [63]. Moreover, the spin-orbit interactions imprinted on the graphene layer acquire sublattice-resolved terms inherited from the noncentrosymmetric TMD layer, namely the spin-valley coupling, as mentioned in the Introduction. All together, in the longwavelenght limit, V ξ SO generally contains three SOC terms compatible with the C 3v point group i.e., where V KM ∝ σ z s z is the intrinsic-type Kane-Mele SOC [64,65], V BR ∝ (σ x s y − σ y s x ) describes the BR effect, [23] V ξ sv ∝ ξs z captures the spin-valley effect which acts on charge carriers as a valley-Zeeman coupling [66], and s is the vector of Pauli matrices acting on the spin space. In practice, the weak z → −z symmetric Kane-Mele SOC can be neglected in comparison to the other effects (see, for example, Ref. [27]) and hence in this work, we approximate the spin-orbit interaction as follows Beyond SOC, the proximity coupling to a TMD also imprints a sublattice-staggered potential H ξ ∆ = ∆ξσ z on the graphene sheet. The staggered on-site energy is believed to be substantially smaller than the SOC energy scales in Eq. (3) (under 0.1 meV according to a recent study [67]) and because it plays a very limited role in both the spin dynamics [42] and coupled spin-charge transport effects [18], we neglect it in the following discussion. We note that more SOC terms are allowed if additional crystal symmetries are broken further reducing the point group of the heterointerface. These include unconventional in-plane-polarized spin-valley and Kane-Mele-type SOC that are symmetry-allowed in graphene coupled to low-symmetry TMDs [29], the implications of which are beyond the scope of the present work. Finally, the last term in Eq.(1) is a random potential produced by scalar impurities, which will be responsible for the extrinsic generation of nonequilibrium spin polarization and spin Hall currents to be discussed in Secs. IV-VI.
The energy-momentum dispersion relation of the lowenergy Dirac bands reads as where M n (k) = [2λ 2 + λ 2 sv + 2n λ 4 + (λ 2 + λ 2 sv )v 2 k 2 ] 1/2 is the spin gap induced by SOC, n, l = ±1 are the band indices and k = |k| (with k measured from a K point). Equation (4) makes manifest the underlying particle-hole symmetry of the Hamiltonian, which results in one or two Dirac bands at the Fermi energy ε depending on the gatetunable charge carrier density (see Fig. 1). To ease the notation, we shall assume , λ sv , λ > 0 in what follows. Inverting Eq.(4) and evaluating the energy-momentum dispersion relation at the Fermi energy ln = ε, one obtains the Fermi momenta, k ln = k( ln = ε), as follows The ± sign depends on the energy range within which lies the Fermi energy ε (Fig. 1). The eigenvectors, ψ ln (x) = e ik·x Φ ln (k), have the following spinorial structure with θ the momentum angle with respect to an arbitrary axis in the graphene plane and α, β and γ given by In this work we specialize to the regime I which should be experimentally accessible in ultra-clean devices with strong interfacial SOC. The Fermi energy therefore lies in the interval , which we call spin-helical regime, where the simplyconnected Fermi surface topology is akin to ideal topological surface states [60]. For that reason, we will drop the band indices and define the Fermi momentum in regime I as k F ≡ k 1,−1 with k 1,−1 = k 1,−1 (ε) as per Eq. (5).
For the transport calculations below, we will also need the density of states at the Fermi level, per valley, From this expression, one easily recovers pristine graphene's well-known expression N 0 = ε/(2πv 2 ) by letting λ sv → 0 first and then λ → 0. We will show below that in this regime, where the electronic states have a well-defined spin helicity, the pseudospin and spin degrees of freedom are constrained in such a way that it becomes possible a full description of the coupled spincharge dynamics in terms of a single transport equation.

III. THE KELDYSH TECHNIQUE
The Keldysh formalism, which goes back to the pioneering works by Schwinger [68] and Keldysh [69], is a powerful generalization of the standard perturbative approach of equilibrium quantum field theory to nonequilibrium problems. Within the Keldysh technique [70][71][72][73][74], the Green's function has the following matrix structurě with G R , G A and G K the respective retarded, advanced and Keldysh components. The Green's function acts on the spin, valley and pseudospin spaces, and thus each block in Eq. (11) can be represented as a 8 × 8 matrix.
G satisfies the left-right subtracted Dyson equation [59], which conveniently gets rid of the delta-function singularity at coinciding space-time points, i.e.
where the square brackets define the commutator. Here the products are to be understood with respect to both the Keldysh and internal spin spaces and to the spacetime coordinates (i.e. convolution products). The kinetic and proximity-induced SOC terms are contained in the bare Green's functionǦ 0 , whereas the quasiparticle self-energyΣ is due to the disorder potential. In the following the self-energy will be obtained by averaging, term by term, the perturbative expansion of the left-right subtracted Dyson equation in the potential U (x) with respect to the disorder configuration. The Green's functioň G entering Eq. (12) is then the disorder-averaged Green's function. In explicit terms, we have and h (x) the Hamiltonian density as derived from the disorder-free part in Eq.
(1). The externally applied (uniform) electric field E is incorporated by means of the standard minimal coupling within the velocity gauge, −i∇ → −i∇ − eE t [59,70], t being the center-of-mass time defined below and e > 0 the unit of electric charge. The aim of this work is to derive a coupled spin-charge transport equation in the spin-helical regime. Following the standard approach, we define the center-of-mass and relative space-time coordinates As customary, we introduce Wigner-mixed coordinates by taking the Fourier transform with respect to the r and τ variables. The key assumption in the derivation of the transport equation is that the center-of-mass space-time variable x, t is a slow variable compared to r, τ . As a result one can perform a gradient expansion (see Appendix A for details) of Eq.(12) obtaining forǦ (x, t, k, ω) ≡Ǧ the following equation of motion where h (k) is the Fourier-transformed bare Hamiltonian density and the curly brackets denote, as usual, the anticommutator. In a compact form, Eq. (15) provides the Figure 2. Disorder self energy in the non-crossing approximation. In the Gaussian case, the self energy consists of a single "rainbow" diagram with two potential insertions (ΣG).
In the T -matrix approach, one effectively resums the full series of single-impurity scattering events (ΣT ), which is then proportional to the impurity areal density ni.
equation of motion for all the components of the Green's function Eq. (11). The Keldysh component in the right hand side term of Eq. (15) is usually named collision integral and, in the spirit of the semiclassical Boltzmann transport theory, it can be divided in a in-and out-term according to The retarded and advanced clean Green's functions are given by where P ln (k) = |Φ(k ln , θ) Φ(k ln , θ)| is the projector onto the eigenstate with indices ln. In the spin-helical regime only the projector P 1−1 (k) ≡ P (k) is relevant and we omit the band indices hereafter for simplicity.
In the presence of random impurities with areal density n i , the retarded and advanced Green's functions are dressed with the corresponding disorder self-energies. The self energy is given by the average over disorder (≺ ... ) of the T -matrix expansion as shown in Fig.  2, that is, In the above is the local retarded (advanced) clean Green's function, i.e. evaluated at coinciding space arguments and the disorder average is defined by and so on and so forth.
In the equation of motion, one determines the Green function self-consistently by replacing G R(A) with the disorder-average one G R(A) . To quadratic order in the disorder potential, in the so-called first Born approximation, one keeps only the first term of the series expansion of Eq. (18). The corresponding retarded/advanced selfenergies read as where from now . . . indicates the integration, normalized to 2π, over the wavevector angle defining the direction of k F . As a result the out−contribution of the collision integral becomes In a self-consistent evaluation of the self-energy, one should replace G R,A 0k with the disorder-averaged Green's function G R,A k in Eq. (19). Provided one is not too close to the Dirac point, one can neglect the disorder correction to the electron spectrum and the self-consistent solution for the self-energy is reasonably approximated by the expression in Eq. (19).
Next, we consider the Keldysh self-energy, which determines the in-contribution of the collision integral. The perturbation expansion of the Keldysh Green's function reads quite different from that of the retarded and advanced Green's functions. After disorder averaging, at each order in u 0 , there are as many terms as positions in which the Keldysh component can be placed with the additional requirement that on the left (right) of the Keldysh component there can be only retarded (advanced) Green's functions. For instance, at first and second orders one has In the end, the Keldysh self-energy then reads where we have introduced the following notation with Interestingly, only the T R,A parts of the T -matrices appear in the Keldysh self-energy upon performing the disorder average as a result of the latter involving impurity insertions in both the retarded and advanced Green's functions at the left and right of the Keldysh Green's function, respectively. This corresponds to the so-called T -matrix insertion in the vertex correction of the Kubo linear response formalism, which effectively resums the infinite set of single-impurity scattering diagrams [48,49]. The "in" contribution to the collision integral Eq.(16) finally reads After some algebra, it is possible to recast the collision integral Eq. (16) in the form (for details see Appendix B) From the point of view of the kinetic equation, the variables k and k represent the momentum before and after a single-impurity scattering event, depending whether one considers the in-and out-contributions. One can easily check that the detailed balance is satisfied. If k ↔ k , the in− and out−terms are interchanged and then the microscopic reversibility of the scattering probability is preserved [75]. The general expression Eq. (25) is one of the main results of this work. The familiar case of a Fermi gas with Gaussian white-noise correlation is recovered by letting T R(A) → 1.

IV. THE EILENBERGER EQUATION
In this section, we solve the kinetic equation presented earlier (cf. Eq.(15) and the following discussion) within the quasiclassical approximation [59,70,71]. The latter has a resemblance to the standard Boltzmann transport theory that lends itself to a physically transparent interpretation of 2D spin-orbit coupled transport effects in the experimentally relevant diffusive regime [25,48,63].
The established framework in which the quasiclassical theory is expressed is that of the real-time formulation of the Keldysh technique, as discussed above. The quasiclassical Green's functionǧ is usually defined by [58] g (x, t, n, ω) .
where we have introduced the variable ξ = 1−1 (k) − ε in terms of the dispersion relation of the spin-helical band at the Fermi energy. Since we assumed ε > 0, all the projectors in what follows are constructed from eigenstates l = 1, n = −1. The integration contour C is taken in such a way to capture the contribution of the Green's function pole, and that the ξ-integration leaves unaffected the angular dependence described by the unit vector n = k/ |k|. The so-called Eilenberger equation is then obtained by applying the ξ-integration to the Eq.(15) by reasonably assuming that the self-energy does not have a further singular behavior, which would add to the Green's function pole. At the leading order in the weak-disorder expansion ετ 1, with τ = 1/(2Im Σ A ) the quasiparticle lifetime, the ξ-integration procedure is not affected by the disorder-induced dressing of the pole. The Eilenberger equation then reads The quasiclassical Green's function has the same triangular matrix structure of the original Green's function, i.e.ǧ and in the clean system, the retarded (advanced) quasiclassical Green's function are where P (k F ) ≡ P (θ) is the projector on the spin-helical band evaluated at the Fermi energy with only the angle dependence remaining. In principle, we have to solve Eq. (27) for all the components ofǧ, but one can considerably simplify the problem by noting that the elastic scattering from scalar impurities merely produces a shift in the pole of the retarded (advanced) component. Hence, as already mentioned, the expression (29) remains unchanged at leading order in powers of (ετ ) −1 . Ultimately g R(A) gives information about the density of states whereas g K ≡ g provides the distribution function. Finally, the Eilenberger equation for the Keldysh component Eq. (27) becomes where the ξ-integration must be performed also in the right hand side of Eq. (15). Explicitly, by applying the ξ-integration to the expression of I given by Eq. (25) and anticipating the definition of the scattering time given below in the Eq. (32), we obtain where, for brevity, we defined g R;A;K ≡ g R,A,K (θ, ω) and its angle averaged g R,A,K = g R,A;K (θ , ω) with θ and θ the angles of k and k , respectively, both momenta evaluated at the Fermi surface. For convenience we have also introduced the basic scattering rate which has the meaning of the inverse quasiparticle lifetime of graphene without the SOC in the Gaussian limit.
By looking at Eq. (30) we note two key differences with respect to the semiclassical Boltzmann transport equation (BTE): (1) there is a commutator term i [h (k F ) , g] missing in the BTEs because there g is a scalar, (2) a general scattering kernel is present which, in principle, holds for any type of impurity (magnetic, etc.), whereas semiclassical BTEs based on a scalar distribution function only applies to scalar impurities. Although not considered in the present paper, the method can be readily extended to impurity potential with a matrix structure in the internal degrees of freedom. In the remainder of this paper, we confine our analysis to the case of a constant and uniform electric field for which the space and time derivatives in the Eq. (30) drop out and the quasiclassical Green's function then depends only on n, ω or equivalently on θ, ω. To leading order in E, one can safely replace g in the force term (in the left hand side) by the equilibrium Keldysh Green's function. The Eilenberger equation (30) then reduces to where g eq (θ, ω) = f (ω)P (θ) is the equilibrium quasiclassical Keldysh Green's function with where T is the temperature (in our natutal units k B = 1).
To proceed further, following Ref. [60], we propose the following ansatz for the quasiclassical Keldysh Green's function g(θ, ω) = g 0 (θ, ω) P (θ) , where g 0 (θ, ω) is a scalar function. Although g is still a matrix, its structure is entirely constrained. The ansatz can be motivated by the following argument. Inspection of the Eq. (30) shows that, at leading order in the dilute regime (ετ 1), the solution must commute with the Hamiltonian and be of order g ∼ τ ∼ τ 0 . The commutator in the left hand side, although of order ετ 1, vanishes because the Hamiltonian density at the Fermi level can be written as h (k F ) = εP (θ) and the remaining terms in the equation are of order unity. We note that this ansatz may not be sufficient when one is dealing with quantum (weak-localization) corrections in the weak-disorder expansion. Using these ingredients and by taking the trace of the Eq.(31) one gets a scalar collision integral for g 0 where is a function of the angle difference that from now we call ϑ ≡ θ − θ .
In the end, once the solution for the quasiclassical Green's function g 0 (θ, ω) is known one can easily obtain the steady-state observables, such as the electric current, the spin polarization and the spin current densities. According with the general recipe in the Abelian case [74,77], using the notation [59] dk if O = σ i ⊗ s j indicates a generic observable, its quantum statistical non-equilibrium average O is given by where O k = Tr [O P (θ)] /2 indicates the equilibrium average and, as before, . . . is the integration over the directions of k = k F n. The trace symbol involves all internal degrees of freedom: sublattice, valley and spin. Because intervalley scattering is neglected for simplicity (the impurity potential U (x) is a scalar quantity), the trace over the valley degree of freedom yields a simple factor of two, which is compensated by the factor of two in the denominator in the definition of the equilibrium average of O k . The relevant observables are the charge current (J = −evσ ⊗ s 0 ), spin polarization (S = 2 σ 0 ⊗ s) and z-polarized spin Hall current (J z s = v 2 σ ⊗ s z ). Here we have reinstated the Planck constant for the sake of clarity. For the C 3v -invariant model with spin-valley coupling, the equilibrium charge current, spin polarization density and (persistent) spin Hall current can be easily evaluated as where α, β and γ depend on the parameters of the Hamiltonian and on the Fermi energy as defined in the Eq.(6). Similar to the well-studied high density regime with two occupied Dirac bands [18,50,51,63], the standard first Born approximation will allows us to obtain the leading semiclassical contribution to the nonequilibrium spin polarization density generated by the inverse spin galvanic effect, namely S i ∝ ij τ E j , with ij the Levi-Civita symbol and i, j = x, y. On the other hand, the calculation of the steady-stateẑ-polarized spin Hall current density, J z i ∝ ij τ E j , will require higher-order corrections to the self-energy due to the pivotal role played by skew scattering mechanism in the extrinsic spin Hall effect.

V. THE EILENBERGER EQUATION IN THE FIRST BORN APPROXIMATION
In this section we consider the self-consistent Born approximation, which amounts to confine to second order in the disorder potential expansion, implying T R(A) −→ 1.
As a result, the matrix transport equation (33) reduces to a simpler scalar effective transport equation for the "charge" component g 0 (θ, ω) of the quasiclassical Green's function where is the scattering kernel in the first Born approximation [c.f. Eq. (37)]. We remind the reader that, unless otherwise stated, the projectors are evaluated at the Fermi surface. Due to the periodicity of the projectors, the scattering kernel can only depend on the difference of the two angles. Furthermore due to the cyclic property of the trace, the expression is invariant under the interchange of the two angles, implying that the effective kernel is an even periodic function of the angle difference ϑ ≡ θ − θ . The transport equation Eq.(43) together with the expressions of the effective kernel Eq.(44) and the velocity vertex Eq. (40), defines the coupled spin-charge response of the projected band at the Fermi energy. This is one of the central results of this paper and can be used to obtain the charge conductivity and spin-galvanic conductivity of the spin-helical band as discussed below. We first illustrate the formalism with the C 6v -invariant Dirac-Rashba model obtained by setting λ sv = 0 in Eq.
(3); the corresponding band structure and spin texture is shown in Fig. 1(a). In this case, the eigenvalue coefficients are α = −β = −ε/ ε(ε + 2λ), γ = −1 and one can obtain simple analytical expressions for all the transport coefficients. In particular the equilibrium averages (40)(41)(42) read as The scattering kernel within the first Born approximation reads as (see Appendix C for details) One sees that the scattering kernel is left-right symmetric (W (ϑ) = W (−ϑ)) to all orders in the scattering potential, i.e. skew-scattering is absent and no extrinsic spin Hall effect (SHE) can be expected in this case. Interestingly, the absence of any form of SHE is an exact property of the model provided that the impurities are non-magnetic. This is a consequence of the in-plane spin texture of the minimal 2D Dirac-Rashba model and can also be understood from the exact Ward identities for the 4-point vertex function [50]. The solution of the transport equation (43) reads By using the explicit expressions of the transport time (50) and of the equilibrium average (45) one obtains By inserting this last result in the Eq.(39) and considering the equilibrium average (46) one gets the following current-induced spin polarization while the charge current reads These relations imply The current-induced spin polarization is in plane and orthogonal to the applied electric field as a consequence of the symmetries of the model. In fact, it is easy to see that the inversion symmetry in the plane of graphene implies S z = 0 and S ·E = 0. This is also the case in the presence of sublattice-staggered interactions (i.e. {λ sv , ∆} = 0) due to the existence of a mirror reflection symmetry in the plane of the heterostructure [63]. The locking of nonequilibrium spin polarization density and charge current in plane at 90 degrees Eq. (54) is thus a general property of non-magnetic graphene heterostructures with C 6v or C 3v point group symmetry. These restrictions are lifted in magnetic graphene heterostructures, where S acquires collinear and out of plane components [63]. We note that Eqs.(52)-(53) coincide exactly with the results for the electric-field-induced spin polarization [76] obtained via the Kubo-Streda linear response theory by Offidani et al [18] and confirm the equivalence of the two approaches.

VI. THE SKEW-SCATTERING MECHANISM
We now consider the self-energy expansion beyond the first Born approximation. This is relevant for models with non-coplanar spin texture (λ sv = 0 and λ = 0), for which skew scattering mechanism is active and thus the model supports an extrinsic SHE with a semiclassical scaling J z i ∝ ij τ E j , in addition to the intrinsic SHE driven by Berry curvature effects [25,50,78]. The extrinsic mechanism is expected to provide the dominant contribution to the SHE in ultra-clean devices with high charge carrier mobility (i.e. τ 1). When the scattering potential is not too strong (|G This approximation corresponds to keeping the first two diagrams in the second line of Fig. 2. We stress that the real part in Eq. (55) does not give contribution to the skew scattering mechanism, but only renormalises the lowest-order scattering amplitude. Mathematically, the skewness in the effective scattering kernel W (ϑ) = W (−ϑ) results from the imaginary part of the retarded/advanced T -matrix, which endows the scattering term in Eq. (43) with an asymmetric contribution, where the magnitude of the effect is controlled by the "coupling constant" The commutator under the trace implies that W ss is odd upon the interchange of θ and θ (in particular, W ss vanishes for θ = θ ). A shift θ + ψ, θ + ψ clearly leaves W ss unchanged because of the periodicity of the projectors with respect to both angles, and hence there can be no dependence on the sum θ + θ . As a result, W ss must be an odd periodic function of ϑ = θ − θ . The solution of the transport equation generalizes Eq.(49) and reads as where are transport times defined through the microscopic rates These rates generalize the well-known transport rates in semiclassical transport theory to include the finite asymmetric rates generated by spin-orbit scattering. The expression is formally equivalent to the solution of the linearized Boltzmann transport equations for 2D fermions with SOC [25,48]. In Appendix D we derive the expression of the microscopic rates τ −1 s and τ −1 a in terms of the parameters α, β and γ defining the generic eigenvector (6) The zero-temperature nonequilibrium averages to be discussed below are computed according to Eq.(39), taking into account the solution (59). After the integration over ω one obtains As a result the nonequilibrium average of a generic observable reads Finally, from the expression Eq.(66) for the physical observables, the Drude conductivity (σ xx = σ yy = J x /E x ), the spin galvanic susceptibility (χ yx = −χ xy = S y /E x ) and the spin Hall conductivity (σ sH ≡ σ z yx = −σ z xy = J z sy /E x ) read respectively : with N F given in Eq.(10) and From the above equations, one sees that σ xx , (ev/ )χ yx and (e/ )σ z yx are expressed in units of conductance quantum. With this choice, the expressions (67-69) depend only on dimensionless combinations of the various parameters. Following Ref. [18], we quantify the charge-to-spin conversion efficiency (CSC) with the the following figure of merit (v F = (1/2π)(k F /N F )) For the SHE, we introduce similarly the spin Hall angle In the numerical analysis to be presented below, we will express all energies in units of λ. The disorder enters through the standard combination ετ / and the coupling constant g ss defined in Eq. (58). and τ ⊥ , which can be obtained by inserting Eqs. (7)(8)(9) in Eqs. (63)(64), is too cumbersome to be presented here, although their numerical evaluation is straightforward.
The skew-scattering rate τ −1 a shows a characteristic sign change upon increasing the Fermi energy. Such a sign change can be recognized by looking at the analytical expression (64). When the Fermi energy matches exactly the BR energy scale, i.e. |ε| = |λ|, the coefficients of the eigenvector (6) acquire a simple expression which implies the vanishing of the factor 1−γ 2 , in the formula (64) for the microscopic rate τ −1 a . Hence, remarkably, the sign change occurs at a well defined energy when the structure of the eigenvector implies the vanishing of the skew scattering already in lowest order for the effective scattering amplitude in the eigenstates (cf. Eq.(C2) and the discussion in Appendix C). In this respect it is illuminating to notice the following equilibrium average i.e. the factor controlling the sign change, is the equilibrium average of the z-axis component of the total internal angular momentum, which includes both spin and pseudo spin. A change of sign as function of the Fermi energy was also observed for the anomalous Hall effect of Dirac fermions in [25], where it was interpreted in terms of the behavior of the equilibrium spin texture of eigenstates. Figure 3 shows the Fermi energy dependence of the CSC efficiency. It is apparent that the figure of merit increases monotonically with the carrier density within the spin-helical regime, attaining a maximum at |ε| = ε + . In regime II, the CSC efficiency decreases monotonically with the Fermi energy |ε| [18] due to the presence of two bands with opposite spin texture. Hence the upper edge ε = ±ε + (for electrons/holes) sets the Fermi energy at which the inverse spin galvanic effect is the most efficient. Upon increasing the spin-valley coupling, there is an interesting effect depending on the charge carrier density. Near the edge of the spin-minority band (ε + ), the spinvalley coupling reduces the CSC efficiency, with an opposite effect observed near the lower edge where the "Mexican hat" feature develops (see Fig. 1 (c)). Such a behavior correlates well with the two different-sign regimes for the SHA (see Fig. 4). The latter changes from positive to negative upon increasing the Fermi energy |ε|, which is a consequence of the sign change in the microscopic scattering rate τ −1 a previously discussed and occurs when the Fermi energy matches the Rashba SOC.
In the spin precession measurements reported in Ref. [45], the spin galvanic nonlocal resistance (R SGE ) in a graphene/WS 2 lateral device indeed shows a maximum, in absolute value, as function of the back-gate voltage (V g ) for both charge carrier polarities. Considering that the charge-neutrality point in the experiment is approximately located at V g ≈ −10 V and the nonlocal signal maximum (for electrons) at V g ≈ −5 V, one may estimate a charge carrier density n ≈ 3.5 × 10 11 cm −2 at the upper edge of the spin-helical regime (see, for example, Ref. [79]). Using the relation between electronic density and Fermi wavevector for regime I (n = k 2 F /2π), we find ε + ≈ 60 meV. By reasonably assuming λ sv λ (in accord with a recent study [67]), one can further estimate λ ≈ 29 meV. This is a somewhat suprising result, given that functional theory calculations predict a proximityinduced BR SOC on the order of 0.4 meV [27,67]. In contrast, our estimate compares well with recent quantum Hall effect measurements indicating λ ≈ 15 meV in graphene/WSe 2 [20]. Strictly speaking, our estimate should be understood as an upper bound to the true BR SOC energy scale since the gate voltage dependence of R SGE also reflects the spin diffusion within the channel, which tends to shift the nonlocal signal maxima towards higher gate voltages |V g | (see also Ref. [47]). Most importantly, the very evidence of a maximum in |R SGE | as a function of carrier density shows that the spin-helical regime is experimentally accessible. Moreover, the qualitative behavior of the CSC efficiency (unlike the spin Hall angle [25]) is little sensitive to microscopic details of the impurity potential [18]. This provides extra confidence that the behavior of the spin galvanic nonlocal signal in Ref. [45] is a fingerprint of emergent spin-helical Dirac fermions at the graphene/TMD interface. On the other hand, a sign change in the SHE signal as the back-gate voltage is swept at fixed carrier polarity as predicted here (see Fig. 4) is not evident in the experimental data [45]. The significant uncertainty in the nonlocal signal at low temperatures renders a quantitative comparison between theory and experiment more challenging. According to our theory, the sign change occurs when |ε| = |λ| ≈ 29 meV. Incidentally, the latter roughly coincides with the energy scale of electron-hole puddles in the experiment δE = v √ πδn (with δn ≈ 2.5 × 10 11 cm −2 the residual carrier density), which could smear out the features of the extrinsic SHE predicted here. Therefore, the sign change in the SHA, if observed in cleaner samples, could provide a direct measure of the BR SOC energy scale. Both features, i.e. the maximum in the CSC efficiency and the sign change in SHA, thus provide valuable transport fingerprints of the spin-helical regime realized in graphene with proximity-induced SOC.
We briefly comment on the validity of our microscopic formulation. Close to the charge neutrality point, a full self-consistent evaluation of the self-energy is required to account for the behavior of the longitudinal charge conductivity. Because the weak-disorder condition implies τ 10, which requires not too low carrier density, one needs also a very large SOC to have the regime I spanning a reasonable energy range. Notwithstanding, the predicted sign change in the SHA at |ε| ≈ |λ| should be generally valid irrespective of the level at which the selfenergy is evaluated and may provide and experimental test of the theory. A detailed self-consistent evaluation to extend our theory closer to the charge neutrality point remains an important development for future work, but is beyond the scope of the present paper.

VII. CONCLUSIONS
In conclusion, we developed a microscopic theory that captures the key features of coupled spin-charge transport in graphene-based heterostructures, in which Dirac fermions experience proximity-induced spin-valley coupling as well as Bychkov-Rashba effect due to interfacial breaking of inversion symmetry. We have restricted ourselves to the spin-helical regime realized at low electronic densities, characterized by the locking of spin and pseudospin degrees of freedom. Such a simplification has allowed us to derive, within the framework of the Eilenberger equation for the quasiclassical Green's function, a single transport equation capturing the low-energy behavior of dirty spin-orbit-coupled Dirac fermions. The Eilenberger equation exhibits a scattering kernel, which is derived within a T-matrix expansion by projecting the disorder potential in the energy eigenstate at the Fermi energy. Both symmetric and skew scattering features have been connected to the spin and pseudo spin texture of the eigenstate through explicit analytical expressions. The spin-charge transport response functions are described in terms of longitudinal (τ ) and transverse (τ ⊥ ) transport times. The transverse transport time is different from zero only in the presence of both Bychkov-Rashba interaction and spin-valley coupling, in accord with the exact covariant conservation laws of the Dirac-Rashba model [50]. Interestingly, the SHE vanishes, then changing sign, when the Fermi energy crosses the Bychkov-Rashba energy scale, irrespective of the value of the spin-valley coupling. This is a consequence of the fact that the equilibrium spin and pseudo spin texture, at this particular value of the Fermi energy, implies a vanishing total spin angular momentum. In the spin-helical energy range, the inverse spin galvanic effect has a maximum at the edge of the spin-minority band. These features then represent the fingerprints of the spin-helical regime and may provide a direct way to estimate the proximity-induced SOC in graphene heterostructures using spin precession measurements in nonlocal geometry that can disentangle SHE and inverse spin galvanic signals [4]. The comparison with recent experimental results [45] provides a reasonable estimate for the Bychkov-Rashba SOC parameter. More importantly, the energy range within which the theory is valid appears to be experimentally accessible, so that the results presented in this paper can be put to a test. R.R. and A.F. are grateful to Sergio O. Valenzuela for useful comments. R.R. would like to thank Cosimo Gorini for useful discussions. A.F. gratefully acknowledges the financial support from the Royal Society, London through a Royal Society University Research Fellowship.
Although the disorder potential does not have a dependence on the scattering angle, its projection on the energy eigenstates gives rise to such a dependence. To this end we define the scattering amplitude for eigenstates at the Fermi surface with wave vector k = k F n and k = k F n where n = (cos θ, sin θ) and similarly for n p(θ, θ ) ≡ Φ(k F , θ)|Φ(k F , θ ) .
Appendix D: Expression of the effective potential and transport times in the C3v-invariant model We want to find the expressions of the scattering kernel in the general case in terms of the eigenvector coefficients defined in Eq. (6). The Born approximation and skew scattering contributions are defined in Eqs. (44) and (57), respectively. For the scattering kernel in the Born approximation we need, apart from the eigenvector normalization factors, Tr (P (θ)P (θ )) = p(θ, θ )p(θ , θ).