Theory of spin-charge-coupled transport in proximitized graphene: An SO(5) algebraic approach

Establishing the conditions under which orbital, spin and lattice-pseudospin degrees of freedom are mutually coupled in realistic nonequilibrium conditions is a major goal in the emergent field of graphene spintronics. Here, we use linear-response theory to obtain a unified microscopic description of spin dynamics and coupled spin-charge transport in graphene with an interface-induced Bychkov-Rashba effect. Our method makes use of an SO(5) extension of the familiar inverse-diffuson approach to obtain a quantum kinetic equation for the single-particle density matrix that treats spin and pseudospin on equal footing and is valid for arbitrary external perturbations. As an application of the formalism, we derive a complete set of drift-diffusion equations for proximitized graphene with scalar impurities in the presence of electric and spin-injection fields which vary slowly in space and time. Our approach is amenable to a wide variety of generalizations, including the study of coupled spin-charge dynamics in layered materials with strong spin-valley coupling and spin-orbit torques in van der Waals heterostructures.


Introduction
There is a current fundamental and technological interest in the harnessing of spin-orbit-coupling (SOC) effects in nonmagnetic media, particularly for the interconversion of charge and spin currents and the generation of nonequilibrium spin polarization [1,2]. A recent trend is the use of two-dimensional (2D) materials to engender electrical control over SOC effects benefiting from their reduced dimensionality and unique opto-electronic properties afforded by atomically thin crystals and their heterostructures [3,4]. With graphene well established as a high-fidelity spin channel material supporting room-temperature spin transport over lengths up to tens of micrometers [5,6,7,8,9,10], an important challenge concerns the manipulation of nonequilibrium spins by pure electrical means for future spinlogic applications. While at first glance this might seem hopeless given the absence of bulk ferromagnetism in graphene [11], not to mention its ultra-low intrinsic SOC [12], several strategies have been proposed to overcome this bottle neck, including adatom engineering [13,14] and proximity effects achieved via van der Waals coupling to high-SOC 2D materials [15,16,17,18,19]. The latter approach has shown great promise because the proximityinduced SOC can be well resolved in energy (i.e., on the order of the quasiparticle broadening), which facilitates the experimental demonstration of SOC effects with reproducibility, e.g. by means of low-field magnetotransport measurements [20,21,22,23,24]. Furthermore, the strong interplay of spin and lattice-pseudospin degrees of freedom in honeycomb layered materials gives rise to fingerprints of unique hallmarks of SOC in spin transport experiments. Most notably, the emergence of spin-helical 2D Dirac fermions in van der Waals heterostructures due to the Bychkov-Rashba (BR) effect [25] has been predicted [26] and demonstrated experimentally [27,28,29,30] to enable efficient spin-charge interconversion at room temperature. Owing to a unique spin-pseudospin entanglement of electronic wavefunctions, the sign and magnitude of the nonequilibrium spin polarization in BR-coupled graphene can be tuned with a back-gate voltage [26], in contrast to spin-galvanic effects generated by topological insulators [31]. Another interesting manifestation of proximity-induced SOC is observed in Hanle-type spin precession experiments [32,33], where a rather unconventional spin dynamics results in spins polarized in the plane of graphene relaxing about ten times faster than out-of-plane spins. This effect, originally predicted by Cummings and co-workers [34], is explained by the combined action of interface-induced spin-valley coupling and intervalley scattering triggered by point defects, which opens an additional Dyakonov-Perel-like relaxation channel for in-plane spins.
Due to the fast experimental progress in the field, the theory is generally lagging but there are notable exceptions.
A microscopic theory of spin dynamics for graphene-based van der Waals heterostructures with C 3v point-group symmetry was put forward in Ref. [35]. Moreover, a controlled diagrammatic approach to calculating linear response functions in the presence of disorder and generic proximity effects was developed in a recent series of works [26,36,37,38,39]. An electrical detection scheme that enables full disentanglement of competing SOC transport effects in diffusive lateral spin-valve devices was also proposed recently in Ref. [40]. These early works highlighted the key role played by symmetry, quantum geometry and impurity scattering in the spin dynamics and spin-orbit-coupled transport phenomena, such as the spin Hall effect (SHE) [41], in honeycomb layers. The diagrammatic approach has proven ideal for the study of the weak-disorder limit relevant for clean samples with large mean free paths, where numerical simulations have traditionally struggled [42,43]. Interestingly, the emergence of noncoplanar k-space spin textures in broken inversion symmetry conditions was shown to allow for a robust skew-scattering-induced SHE that dominates over the intrinsic (spin-Berry-curvature) contribution in the clean limit, while not requiring spin-dependent disorder potentials [36]. On the other hand, tight-binding methods have provided useful insights into the highly disordered limit (e.g. via simulations of extrinsic SHE efficiency in samples with a high coverage of adatoms [44]). However, a unified semiclassical description of spin-orbit-coupled transport effects in the presence of generic time-and space-dependent perturbations is still lacking, even for the simplest case of C 6v -symmetric graphene heterostructures. The aim of this paper is to fill this gap by developing a theoretical framework that encompasses all the known phenomenology and has the potential to provide new predictions for the numerous opto-spintronic phenomena supported by van der Waals materials [40,42,45,46,47,48]. To that end, we devise a Green's function inverse diffuson approach that respects the SO(5) algebraic structure of spinorbit-coupled 2D Dirac fermions (rather than projecting out the sublattice degree of freedom to obtain a simplified description in terms of Bloch-type equations [49]), which can be used to derive quantum kinetic equations in a simple and mathematically transparent fashion. This extension of the powerful diffuson approach originally developed for 2D electron gases [50,51] will allow us not only to keep track of the entangled dynamics of pseudospin and spin observables, but also to obtain useful analytical expressions for linear response functions to generalized fields of experimental relevance, such as Zeeman-type spin-injection fields.
This paper is organized as follows. In Sec. 1.1, we introduce the 2D Dirac-Rashba model capturing the lowenergy dynamics of C 6v -symmetric graphene heterostructures. Sec. 1.2 describes the diagrammatic framework and the semiclassical approximation used in treating scattering effects. In Sec. 2 we derive the semiclassical driftdiffusion transport equations for the experimentally accessible macroscopic observables and discuss various features and applications of the formalism. Sec. 3 presents our conclusions.

The 2D Dirac-Rashba model
For our theoretical study of spin-charge-coupled transport in spin-orbit-coupled graphene, we use a model of noninteracting electrons subject to a random impurity potential. Because we are interested in emergent phenomena stemming from interfacial breaking of inversion symmetry, we assume that the BR interaction is the dominant SOC effect (see Figs. 1(a) and (b)) [25,52]. The effective Hamiltonian at the K valley can be written in terms of tensor products of Pauli matrices σ a ⊗ s b (a, b = x, y, z) acting on the pseudospin-spin space as [53] where v 10 6 m/s is the bare Fermi velocity of the massless Dirac fermions, p µ = (−ε/v, −ı ∇) is the energymomentum 3-vector, A µ = a=x,y,z A µ a s a (µ = 0, x, y, z) is the non-Abelian SU(2) gauge field capturing all spin-dependent effects and V (x) is the impurity potential. Here, summation over repeated indices is implied with σ 0 denoting the identity matrix. The BR effect with coupling strength λ is encoded in the gauge field components A x y = −A y x = λ/v, which generate the corresponding spin-orbit interaction via minimal coupling [H BR = λ (σ×s) z ]. Equation (1) has been dubbed the 2D Dirac-Rashba model in recent literature [26,36]. We are not concerned with purely extrinsic transport phenomena induced by spin-orbit-active impurities (i.e. random SOC), which have been the object of detailed microscopic analysis in early work [13,54,55,56,57,58]. Rather, our focus here is on the transport properties of spin-helical 2D Dirac fermions realized in graphene-based heterostructures with a dilute The dispersion relation of the minimal Dirac-Rashba model, H 0k = v σ · k + H BR , reads as where µ(ν) = ±1 labels, respectively, the spin-helicity and polarity of charge carriers ( Fig. 1(b)). The BR effect lifts the spin degeneracy and locks the spin polarization and wavevector at right angles, leading to a spin-helical configuration in k-space ( Fig. 1(c)). The Dirac nature of the carriers enables a low electronic density regime characterized by a simply-connected Fermi surface with well defined spin-helicity (i.e. µ = −1), akin to the surface states of topological insulators, with interesting consequences for spin-charge conversion effects [35].
The equilibrium spin polarization associated with the Bloch eigenstates of H 0k is easily computed as where σ µνk = (1/ v)∇ k µν (k) is the expectation value of the pseudospin polarization vector. The spin winding of the Fermi surface is shared with other surfaces possessing broken inversion symmetry, but unlike conventional spin-helical states, its energy dependence reflects a spin-angular momentum transfer between spin and pseudospin channels. Such a spin-pseudospin coupling is responsible for the conspicuous wavevector dependence of the equilibrium spin texture [Eq. (3)]. Indeed, the spin-polarization magnitude is vanishing at the corners of the first Brillouin zone, where the band velocity is vanishing ( σ µνk→0 = 0). Away from the zone corners K and K , the spin texture magnitude increases monotonically with the Fermi wavevector until it saturates away from the spin gap (i.e. for |ε| 2|λ|). As noted by Rashba [25], the strong momentum-dependence of the spin texture at low energies is a unique fingerprint of BR-coupled 2D Dirac fermions.
The random potential in Eq. (1), which in this work is assumed to be scalar, affects the spin dynamics of charge carriers by inducing scattering between electronic states with different effective Larmor fields Ω µνk = λ s µνk ≈ −µνλk ×ẑ for |ε| |λ|. Due to the random change of precession axis, an initial nonequilibrium spin polarization will decay exponentially with time. In the standard weak-SOC regime (realized in systems with a small spin splitting compared to the disorder-induced broadening [59,60,61]), the spin-relaxation rate is τ −1 s ∝ λ 2 τ , where τ is the elastic scattering time ( Fig. 1(d)). Interestingly, the spin-relaxation rate of out-of-plane spins (i.e. polarized along theẑ-axis) is twice that of in-plane spins (τ s,⊥ /τ s, = 1/2), which could provide an experimentally detectable signature of BR effect using Hanle-type spin precession measurements [62,63]. Impurity scattering also plays a key role in coupled spin-charge transport phenomena, profoundly affecting the efficiency of spin Hall and spin-charge conversion (spin-galvanic) effects, even in the clean limit with |ε|τ 1, as we shall see below.

Theoretical framework
To derive a rigorous microscopic picture of coupled spin-charge transport, we evaluate the density matrix response function employing many-body perturbation theory methods [64,65,66]. Our aim is to generalize the familiar diffuson approach for 2D electron gases [50,51] to accommodate the enlarged SO(5) Clifford algebra of 2D Dirac fermions. Because we are interested in the diffusive regime realized in weakly disordered samples with |ε|τ 1, we neglect quantum corrections arising from weak localization and higher-order spin-orbit scattering effects, such as quantum side jumps and diffractive skew scattering [57,58,67]. In the diagrammatic language, such semiclassical approximation amounts to discarding crossing diagrams encoding coherent multiple impurity scattering events, as well as higher-order noncrossing diagrams with a perturbation parameter 1/(|ε|τ ) 1 [57]. Unless stated otherwise, we work in natural units with ≡ 1 ≡ e. Additionaly, for ease of notation, we assume throughout that ε, λ > 0.
The central object in our approach is the real-time retarded(R)/advanced(A) single-particle Green's function where G a 0k (ε) = [ε − vσ · k − λ (σ × s) ·ẑ + ı a 0 + ] −1 is the Fourier transform of the clean Green's function and is the quasiparticle self energy in the noncrossing approximation. For uncorrelated short-range impurity potentials, the self energy is k-independent, and so we set Σ a k (ε) ≡ Σ a (ε) from here onwards. In order to keep our discussion as simple as possible, we compute the self energy under the assumption of weak Gaussian disorder. This will simplify our analytical treatment significantly while capturing the essential physics [68]. Diagrammatically, the weak disorder approximation amounts to retaining only the contribution from the areal density and u 0 is the potential strength, one arrives at where 1/2τ ≡ η = n i u 2 0 ε/4v 2 is the quasiparticle broadening and γ r ≡ (σ x s y − σ y s x )/2. The existence of two distinct transport regimes at low energies (i.e. ε > 2λ and ε < 2λ) is a unique feature of the 2D Dirac-Rashba model, which is responsible for the existence of a maximum in current-induced spin polarization efficiency when the (gate-tunable) Fermi energy lies precisely at the spin-gap edge [26].
In this paper, we are primarily interested in the diffusive coupled spin-charge dynamics which occur in graphene flakes with weak proximity-induced SOC at moderate-high charge carrier densities (ετ 1 λτ ) [37]. The condition ε λ implies that the BR-slit bands with opposite spin helicities are occupied at the Fermi level. To study the behavior of the system away from equilibrium, it is convenient to introduce the generalized one-particle density operatorρ where γ αβ = σ α ⊗ s β (with α, β = 0, x, y, z) span the vector space of Hermitian 4 × 4 matrices. The density matrix characterizing a given nonequilibrium state can be expanded as a linear combination of the Clifford algebra basis elements, such that where ... denotes quantum and disorder averages and d = dim H ≡ 4 is a normalization factor. The expectation value of a generic local observable, where tr indicates the trace over internal degrees of freedom and ρ αβ (x, t) ≡ ρ αβ (x, t) .
In this work we are concerned with the semiclassical dynamics of typical spin transport observables, such as the spin-polarization density and the spin current density. Thus, it is more convenient to work directly with the deviation from equilibrium of the expectation values, i.e. δO(x, t) Here, ρ 0 αβ denotes the equilibrium part of the (disorder-averaged) density matrix. The macroscopic observables of interest to us are the charge density, N , spin polarization density, S a (a = x, y, z), charge current density, J i (i = x, y), and spin current density, J a i . The corresponding expectation values away from equilibrium are defined as where we have temporarily reinstated and e to distinguish between charge and spin currents.
The interaction Hamiltonian is where A ext αβ (x, t) (α, β = 0, x, y, z) spans the Clifford algebra and thus describes any type of charge-spin perturbation applied to the system. In Sec. 2.1, we shall show that the linear response density matrix δρ(x, t), when properly coarse-grained over typical length and time scales (i.e., |x| l ≡ vτ and t τ ), is governed by an enlarged 16 × 16 diffuson Hamiltonian that is the SO(5) analogue of the familiar inverse density fluctuation propagator of 2D electron gases [50,51,69] and topological insulators [70,71]. In Sec. 2.2, the linear response machinery will be applied to derive the full set of drift-diffusion transport equations for the variables δN , δS a , δJ i and δJ a i , and thus establish a rigorous microscopic picture for the coarse-grained dynamics of the problem. In the following, we define N ≡ δN , S a ≡ δS a , J i ≡ δJ i and J a i ≡ δJ a i for ease of notation.

Diffuson Dirac Hamiltonian
We start by setting up the formalism needed to derive a quantum kinetic equation for BR-coupled 2D Dirac fermions.
Let us briefly discuss the general structure of H D (q, ω) in the long wavelength limit of interest to us. Zero entries in the 16 × 16 bubble matrix M (0, ω) can be readily identified by applying the following C 6v point-group operations [72]: (i) C 2 rotation exchanging sublattices and (ii) mirror-reflection R x leaving sublattices invariant.
We now turn to the spin-charge eigenmodes sustained by the system. The first step is to compute the gradient expansion of Eq. (21). Working in the diffusive regime ( τ 1 λτ ) greatly simplifies matters due to many bubbles being parametrically small. To leading order in v|q| and ωτ , we find where Q, P i and L are 16 × 16 matrices given by Q = −ı(∇ ω H D (0, ω)) ω=0 , P = ı(∂ q H D (q, 0)) q=0 and L = H D (0, 0). The calculation of these matrices is rather cumbersome, yielding unwieldy expressions for the matrices Q = Q(ε, λ, τ ), P i = P i (ε, λ, τ ) and L = L(ε, λ, τ ). We therefore provide the explicit form of the Eq. (23) in the Appendix, from which Q, P i and L can be inferred.
The diffuson Hamiltonian in Eq. (23) is linear in the wavevector q due to the Dirac nature of the lowlying excitations in graphene heterostructures. The diffuson spectrum at low energies is shown in Fig. 3. The quadratic dispersion of the gapless mode follows from the diffusive pole structure of the density-density response, i.e. [R 00,00 (q, 0)] −1 ∼ (vτ |q|) 2 . For q = 0, this mode is an admixture of charge N , spin current J a i and spin polarization S x,y fluctuations. The gapped states, on the other hand, describe eigenmodes of the nonequilibrium spin polarization S a [73]. Note that two of these states (S x,y ) are degenerate at q = 0 due to the rotational (C v∞ ) symmetry of the single-particle Hamiltonian [Eq. (1)]. Their gaps at q = 0 are given by ∆ s ≡ ∆ (λ/η) 2 for the S x,y -and ∆ ⊥ 2∆ for the S z -mode. The twice as fast dephasing of out-of-plane spin fluctuations is a fingerprint of BR SOC, a feature which is known to survive even in the strong (unitary) scattering regime [37].
The low-lying modes shown in Fig. 3 are remarkably similar to those of a BR-coupled 2D electron gas [51].
We verified that the diffuson spectra of the two Rashba models can be exactly mapped onto each other, despite A comment is in order regarding the validity of our assumptions in the light of recent findings. In the example of Fig. 3, we considered a small BR coupling of only 0.1 meV, which is in line with density functional theory calculations for clean graphene/group VI dichalcogenide heterostructures [16,18]. On the other hand, the semi-empirical Slater-Koster parametrization of Ref. [74], as well as early magnetotransport measurements [20,21,22,23,24], have found much higher proximity-induced SOC (up to ≈ 10 meV). Furthermore, a recent joint theory-experiment study suggests that the interfacial Rashba coupling can be made as large as 100 meV, by placing a graphene flake on top of a metallic substrate with a suitable work function mismatch [29]. The theory developed in this work is expected to remain accurate provided that the Dirac bands remain intact (so that the low-energy picture in Fig. 1 is justified) and the system is sufficiently disordered so that λτ 1.

Unified coupled spin-charge drift-diffusion equations
The quantum kinetic equation governing the one-particle reduced density matrix in the large distance and long time limits is obtained after an inverse Fourier transform (−ıω → ∂ t and ıq i → ∇ i ) of Eq. (23) as Transport equations for the coarse-grained variables N ( can now be derived by replacing the 16 × 16 matrices Q = Q(ε, λ, τ ), P i = P i (ε, λ, τ ) and L = L(ε, λ, τ ) with their explicit forms (see Appendix). To leading order in 1/ετ and λτ , we obtain where E i = −∂ t A ext i0 is the externally applied electric field and B a = (1/v) ∂ t A ext 0a (x, t) is the external Zeeman field ("spin injection field") that induces a nonequilibrium spin density [75]. The coefficients Λ ab c , Γ ab ic , Ω a ic and Υ ab ic are listed in Table 2 and ij denotes the rank-2 Levi-Civita symbol (i, j = x, y). Different perturbations (e.g. a spin-dependent electric field E a i = −∂ t A ia ) can be easily incorporated via a suitable parameterization of the generalized external vector potential entering Eq. (24).
The drift-diffusion equations (25)- (28) are the main results of this work. Equations (25) and (26) where D = v 2 τ is the diffusion constant [76], hold to good accuracy insofar as ω λ. Sublattice-staggered charge N s and spin densities S a s (see Table 1) are conspicuously absent from these relations. In terms of the underlying diffuson Hamiltonian, these observables are linked to dispersionless modes with large gaps and are thus effectively decoupled from the low-energy dynamics. We expect such terms to play a role in graphene heterostructures with broken sublattice symmetry [34,35,36], which is beyond the scope of the this article.

Spin Hall and spin-galvanic effects: the DC regime
Equations (25)-(28) provide a physically transparent scheme to interpret and predict a variety of SOC phenomena of fundamental and technological relevance. Before discussing new applications of the formalism, we briefly revisit two well established results for BR-coupled graphene [26,36]. We start with the SHE [41], i.e. the appearance of a transverse spin current upon application of a DC charge current, first observed in semiconductors [77,78]. In graphene with random SOC (e.g. induced by dilute adatoms), a robust SHE can be induced via resonant skew scattering [13,14,54,55,56]. In contrast, for graphene systems with a spatially uniform BR effect, the SHE is strictly vanishing unless supplemented with proximity-induced spin-valley coupling or spin-dependent disorder.
Detailed information on the SHE can be obtained from the drift-diffusion equations with little effort. The so-called intrinsic spin Hall angle θ int sH = 2τ λ 2 /ε, which appears explicitly in Eqs. (29)-(30), diverges in the clean limit (τ → ∞) and does not correspond per se to a steady-state transport quantity. The actual spin Hall angle, defined as the ratio between near-equilibrium spin Hall current and applied charge current, is obtained by solving the system of coupled equations (25)- (28) and receives important disorder corrections even in the clean limit. The "SHE cancellation" in the DC limit θ dis sH (0) = −θ int sH is a fundamental consequence of SU (2)  spin covariance of pure Rashba models as shown by Dimitrova [79]. This result can also be interpreted as the unavoidable outcome for a 2D system with an isotropic and fully in-plane spin texture. Because the electronic states are admixtures of orthogonal spin states, phase shifts experienced by the spin-up and spin-down components of scattered wavefunctions from scalar impurities cannot be distinguished, implying the absence of skew scattering [36]. As discussed in Sec. 2.4, a robust SHE nevertheless takes place at finite frequencies (i.e. an optical SHE) or when the system is perturbed by a spin-injection field.
Also of interest is the inverse spin-galvanic effect (ISGE), whereby an applied current magnetizes the conduction electrons, thus generating a net spin polarization density [80,81]. Its microscopic origin lies in the spin-momentum locking of Bloch eigenstates caused by the BR effect [ Fig. 1 (c)]. While the equilibrium spin polarization averaged over the Fermi surface in a nonmagnetic system must vanish identically [Eq. (3)], an external electric field effectively breaks the time-reversal symmetry, by causing an imbalance the occupation of states with opposite momenta, which allows the build up of a net transverse spin polarization. The ISGE efficiency can be easily read out from Eq. (30) by replacing the pure spin current by its steady-state value in the minimal model, i.e. J a i = 0. The charge-to-spin conversion efficiency is obtained as: This relation (first derived in Ref. [26]) discloses an optimal spin-charge conversion efficiency at the spin-gap edge, i.e. κ xy (ε = 2λ) = 1. The ISGE efficiency parameter decays algebraically with the energy of charge carriers, which makes the effect detectable at room temperature over a wide range of charge carrier densities. The robust ISGE in graphene with BR effect has been observed in a recent series of experiments [27,29,28,30].

Application: Optical spin Hall and spin galvanic effects
As a novel application of the formalism, we derive the optical response of BR-coupled graphene. To this end, we solve the charge-spin drift-diffusion equations [Eqs. (26)-(28)] subject to time-dependent electric and spin-injection fields, with Fourier transforms E i (ω) (i = x, y) and B a (ω) (a = x, y, z), respectively. In the Dyakonov-Perel regime with ωτ λτ 1 τ , the macroscopic observables of interest are found, after tedious but straightforward calculations, to be where the factor of g v = 2 accounts for valley degeneracy. Furthermore, zij denotes the rank-3 Levi-Civita symbol, where subleading corrections of order ωτ have been neglected for simplicity. We note that the expression for σ xx (ω) coincides with the familiar Drude model result, which is valid in the semiclassical regime with ετ 1 [82]. More interestingly, the optical spin-galvanic susceptibility [Eq. (38)] generalizes the findings of Ref. [26] to an external electric field with nonzero frequency. Here, ωτ emerges as an important parameter that governs the imaginary part of the response function, whereas the physics of the DC regime reflects the average spin-precession angle experienced between consecutive scattering events (i.e. θ p = λτ ).

Summary and outlook
We derived a quantum kinetic equation and the associated set of coupled spin-charge linear transport equations that govern the dynamics of Rashba-coupled 2D Dirac fermions in graphene proximitized by high-SOC materials.
These equations, which are valid in the presence of arbitrary external fields, provide a quantitative description of rich interlinked spin-orbit scattering phenomena characteristic of 2D systems with broken inversion symmetry, including Dyakonov-Perel-type spin relaxation, direct and inverse DC spin-galvanic and optical spin Hall effects.
The distinctive feature of the SO(5) algebraic approach formulated in this work is that the exact large distance and long time behavior of the linear response one-particle density matrix, and thus also the expectation value of any local observable, is uniquely determined by a generalized inverse diffuson matrix (i.e. a diffuson Hamiltonian) that spans the full vector space of a 16-dimensional Clifford algebra. This is to be contrasted with the familiar SU (2) approach for two-dimensional electron gases [50,83], whose Fourier-space diffuson operators are 4 × 4 matrices restricted to the space of charge and spin-polarization densities. The enlarged (16 × 16) diffuson Hamiltonian derived here emphasizes the manifestation of entanglement between the spin and pseudospin (sublattice) degrees of freedom that is ubiquitous across van der Waals materials. Furthermore, it provides direct access to the time evolution of all thermodynamic macroscopic observables, including spin-current density, spin-polarization density and sublattice-staggered densities.
This work opens up a number of avenues that can be explored in the framework introduced here. These range from the exploration of nonequilibrium opto-spintronic phenomena in group-VI dichalcogenide monolayers and van der Waals heterostructures with sizable spin-valley coupling, to the role played by spin-orbit-active impurities in spin dynamics and spin-charge conversion effects. In particular, asymmetric scattering precession [56] and skew scattering effects can be systematically explored by means of a nonperturbative T -matrix ladder scheme that resums all single-impurity scattering diagrams [26,37,39]. Such an extension of our diagrammatic treatment would be helpful in understanding the emergent transport physics of 2D van der Waals materials with broken sublattice symmetry, where the presence of non-coplanar k-space spin textures at low energies is known to enable robust spin Hall effects irrespective of the type of impurities [36,38], in addition to strongly modifying the spin dynamics [34,35]. Moreover, our formalism could be employed to investigate how proximity-induced SOC affects the nonlocal resistance in Hanle-type spin precession experiments beyond its impact on the spin lifetimes. This could be examined by deriving a generalized spin diffusion equation accounting for the interplay of spin-valley coupling, intervalley scattering and the Bychkov-Rashba effect.

Data statement
This publication is theoretical work that does not require supporting research data.