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

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 figures 1(a) and (b)) [25, 52]. The effective Hamiltonian at the K valley can be written in terms of tensor products of Pauli matrices $\sigma_{a}\otimes s_{b}$ ($a,b = x,y,z$) acting on the pseudospin–spin space as [53]

$$\begin{equation} H = \int d^{2}\mathbf{x}\,\Psi^{\dagger}(\mathbf{x})\,\left[\,v\sigma_{\mu}\left(p^{\mu}+\boldsymbol{\mathcal{A}}^{\mu}\right)+V(\mathbf{x})\,\right]\Psi(\mathbf{x})\,, \end{equation} \tag{ 1 }$$

where $v\simeq10^{6}$ m s⁻¹ is the bare Fermi velocity of the massless Dirac fermions, $p^{\mu} = (-\varepsilon/v,-\imath\hbar\nabla)$ is the energy–momentum 3-vector, $\boldsymbol{\mathcal{A}}^{\mu} = \sum_{a = x,y,z}\mathcal{A}_{a}^{\mu}\,s_{a}$ ($\mu = 0,x,y,z)$ is the non-Abelian SU(2) gauge field capturing all spin-dependent effects and $V(\mathbf{x})$ is the impurity potential. Here, summation over repeated indices is implied with σ₀ denoting the identity matrix. The BR effect with coupling strength λ is encoded in the gauge field components $\mathcal{A}_{y}^{x} = -\mathcal{A}_{x}^{y} = \lambda/v$, which generate the corresponding spin–orbit interaction via minimal coupling [$\mathcal{H}_{\textrm{BR}} = \lambda\,(\boldsymbol{\sigma} \times \boldsymbol{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–58]. Rather, our focus here is on the transport properties of spin–helical 2D Dirac fermions realized in graphene-based heterostructures with a dilute concentration of spin-transparent impurities and a well-established (i.e. spatially uniform) BR effect.

Zoom In Zoom Out Reset image size

The dispersion relation of the minimal Dirac–Rashba model, $\mathcal{H}_{0\mathbf{k}} = \hbar v\,\boldsymbol{\sigma}\cdot\mathbf{k}+\mathcal{H}_{\textrm{BR}}$, reads as

$$\begin{equation} \varepsilon_{\mu\nu}(\mathbf{k}) = \mu\lambda+\nu\sqrt{\lambda^{2}+\hbar^{2}v^{2}k^{2}}\,, \end{equation} \tag{ 2 }$$

where $\mu(\nu) = \pm1$ labels, respectively, the spin-helicity and polarity of charge carriers (figure 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 (figure 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. $\mu = -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 $\mathcal{H}_{0\mathbf{k}}$ is easily computed as

$$\begin{equation} \langle\mathbf{S}\rangle_{\mu\nu\mathbf{k}} = \frac{\hbar}{2}\langle\boldsymbol{s}\rangle_{\mu\nu\mathbf{k}} = -\frac{\hbar}{2}\mu\,\langle\boldsymbol{\sigma}\rangle_{\mu\nu\mathbf{k}}\times\hat{\mathbf{z}}\,, \end{equation} \tag{ 3 }$$

where $\langle\boldsymbol{\sigma}\rangle_{\mu\nu\mathbf{k}} = (1/\hbar v)\nabla_{\mathbf{k}}\epsilon_{\mu\nu}(\mathbf{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 (equation (3)). Indeed, the spin-polarization magnitude is vanishing at the corners of the first Brillouin zone, where the band velocity is vanishing ($\langle\boldsymbol{\sigma}\rangle_{\mu\nu\mathbf{k}\rightarrow0} = 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 $|\varepsilon|\gtrsim 2|\lambda|$). 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 equation (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 $\boldsymbol{\Omega}_{\mu\nu\mathbf{k}} = \lambda\langle\mathbf{s}\rangle_{\mu\nu\mathbf{k}}\approx-\mu\nu\lambda\,\hat{\mathbf{k}}\times\hat{z}$ for $|\varepsilon|\gg|\lambda|$. 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–61]), the spin-relaxation rate is $\tau_{s}^{-1}\propto\lambda^{2}\tau$, where τ is the elastic scattering time (figure 1(d)). Interestingly, the spin-relaxation rate of out-of-plane spins (i.e. polarized along the $\hat{z}$-axis) is twice that of in-plane spins ($\tau_{s,\perp}/\tau_{s,\parallel} = 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 $|\varepsilon|\tau\gg1$, as we shall see below.

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–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 $|\varepsilon|\tau\gg1$, 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/(|\varepsilon|\tau)\ll1$ [57]. Unless stated otherwise, we work in natural units with $\hbar\equiv1\equiv e$. Additionaly, for ease of notation, we assume throughout that $\varepsilon,\lambda\gt0$.

The central object in our approach is the real-time retarded(R)/advanced(A) single-particle Green's function ($a = A,R\equiv-,+$) defined as

$$\begin{equation} G^{a}(\mathbf{x},\mathbf{x}^{^{\prime}};t-t^{^{\prime}}) = \mp\imath\left\langle 0|T\left[\Psi(\mathbf{x},t),\Psi^{\dagger}(\mathbf{x}^{^{\prime}},t^{^{\prime}})\right]|0\right\rangle \theta(\pm t\mp t^{^{\prime}}), \end{equation} \tag{ 4 }$$

where $\Psi^{\dagger}(\mathbf{x},t)$ [$\Psi(\mathbf{x},t)$] are 4-component field operators creating (annihilating) a 2D Dirac fermion at position x and time t, T is the time-ordering symbol and $\theta(.)$ is the Heaviside step function. After disorder averaging (indicated below by an overline), the Green's function in momentum-frequency space acquires the familiar form

$$\begin{equation} \mathcal{G}_{\mathbf{k}}^{a}(\varepsilon) = \frac{1}{[\mathcal{G}_{0\mathbf{k}}^{a}(\varepsilon)]^{-1}-\Sigma_{\mathbf{k}}^{a}(\varepsilon)}\,, \end{equation} \tag{ 5 }$$

where $\mathcal{G}_{0\mathbf{k}}^{a}(\varepsilon) = [\varepsilon-v\boldsymbol{\sigma}\cdot\mathbf{k}-\lambda\,\left(\boldsymbol{\sigma}\times\mathbf{s}\right)\cdot\hat{z}+\imath\,a\,0^{+}]^{-1}$ is the Fourier transform of the clean Green's function and

$$\begin{equation} \Sigma_{\mathbf{k}}^{a}(\varepsilon) = \int d\mathbf{(x-x^{^{\prime}})}\,e^{-\imath\mathbf{k}(\mathbf{x}-\mathbf{x}^{^{\prime}})}\overline{\langle\mathbf{x}^{^{\prime}}|V\frac{1}{1-G_{0}^{a}(\varepsilon)V}|\mathbf{x}\rangle}\, \end{equation} \tag{ 6 }$$

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 $\Sigma_{\mathbf{k}}^{a}(\varepsilon)\equiv\Sigma^{a}(\varepsilon)$ 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 'rainbow' diagram with two impurity potential lines; see figures 2(a) and (b). Replacing the explicit form of the two-point correlator $\overline{V(\mathbf{x})V(\mathbf{x}^{^{\prime}})} = n_{i}u_{0}^{2}\,\delta(\mathbf{x}-\mathbf{x}^{^{\prime}})$ in the Dyson expansion of the self energy, where n_i is the impurity areal density and u₀ is the potential strength, one arrives at

$$\begin{equation} \Sigma^{\pm}(\varepsilon) = \mp\,\left(\imath\,/2\tau\right)\left[\theta(\varepsilon-2\lambda)+\theta(2\lambda-\varepsilon)\left(1+\sigma_{z}s_{z}-\frac{\varepsilon}{\lambda}\gamma_{r}\right)\frac{\lambda}{\varepsilon}\right]\,, \end{equation} \tag{ 7 }$$

where $1/2\tau\equiv\eta = n_{i}\,u_{0}^{2}\varepsilon/4v^{2}$ is the quasiparticle broadening and $\gamma_{r}\equiv(\sigma_{x}s_{y}-\sigma_{y}s_{x})/2$. The existence of two distinct transport regimes at low energies (i.e. $\varepsilon\gt2\lambda$ and $\varepsilon\lt2\lambda$) 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].

Zoom In Zoom Out Reset image size

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 ($\varepsilon\tau\gg1\gg\lambda\tau$) [37]. The condition $\varepsilon\gg\lambda$ 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

$$\begin{equation} \hat{\rho}_{\alpha\beta}(\mathbf{x};t) = \Psi^{\dagger}(\mathbf{\mathbf{x}},t)\,\mathcal{\gamma}_{\alpha\beta}\,\Psi(\mathbf{x},t)\,, \end{equation} \tag{ 8 }$$

where $\mathcal{\mathcal{\gamma}}_{\alpha\beta} = \sigma_{\alpha}\otimes s_{\beta}$ (with $\alpha,\beta = 0,x,y,z$) span the vector space of Hermitian $4\times4$ matrices. The density matrix characterizing a given nonequilibrium state can be expanded as a linear combination of the Clifford algebra basis elements, such that

$$\begin{equation} \rho(\mathbf{x},t) = \frac{1}{d}\sum_{\alpha,\beta = 0,x,y,z}\gamma_{\alpha\beta}\,\langle\hat{\rho}_{\alpha\beta}(\mathbf{x},t)\rangle, \end{equation} \tag{ 9 }$$

where $\langle{\ldots}\rangle$ denotes quantum and disorder averages and $d = \textrm{dim}\,H\equiv4$ is a normalization factor. The expectation value of a generic local observable, $\mathcal{O} = \sum_{\alpha\beta}\mathcal{O}_{\alpha\beta}\gamma_{\alpha\beta}$, is obtained according to

$$\begin{equation} \langle\mathcal{O}(\mathbf{x},t)\rangle: = \textrm{tr}\left[\mathcal{O}\,\rho(\mathbf{x},t)\right] = \sum_{\alpha,\beta}\mathcal{O}_{\alpha\beta}\,\rho_{\alpha\beta}(\mathbf{x},t)\,, \end{equation} \tag{ 10 }$$

where $\textrm{tr}$ indicates the trace over internal degrees of freedom and $\rho_{\alpha\beta}(\mathbf{x},t)\equiv\langle\hat{\rho}_{\alpha\beta}(\mathbf{x},t)\rangle$.

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. $\langle\delta\mathcal{O}(\mathbf{x},t)\rangle: = \textrm{tr}\left[\mathcal{O}\,\delta\rho(\mathbf{x},t)\right]$ with $\delta\rho_{\alpha\beta}(\mathbf{x},t): = \rho_{\alpha\beta}(\mathbf{x},t)-\rho_{\alpha\beta}^{0}$. Here, $\rho_{\alpha\beta}^{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, $\mathcal{J}_{i}^{a}$. The corresponding expectation values away from equilibrium are defined as

$$\begin{equation} \{\delta N,\delta S^{a},\delta J_{i},\delta\mathcal{J}_{i}^{a}\} = \{-e\delta\rho_{00},\,\frac{\hbar}{2}\delta\rho_{0a},\,-ev\delta\rho_{i0},\,\frac{\hbar v}{2}\delta\rho_{ia}\}\,, \end{equation} \tag{ 11 }$$

where we have temporarily reinstated $\hbar$ and e to distinguish between charge and spin currents.

The interaction Hamiltonian is $V(t) = \int d^{2}\mathbf{x}\,\Psi^{\dagger}(\mathbf{x})\,\left[\mathcal{H}_{\textrm{int}}(\mathbf{x},t)\right]\Psi(\mathbf{x})$ with

$$\begin{equation} \mathcal{H}_{\textrm{int}}(\mathbf{x},t) = v\sum_{\alpha,\beta}\gamma_{\alpha\beta}\,\mathcal{A}_{\alpha\beta}^{\textrm{ext}}(\mathbf{x},t)\,, \end{equation} \tag{ 12 }$$

where $\mathcal{A}_{\alpha\beta}^{\textrm{ext}}(\mathbf{x},t)$ ($\alpha,\beta = 0,x,y,z$) spans the Clifford algebra and thus describes any type of charge–spin perturbation applied to the system. In section 2.1, we shall show that the linear response density matrix $\delta\rho(\mathbf{x},t)$, when properly coarse-grained over typical length and time scales (i.e. $|\mathbf{x}|\gg l\equiv v\tau$ and $t\gg\tau$), is governed by an enlarged $16\times16$ 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 section 2.2, the linear response machinery will be applied to derive the full set of drift–diffusion transport equations for the variables $\delta N$, $\delta S^{a}$, $\delta J_{i}$ and $\delta\mathcal{J}_{i}^{a}$, and thus establish a rigorous microscopic picture for the coarse-grained dynamics of the problem. In the following, we define $N\equiv\delta N$, $S^{a}\equiv\delta S^{a}$, $J_{i}\equiv\delta J_{i}$ and $\mathcal{J}_{i}^{a}\equiv\delta\mathcal{J}_{i}^{a}$ for ease of notation.

We start by setting up the formalism needed to derive a quantum kinetic equation for BR-coupled 2D Dirac fermions. From standard linear response theory, the zero temperature density matrix is given by

$$\begin{equation} \delta\rho_{\alpha\beta}(\mathbf{x},t) = \sum_{\bar{\alpha},\bar{\beta} = 0,x,y,z}\int d{\textbf{x}}^{^{\prime}}\int dt^{^{\prime}}\:\mathcal{R}_{\alpha\beta,\bar{\alpha}\bar{\beta}}({\textbf{x}}-{\textbf{x}}^{^{\prime}},t-t^{^{\prime}})\,v\mathcal{\mathcal{A}}_{\bar{\alpha}\bar{\beta}}^{\textrm{ext}}({\textbf{x}}^{^{\prime}},t^{^{\prime}})\,, \end{equation} \tag{ 13 }$$

where

$$\begin{equation} \mathcal{R}_{\alpha\beta,\bar{\alpha}\bar{\beta}}({\textbf{x}}-{\textbf{x}}^{^{\prime}},t-t^{^{\prime}}) = -\imath\theta(t-t^{^{\prime}})\langle\overline{[\hat{\rho}_{\alpha\beta}(\mathbf{x},t),\hat{\rho}_{\bar{\alpha}\bar{\beta}}(\mathbf{x}^{^{\prime}},t^{^{\prime}})]}\rangle\, \end{equation} \tag{ 14 }$$

is the retarded response function. Equation (14) is best evaluated in momentum–frequency space using the Green's-function method [66]. A summation of noncrossing two-particle diagrams, as depicted in figures 2(c) and (d), leads to

$$\begin{equation} \mathcal{R}_{\alpha\beta,\gamma\delta}(\mathbf{q},\omega)\simeq\frac{\omega}{2\pi i}\,\sum_{\mathbf{k}}\,\textrm{tr}\left\{\gamma_{\alpha\beta}\,\mathcal{G}_{\mathbf{k+q}}^{R}(\varepsilon+\omega)\,\tilde{\gamma}_{\gamma\delta}(\mathbf{q},\omega)\,\mathcal{G}_{\mathbf{k}}^{A}(\varepsilon)\right\} \,, \end{equation} \tag{ 15 }$$

where the renormalized vertex operator $\tilde{\gamma}_{\alpha\beta}(\mathbf{q},\omega)$ satisfies the Bethe–Salpeter equation

$$\begin{equation} \tilde{\gamma}_{\alpha\beta}(\mathbf{q},\omega) = \gamma_{\alpha\beta}+n_{i}u_{0}^{2}\,\sum_{{\textbf{k}}}\,\mathcal{G}_{\mathbf{k+\mathbf{q}}}^{R}(\varepsilon+\omega)\,\tilde{\gamma}_{\alpha\beta}(\mathbf{q},\omega)\,\mathcal{G}_{\mathbf{k}}^{A}(\varepsilon)\,. \end{equation} \tag{ 16 }$$

Next, we project both sides of equation (16) onto the Dirac matrices and define

$$\begin{align} \tilde{\boldsymbol{\gamma}}_{\alpha\beta} & (\mathbf{q},\omega) = (\tilde{\gamma}_{\alpha\beta00}(\mathbf{q},\omega),...,\tilde{\gamma}_{\alpha\beta zz}(\mathbf{q},\omega))^{T}\,, \end{align} \tag{ 17 }$$

with $\tilde{\gamma}_{\alpha\beta\varrho\varsigma}(\mathbf{q},\omega) = (1/d)\,\textrm{tr}[\gamma_{\varrho\varsigma}\,\tilde{\gamma}_{\alpha\beta}(\mathbf{q},\omega)]$, to obtain the renormalized vertex in a suitable closed form

$$\begin{equation} \tilde{\mathbf{\boldsymbol{\gamma}}}_{\alpha\beta}(\mathbf{q},\omega) = (I-[M(\mathbf{q},\omega)]^{T})^{-1}\boldsymbol{\gamma}_{\alpha\beta}\,, \end{equation} \tag{ 18 }$$

where $\boldsymbol{\gamma}_{\alpha\beta} = (0,{\ldots},1,{\ldots},0)^{T}$ is an auxiliary vector with nonzero component $(\boldsymbol{\gamma}_{\alpha\beta})_{\alpha\beta} = 1$, I is the identity matrix, T denotes the transpose operation and M is a square matrix of 'bubbles' with elements

$$\begin{equation} M_{\mu\nu,\varrho\varsigma}(\mathbf{q},\omega) = \frac{n_{i}u_{0}^{2}}{d}\sum_{{\textbf{k}}}\textrm{tr}\left[\mathcal{G}_{\mathbf{k+\mathbf{q}}}^{R}(\varepsilon+\omega)\gamma_{\mu\nu}\mathcal{G}_{\mathbf{k}}^{A}(\varepsilon)\gamma_{\varrho\varsigma}\right]\,. \end{equation} \tag{ 19 }$$

We may now recast the Fourier space response function into a more compact form $\mathcal{R}_{\alpha\beta,\bar{\alpha}\bar{\beta}}(\mathbf{q},\omega) =$ $-\imath\omega\nu_{0}\mathcal{D}_{\alpha\beta,\bar{\alpha}\bar{\beta}}(\mathbf{q},\omega)$, where $\nu_{0} = \varepsilon/\pi v^{2}$ is the density of states per spin and

$$\begin{equation} \mathcal{D}_{\alpha\beta,\bar{\alpha}\bar{\beta}}(\mathbf{q},\omega) = \tau\left(\tilde{\gamma}_{\bar{\alpha}\bar{\beta}\alpha\beta}(\mathbf{q},\omega)-\delta_{\alpha\bar{\alpha}}\delta_{\beta\bar{\beta}}\right)\, \end{equation} \tag{ 20 }$$

is the so-called diffuson (here, $\delta_{\alpha\beta}$ is the Kronecker delta symbol). Its inverse, the diffuson Hamiltonian $\mathcal{\mathcal{H}_{\mathcal{D}}}\equiv\mathcal{D}^{-1}$,

$$\begin{align} \mathcal{H}_{\mathcal{D}}(\mathbf{q},\omega) & = \frac{1}{\tau}\left(\left(I-[M(\mathbf{q},\omega)]^{T}\right)^{-1}-I\right)^{-1}, \end{align} \tag{ 21 }$$

provides the kernel of the linear response quantum kinetic equation

$$\begin{equation} \mathcal{H}_{\mathcal{D}}(\mathbf{q},\omega)\cdot\delta\vec{\rho}(\mathbf{q},\omega) = -\imath v\omega\nu_{0}\vec{\mathcal{A}}_{\textrm{ext}}(\mathbf{q},\omega)\,, \end{equation} \tag{ 22 }$$

where $\delta\vec{\rho}(\mathbf{q},\omega)\equiv(\delta\rho_{00}(\mathbf{q},\omega),{\ldots},,\delta\rho_{zz}(\mathbf{q},\omega))^{T}$ and $\vec{\mathcal{A}}_{\textrm{ext}}(\mathbf{q},\omega)\equiv(\mathcal{A}_{00}^{\textrm{ext}}(\mathbf{q},\omega),{\ldots},A_{zz}^{\textrm{ext}}(\mathbf{q},\omega))^{T}$ are the Fourier-space components of the one-particle density matrix and generalized external vector potential, respectively.

Let us briefly discuss the general structure of $\mathcal{\mathcal{H}_{\mathcal{D}}}(\mathbf{q},\omega)$ in the long wavelength limit of interest to us. Zero entries in the $16\times16$ bubble matrix $M(0,\omega)$ can be readily identified by applying the following $C_{6v}$ point-group operations [72]: (i) C₂ rotation exchanging sublattices and (ii) mirror-reflection R_x leaving sublattices invariant. For example, the bubbles $M_{00,0a}(0,\omega)$, encoding charge-density–spin-density-type responses ($a = x,y,z$), vanish identically because the associated vertices ('00' and '0a') transform differently under at least one unitary symmetry. Overall, there are 64 non-zero bubbles in the long-wavelength $\mathbf{q}\rightarrow0$ limit. The symmetry classification is summarized in table 1.

Table 1. Classification of Dirac matrices and corresponding observables under C₂-rotation, mirror-reflection R_x and time-reversal operation. Sublattice-staggered densities are indicated with subscript 's'.

	γ00	γ01	γ02	γ03	γ10	γ11	γ12	γ13	γ20	γ21	γ22	γ23	γ30	γ31	γ32	γ33
$\mathcal{O}$	N	Sx	Sy	Sz	Jx	$\mathcal{J}_{x}^{x}$	$\mathcal{J}_{x}^{y}$	$\mathcal{J}_{x}^{z}$	$J_{y}^ {}$	$\mathcal{J}_{y}^{x}$	$\mathcal{J}_{y}^{y}$	$\mathcal{J}_{y}^{z}$	Ns	$S_{s}^{x}$	$S_{s}^{y}$	$S_{s}^{z}$
C2	+1	−1	−1	+1	−1	+1	+1	−1	−1	+1	+1	−1	+1	−1	−1	+1
Rx	+1	−1	+1	−1	+1	−1	+1	−1	−1	+1	−1	+1	−1	+1	−1	+1
$\mathcal{T}$	+1	−1	−1	−1	−1	+1	+1	+1	−1	+1	+1	+1	−1	+1	+1	+1

We now turn to the spin–charge eigenmodes sustained by the system. The first step is to compute the gradient expansion of equation (21). Working in the diffusive regime ($\epsilon\tau\gg1\gg\lambda\tau$) greatly simplifies matters due to many bubbles being parametrically small. To leading order in $v|\mathbf{q}|$ and ωτ, we find

$$\begin{equation} \mathcal{H}_{\mathcal{D}}(\mathbf{q},\omega)\simeq\imath\omega Q-\imath\mathbf{P}\cdot\mathbf{q}+L\,, \end{equation} \tag{ 23 }$$

where Q, P_i and L are $16 \times 16$ matrices given by $Q = -\imath(\nabla_{\omega}\mathcal{H}_{\mathcal{D}}(0,\omega))_{\mathbf{\omega} = 0}$, $\mathbf{P} = \imath(\partial_{\mathbf{q}}\mathcal{H}_{\mathcal{D}}(\mathbf{q},0))_{\mathbf{q} = 0}$ and $L = \mathcal{H}_{\mathcal{D}}(0,0)$. The calculation of these matrices is rather cumbersome, yielding unwieldy expressions for the matrices $Q = Q(\varepsilon,\lambda,\tau)$, $P_{i} = P_{i}(\varepsilon,\lambda,\tau)$ and $L = L(\varepsilon,\lambda,\tau)$. We therefore provide the explicit form of the equation (23) in the appendix, from which Q, P_i and L can be inferred.

The diffuson Hamiltonian in equation (23) is linear in the wavevector q due to the Dirac nature of the low-lying excitations in graphene heterostructures. The diffuson spectrum at low energies is shown in figure 3. The quadratic dispersion of the gapless mode follows from the diffusive pole structure of the density–density response, i.e. $[\mathcal{R}_{00,00}(\mathbf{q},0)]^{-1}\sim(v\tau|\mathbf{q}|)^{2}$. For $\mathbf{q}\neq0$, this mode is an admixture of charge N, spin current $\mathcal{J}_{i}^{a}$ 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 $\mathbf{q} = 0$ due to the rotational ($C_{v\infty}$) symmetry of the single-particle Hamiltonian (equation (1)). Their gaps at $\mathbf{q} = 0$ are given by $\Delta_{s}\equiv\Delta_{\parallel}\simeq(\lambda/\eta)^{2}$ for the $S^{x,y}$- and $\Delta_{\perp}\simeq2\Delta_{\parallel}$ 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].

Zoom In Zoom Out Reset image size

The low-lying modes shown in figure 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 their distinctly different diffuson Hamiltonians (i.e. $4 \times 4$ Schrödinger-like for 2D electron gases and 1$6 \times 16$ Dirac-like for graphene). Defining $\xi = 2v|\mathbf{q}|\tau$, the eigenvalues of $\mathcal{\mathcal{H}_{\mathcal{D}}}(\mathbf{q},0)$ in the limit $v|\mathbf{q}|\ll1$ are given by $\epsilon_{D}^{0} = \xi^{2}$, $\epsilon_{D}^{1} = \xi^{2}+\Delta_{s}$ and $\epsilon_{D}^{\pm} = \xi^{2}+\Delta_{s}\left(\frac{3}{2}\pm\frac{1}{2}\sqrt{1+16\xi^{2}/\Delta_{s}}\right)$. The Rashba diffuson eigenvalues derived by Wenk et al [51] are recovered by letting $\lambda/v\rightarrow2m_\mathrm{e}\alpha_{2}$, where $m_\mathrm{e}$ is the effective electron mass and α₂ the Rashba parameter. The existence of such a mapping reflects the same basic spin-relaxation (Dyakonov–Perel) mechanism at work. Indeed, our findings put on a firm ground previous heuristic arguments [36, 37] for the equivalence of two models in the weak SOC regime. The remaining $16-4 = 12$ modes in the Dirac–Rashba model are characterized by very large gaps ($\Delta_{s}^{+}\gtrsim1/2\gg\Delta_{s}$), and as such play no role in the diffusive regime.

A comment is in order regarding the validity of our assumptions in the light of recent findings. In the example of figure 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 [74], as well as early magnetotransport measurements [20–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]. We note that in any case, 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 figure 1 is justified) and the system is sufficiently disordered so that $\lambda\tau\ll1$.

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 ($-\imath\omega\rightarrow\partial_{t}$ and $\imath q_{i}\rightarrow\nabla_{i}$) of equation (23) as

$$\begin{equation} \mathcal{H}_{\mathcal{D}}\cdot\delta \vec{\rho}(\mathbf{x},t) = \left(-Q\frac{\partial}{\partial t}-\sum_{i = x,y}P_{i}\nabla_{i}+L\right)\cdot\delta\vec{\rho}(\mathbf{x},t) = v\nu_{0}\partial_{t}\vec{\mathcal{A}}_{\textrm{ext}}(\mathbf{x},t)\,. \end{equation} \tag{ 24 }$$

Transport equations for the coarse-grained variables $N(\mathbf{x},t) = -\delta\rho_{00}(\mathbf{x},t)$, $S^{a}(\mathbf{x},t) = \frac{1}{2}\delta\rho_{0a}(\mathbf{x},t)$, $J_{i}(\mathbf{x},t) = -v\delta\rho_{i0}(\mathbf{x},t)$ and $\mathcal{J}_{i}^{a}(\mathbf{x},t) = \frac{v}{2}\delta\rho_{ia}(\mathbf{x},t)$ can now be derived by replacing the $16 \times 16$ matrices $Q = Q(\varepsilon,\lambda,\tau)$, $P_{i} = P_{i}(\varepsilon,\lambda,\tau)$ and $L = L(\varepsilon,\lambda,\tau)$ with their explicit forms (see appendix). To leading order in $1/\varepsilon\tau$ and λτ, we obtain

$$\begin{equation} \partial_{t}N+\nabla\cdot J = 0\,, \end{equation} \tag{ 25 }$$

$$\begin{equation} \partial_{t}S^{a}+\nabla\cdot\mathcal{J}^{a} = \frac{2\lambda}{v}\,\Lambda_{c}^{ab}\mathcal{J}_{b}^{c}+\frac{\nu_{0}}{2}\mathcal{B}^{a}\,, \end{equation} \tag{ 26 }$$

$$\begin{equation} \partial_{t}J_{i}+\frac{v^{2}}{2}\nabla_{i}N = -\frac{J_{i}}{2\tau}-\frac{4\lambda^{2}}{\epsilon}\epsilon_{ij}\mathcal{J}_{j}^{z}+\frac{\nu_{0}}{2}v^{2}E_{i}\,, \end{equation} \tag{ 27 }$$

$$\begin{equation} \partial_{t}\mathcal{J}_{i}^{a}+\frac{v^{2}}{2}\nabla_{i}S^{a}+v\lambda\tau\,\Gamma_{ic}^{ab}\nabla_{b}\mathcal{J}_{i}^{c} = -\frac{\mathcal{J}_{i}^{a}}{2\tau}+v\lambda\varOmega_{ic}^{a}S^{c}+\lambda^{2}\tau\,\varUpsilon_{ic}^{ab}\mathcal{J}_{b}^{c}-\delta_{az}\epsilon_{ij}\frac{\lambda^{2}}{\varepsilon}J_{j}\,, \end{equation} \tag{ 28 }$$

where $E_{i} = -\partial_{t}\mathcal{A}_{i0}^{\textrm{ext}}$ is the externally applied electric field and $\mathcal{B}^{a} = (1/v)\thinspace\partial_{t}\mathcal{A}_{0a}^{\textrm{ext}}(\mathbf{x},t)$ is the external Zeeman field ('spin injection field') that induces a nonequilibrium spin density [75]. The coefficients $\Lambda_{c}^{ab},\Gamma_{ic}^{ab},\varOmega_{ic}^{a}$ and $\varUpsilon_{ic}^{ab}$ 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 $\mathcal{E}_{i}^{a} = -\partial_{t}\mathcal{A}_{ia}$) can be easily incorporated via a suitable parameterization of the generalized external vector potential entering equation (24).

Table 2. Coefficients in the coupled spin–charge drift–diffusion equations. Only nonzero components are listed.

Λ	Spin density precession	$\Lambda_{x}^{zx} = \Lambda_{y}^{zy} = -1\,$,  $\Lambda_{z}^{xx} = \Lambda_{z}^{yy} = 1$
Γ	Spin current transfer	$\Gamma_{xz}^{xx} = \Gamma_{yz}^{yy} = 1$,  $\Gamma_{xx}^{zx} = \Gamma_{xy}^{zy} = \Gamma_{yx}^{zx} = \Gamma_{yy}^{zy} = -1,$ $\Gamma_{xz}^{yy} = \Gamma_{yz}^{xx} = 1$
$\varOmega$	Spin current precession	$\varOmega_{xz}^{x} = \varOmega_{yz}^{y} = 1$,  $\varOmega_{xx}^{z} = \varOmega_{yy}^{z} = -1$
$\varUpsilon$	Spin current swapping	$\varUpsilon_{xy}^{xy} = \varUpsilon_{yx}^{yx} = 1$,   $\varUpsilon_{xx}^{yy} = \varUpsilon_{yy}^{xx} = -1$

The drift–diffusion equations (25)–(28) are the main results of this work. Equations (25) and (26) are generalized continuity relations and equations (27) and (28) express the time evolution of the currents as a sum of drift, diffusion, spin precession and spin–charge conversion processes. Standard constitutive relations

$$\begin{equation} J_{i} = -D\,\nabla_{i}N-\frac{8\tau\lambda^{2}}{\varepsilon}\epsilon_{ij}\mathcal{J}_{j}^{z}+\nu_{0}\tau v^{2}E_{i}\,, \end{equation} \tag{ 29 }$$

$$\begin{equation} \mathcal{J}_{i}^{a} = -D\nabla_{i}S^{a}-2v\lambda\tau^{2}\,\Gamma_{ic}^{ab}\nabla_{b}\mathcal{J}_{i}^{c}+2v\lambda\tau\varOmega_{ic}^{a}S^{c}+2\lambda^{2}\tau^{2}\,\varUpsilon_{ic}^{ab}\mathcal{J}_{b}^{c}-\delta_{a,z}\epsilon_{ij}\frac{2\tau\lambda^{2}}{\varepsilon}J_{j}\,, \end{equation} \tag{ 30 }$$

where $D = v^{2}\tau$ is the diffusion constant [76], hold to good accuracy insofar as $\omega\ll\lambda$. Sublattice-staggered charge $N_\mathrm{s}$ and spin densities $S_{s}^{a}$ (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–36], which is beyond the scope of the this article.

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–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 $\theta_{\textrm{sH}}^{\textrm{int}} = 2\tau\lambda^{2}/\varepsilon$, which appears explicitly in equations (29)–(30), diverges in the clean limit ($\tau\rightarrow\infty$) 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,

$$\begin{equation} \theta_{\textrm{sH}}(\omega): = 2\epsilon_{ij}\left.\frac{\mathcal{J}_{i}^{z}}{J_{j}}\right|_{\mathbf{q} = 0,\mathcal{B}^{a} = 0} = \;\theta_{\textrm{sH}}^{\textrm{int}}+\theta_{\textrm{sH}}^{\textrm{dis}}(\omega)\,, \end{equation} \tag{ 31 }$$

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 $\theta_{\textrm{sH}}^{\textrm{dis}}(0) = -\theta_{\textrm{sH}}^{\textrm{int}}$ 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 section 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 (figure 1(c)). While the equilibrium spin polarization averaged over the Fermi surface in a nonmagnetic system must vanish identically (equation (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 equation (30) by replacing the pure spin current by its steady-state value in the minimal model, i.e. $\mathcal{J}_{i}^{a} = 0$. The charge-to-spin conversion efficiency is obtained as:

$$\begin{equation} \kappa_{ij}: = 2v\left.\frac{S^{i}}{J_{j}}\right|_{\mathbf{q} = 0,\mathcal{B}^{a} = 0} = \epsilon_{ij}\frac{1}{\lambda\tau}\theta_{\textrm{sH}}^{\textrm{int}} = \epsilon_{ij}\frac{2\lambda}{\varepsilon}\,. \end{equation} \tag{ 32 }$$

This relation (first derived in [26]) discloses an optimal spin-charge conversion efficiency at the spin-gap edge, i.e. $\kappa_{xy}(\varepsilon = 2\lambda) = 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–30].

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 (equations (26)–(28)) subject to time-dependent electric and spin-injection fields, with Fourier transforms $E_{i}(\omega)$ ($i = x,y$) and $\mathcal{B}_{a}(\omega)$ ($a = x,y,z$), respectively. In the Dyakonov–Perel regime with $\omega\tau\ll\lambda\tau\ll1\ll\epsilon\tau$, the macroscopic observables of interest are found, after tedious but straightforward calculations, to be

$$\begin{align} J_{i}(\omega) & = \frac{g_{v}\varepsilon\tau}{\pi}\,\left[\frac{E_{i}(\omega)}{1-2\imath\omega\tau}+\epsilon_{ij}\frac{2\lambda}{v\varepsilon}\frac{\mathcal{B}_{j}(\omega)}{\left(1-\imath\omega\tau_{\parallel}\right)\left(1-2\imath\omega\tau\right)}\right]\,, \end{align} \tag{ 33 }$$

$$\begin{align} \nonumber \\ S_{a}(\omega) & = \frac{g_{v}\varepsilon\tau_{s,a}}{2\pi v^{2}}\,\left[\epsilon_{jaz}\,\frac{8v\lambda^{3}\tau^{2}}{\varepsilon}\frac{E_{j}(\omega)}{\left(1-2\imath\omega\tau\right)\left(1-\imath\omega\tau_{s,a}\right)}+\frac{\mathcal{B}_{a}(\omega)}{1-\imath\omega\tau_{s,a}}\right]\,, \end{align} \tag{ 34 }$$

$$\begin{align} \nonumber \\ \mathcal{J}_{i}^{z}(\omega) & = \frac{g_{v}\varepsilon\tau_{ss}^{z}}{2\pi v}\,\left[\imath v\epsilon_{ij}\frac{\omega\tau}{\varepsilon\tau_{ss}^{z}}\frac{E_{j}(\omega)}{\left(1-2\imath\omega\tau\right)\left(1-\imath\omega\tau_{\parallel}\right)}-\frac{\mathcal{B}_{i}(\omega)}{1-\imath\omega\tau_{\parallel}}\right]\,, \end{align} \tag{ 35 }$$

$$\begin{align} \nonumber \\ \mathcal{J}_{i}^{i}(\omega) & = \frac{g_{v}\varepsilon\tau_{ss}^{i}}{2\pi v}\,\frac{\mathcal{B}_{z}(\omega)}{1-\imath\omega\tau_{\perp}}\,, \end{align} \tag{ 36 }$$

where the factor of $g_{v} = 2$ accounts for valley degeneracy. Furthermore, _zij denotes the rank-3 Levi-Civita symbol, $\tau_{ss}^{i = x,y} = \tau_{ss}^{z}/2$ with $\tau_{ss}^{z} = 1/2\lambda$, and $\tau_{s,a}$ is the Dyakonov–Perel relaxation time introduced earlier [here, $\tau_{s,(x,y)}\equiv\tau_{\parallel}$ and $\tau_{s,z}\equiv\tau_{\perp} = \tau_{\parallel}/2$, with $\tau_{\parallel} = (4\lambda^{2}\tau)^{-1}$].

Equations (33)–(36) reveal several interesting features of the Dirac–Rashba model. First, a finite SHE is established in the presence of a time-dependent electric field (equation (35)). Second, a spin-injection field gives rise to a charge current (i.e. direct spin galvanic effect) (equation (33)) and a nonequilibrium spin polarization (equation (34)). Moreover, a pure spin current is induced by an applied spin-injection field (equations (35)–(36)). The induced spin-current density is polarized transversely to the spin-injection field. With the exception of the SHE, all such effects are present in the DC limit.

The linear response functions to external fields can be readily obtained from equations (33)–(36). For example, the optical conductivity [$\sigma_{ij}(\omega) = \delta J_{i}(\omega)/\delta E_{j}(\omega)$], spin-galvanic susceptibility [$\chi_{ij}(\omega) =$ $\delta S_{i}(\omega)/\delta E_{j}(\omega)$] and spin Hall conductivity [$\sigma_{yx}^{z}(\omega) = \delta\mathcal{J}_{y}^{z}(\omega)/\delta E_{x}(\omega)$] read as

$$\begin{equation} \sigma_{ii}(\omega) = 2g_{v}\,\frac{\varepsilon\tau}{1-2\imath\omega\tau}\,,\qquad(\mathrm{units\,of\,}e^{2}/h) \end{equation} \tag{ 37 }$$

$$\begin{equation} \chi_{ij}(\omega)\simeq-\frac{2g_{v}\epsilon_{ij}}{v}\,\frac{\lambda\tau}{1-\imath\omega\tau_{\parallel}}\,,\qquad(\mathrm{units\,of\,}e/2\pi) \end{equation} \tag{ 38 }$$

$$\begin{equation} \sigma_{ij}^{z}(\omega)\simeq-\imath\epsilon_{ij}g_{v}\,\frac{\omega\tau}{1-\imath\omega\tau_{\parallel}}\,,\qquad(\mathrm{units\,of\,}e/2\pi) \end{equation} \tag{ 39 }$$

where subleading corrections of order ωτ have been neglected for simplicity. These results deserve a few comments. The expression for $\sigma_{xx}(\omega)$ coincides with the familiar Drude model result, which is valid in the semiclassical regime with $\varepsilon\tau\gg1$ [82]. More interestingly, the optical spin-galvanic susceptibility (equation (38)) generalizes the findings of [26] to an external electric field with nonzero frequency. Here, $\omega\tau_{\parallel}$ 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. $\theta_{p} = \lambda\tau$).