\section*{Appendix A: Electronic structure and pseudospin-spin texture}

The low-energy Hamiltonian of $C_{3v}$-invariant vdW monolayers reads
as
\begin{equation}
H_{\xi}=\Sigma_{\mu}\,\left(vp+\mathcal{A}_{\textrm{SO}}+\mathcal{A}_{\textrm{xc}}+\mathcal{A}_{\textrm{orb}}\right)^{\mu}\,,\label{eq:1}
\end{equation}
where $p^{\mu}\equiv(-\epsilon/v,p_{x},p_{y})$ is the 3-momentum
of the interface, $\Sigma_{\mu}=(\Sigma_{0},\vec{\Sigma})$ and 
\begin{align}
 & \mathcal{A}_{\textrm{SO}}=\lambda\,(s_{y}\hat{x}-s_{x}\hat{y})+\lambda_{\textrm{KM}}s_{z}\hat{z}+\xi\lambda_{\textrm{sv}}s_{z}\hat{t}\,,\label{eq:A_SO}\\
 & \mathcal{A}_{\textrm{ex}}+\mathcal{A}_{\textrm{orb}}=(\boldsymbol{\mathfrak{m}}\cdot\mathbf{s})\hat{t}+\xi\Delta\,\hat{z}\,,\label{eq:A_2}
\end{align}
are non-Abelian gauge fields capturing all symmetry-allowed SOCs {[}Eq.
(\ref{eq:A_SO}){]}, on-site staggered potential and interfacial exchange
coupling {[}Eq. (\ref{eq:A_2}){]}. Here, $\boldsymbol{\mathfrak{m}}=-\Delta_{\textrm{xc}}\,\mathbf{m}\equiv\mathfrak{m}(\sin\phi\,\hat{\mathbf{x}}+\cos\phi\,\hat{\mathbf{z}})$
with $\pi\ge\phi\ge0$ parametrizes the exchange field. For brevity,
in this supplementary information, the analytical expressions are
provided for the strong SOC regime with $\lambda>\mathfrak{\tilde{m}}_{z}\gg\mathfrak{m}_{x}$
and $\lambda_{\textrm{KM}}=0$, where $\mathfrak{\tilde{m}}_{z}\equiv\mathfrak{m}_{z}+\xi\lambda_{\textrm{sv}}$.
Figure \ref{fig:spectrum} shows the low-energy spectrum for two representative
systems: (a) TMD/graphene/FM and (b) TMD/FM. In panel (a) only low-energy
states (within the TMD gap) are shown. In-plane magnetization ($\mathfrak{m}_{x}$)
breaks the $C_{v\infty}$ symmetry of the continuum model, rending
the Fermi surface anisotropic.

The electronic structure comprises four distinct spectral regions:

\bigskip{}

\begin{itemize}
\item Regime Ia: Low energy regime where the Fermi level crosses an electron/hole
pocket for 
\begin{equation}
\epsilon_{\text{Ia}}<|\epsilon|<\epsilon_{\text{\text{Ib}}};
\end{equation}
\item Regime Ib: Very narrow energy range where the Fermi level crosses
two different Fermi rings both belonging to the spin majority band.
This happens for
\begin{equation}
\epsilon_{\text{\text{Ib}}}<|\epsilon|<\epsilon_{\text{\text{Ic}}};
\end{equation}
\item Regime Ic: Intermediate regime where the Fermi level crosses only
the spin majority band, hinting at stronger spin density responses
for 
\begin{equation}
\epsilon_{\text{Ic}}<|\epsilon|<\epsilon_{\text{II}};
\end{equation}
\item Regime II: Typical high-electronic density regime in the experiments.
Here, we have $\epsilon>\epsilon_{\text{II}}$ and the Fermi level
crosses two Fermi rings with opposite spin textures;
\end{itemize}
\begin{figure}[b]
\includegraphics[width=0.9\columnwidth]{fig_01_supp.pdf}

\caption{\label{fig:spectrum}Electronic structure of a TMD|graphene|FM (a)
and TMD|FM (b) heterostructure plotted along $\mathbf{p}=(0,p_{y})^{\textrm{t}}$
(full line) and $\mathbf{p}=(p_{x},0)^{\textrm{t}}$ (dashed lines)
directions. Shaded areas highlight the spectral region II. Parameters:
for (a-b) $\lambda=40$ meV, $\mathfrak{m}_{x}=20$ meV, $\mathfrak{m}_{z}=30$
meV and for (b) $\Delta=0.4$ eV, $\lambda_{\textrm{sv}}(\epsilon<0)=150$
meV and $\lambda_{\textrm{sv}}(\epsilon>0)=5$ meV.}
\end{figure}
The expressions for the different limits are given in Table \ref{tab:limits}.
In regime II, the spin operators have the following equilibrium average
values at the Fermi energy, in the asymptotic limit $\epsilon\gg\mathfrak{m}_{z}\equiv\tilde{\mathfrak{m}}_{z}\gg\mathfrak{m}_{x}$
(here, $\theta$ is the wavevector angle with respect to $\hat{x}$
axis):
\begin{align}
\left\langle s_{x}\right\rangle  & \cong\frac{\lambda}{\sqrt{\lambda^{2}+\mathfrak{m}_{z}^{2}}}\left(1-\frac{\left(\Delta\mathfrak{m}_{z}+\lambda^{2}\right){}^{2}}{2\epsilon^{2}\left(\lambda^{2}+\mathfrak{m}_{z}^{2}\right)}\right)\sin\theta\nonumber \\
 & +\frac{\mathfrak{m}_{x}}{2\sqrt{\lambda^{2}+\mathfrak{m}_{z}^{2}}}\left(1+\frac{\mathfrak{m}_{z}^{2}+\lambda^{2}\cos2\theta}{\lambda^{2}+\mathfrak{m}_{z}^{2}}\right),
\end{align}

\begin{align}
\left\langle s_{y}\right\rangle  & \cong-\frac{\lambda}{\sqrt{\lambda^{2}+\mathfrak{m}_{z}^{2}}}\left(1-\frac{\left(\Delta\mathfrak{m}_{z}+\lambda^{2}\right){}^{2}}{2\epsilon^{2}\left(\lambda^{2}+\mathfrak{m}_{z}^{2}\right)}\right)\cos\theta\nonumber \\
 & +\frac{\mathfrak{m}_{x}}{2\sqrt{\lambda^{2}+\mathfrak{m}_{z}^{2}}}\frac{\lambda^{2}\sin2\theta}{\lambda^{2}+\mathfrak{m}_{z}^{2}},
\end{align}

\begin{align}
\left\langle s_{z}\right\rangle  & \cong\frac{\mathfrak{m}_{z}}{\sqrt{\lambda^{2}+\mathfrak{m}_{z}^{2}}}\left(1+\frac{\lambda^{2}}{2\epsilon^{2}}\frac{\Delta^{2}-\mathfrak{m}_{z}^{2}}{\lambda^{2}+\mathfrak{m}_{z}^{2}}-\frac{\mathfrak{m}_{x}\lambda\sin\theta}{\lambda^{2}+\mathfrak{m}_{z}^{2}}\right)\nonumber \\
 & +\frac{\lambda^{4}}{\epsilon^{2}}\frac{\Delta-\mathfrak{m}_{z}}{\left(\lambda^{2}+\mathfrak{m}_{z}^{2}\right){}^{3/2}},
\end{align}
for the spin majority band, the other band has opposite polarity.
The pseudospin texture, on the other hand, is:

\begin{equation}
\left\langle \Sigma_{x}\right\rangle \cong\left(1-\frac{\lambda^{2}+\Delta^{2}}{2\epsilon^{2}}\right)\cos\theta\pm\frac{\lambda\mathfrak{m}_{x}}{2\epsilon}\frac{\sin2\theta}{\sqrt{\lambda^{2}+\mathfrak{m}_{z}^{2}}},
\end{equation}
\begin{equation}
\left\langle \Sigma_{y}\right\rangle \cong\left(1-\frac{\lambda^{2}+\Delta^{2}}{2\epsilon^{2}}\right)\sin\theta\mp\frac{\alpha\mathfrak{m}_{x}}{2\epsilon}\frac{1+\cos2\theta}{\sqrt{\lambda^{2}+\mathfrak{m}_{z}^{2}}},
\end{equation}

\begin{equation}
\left\langle \Sigma_{z}\right\rangle \cong\frac{\Delta}{\epsilon}\mp\frac{\mathfrak{m}_{z}}{\sqrt{\lambda^{2}+\mathfrak{m}_{z}^{2}}}\frac{\Delta\mathfrak{m}_{z}+\lambda^{2}}{\epsilon^{2}}\left(1-\frac{\lambda\mathfrak{m}_{x}\sin\theta}{\lambda^{2}+\mathfrak{m}_{z}^{2}}\right),
\end{equation}
for the spin majority ($+$)/minority ($-$) bands. Figure (\ref{fig:spin_text})
shows the spin and pseudospin profiles along de $\hat{x}$ and $\hat{y}$
directions in momentum space. The orbital mass ($\Delta$) broadens
up the $\mathbf{p}$-space spin texture dramatically, which boosts
the generation of out-of-plane spin polarization in applied current.

\begin{table*}[t]
\begin{tabular}{|c|c|c|c|c|}
\hline 
 & $\epsilon_{\text{II}}$ & $\epsilon_{\text{\text{Ic}}}$ & $\epsilon_{\text{\text{Ib}}}$ & $\epsilon_{\text{Ia}}$\tabularnewline
\hline 
\hline 
$\,$$\,$$\Delta=0$$\,$$\,$ & $\sqrt{\mathfrak{\tilde{m}}_{z}^{2}+4\lambda^{2}}$ & $\mathfrak{\tilde{m}}_{z}+\frac{\mathfrak{m}_{x}^{2}}{2\mathfrak{\tilde{m}}_{z}}$ & $\frac{\mathfrak{\tilde{m}}_{z}\lambda}{\sqrt{\mathfrak{\tilde{m}}_{z}^{2}+\lambda^{2}}}+\frac{\lambda\mathfrak{m}_{x}\mathfrak{\tilde{m}}_{z}\sqrt{2\lambda^{2}+\mathfrak{\tilde{m}}_{z}^{2}}}{\left(\lambda^{2}+\mathfrak{\tilde{m}}_{z}^{2}\right)^{3/2}}$ & $\frac{\tilde{\mathfrak{m}}_{z}\lambda}{\sqrt{\mathfrak{\tilde{m}}_{z}^{2}+\lambda^{2}}}-\frac{\lambda\mathfrak{m}_{x}\mathfrak{\tilde{m}}_{z}\sqrt{2\lambda^{2}+\mathfrak{\tilde{m}}_{z}^{2}}}{\left(\lambda^{2}+\mathfrak{\tilde{m}}_{z}^{2}\right)^{3/2}}$\tabularnewline
\hline 
$\,$$\,$$\Delta\neq0$$\,$$\,$ & $\mathfrak{\tilde{m}}_{z}+\Delta$ & $\sqrt{\left(\Delta-\mathfrak{\tilde{m}}_{z}\right)^{2}+4\lambda^{2}}$ & $\frac{\lambda(\Delta+\tilde{\mathfrak{m}}_{z})}{\sqrt{\lambda^{2}+\mathfrak{\tilde{m}}_{z}^{2}}}+\frac{\lambda\mathfrak{m}_{x}\sqrt{\mathfrak{\tilde{m}}_{z}(\Delta+\mathfrak{\mathfrak{\tilde{m}}}_{z})\left(2\lambda^{2}+\mathfrak{\tilde{m}}_{z}^{2}-\Delta\mathfrak{\mathfrak{\tilde{m}}}_{z}\right)}}{\left(\lambda^{2}+\mathfrak{\mathfrak{\tilde{m}}}_{z}^{2}\right)^{3/2}}$ & $\frac{\lambda(\Delta+\tilde{\mathfrak{m}}_{z})}{\sqrt{\lambda^{2}+\mathfrak{m}_{z}^{2}}}-\frac{\lambda\mathfrak{m}_{x}\sqrt{\tilde{\mathfrak{m}}_{z}(\Delta+\tilde{\mathfrak{m}}_{z})\left(2\lambda^{2}+\tilde{\mathfrak{m}}_{z}^{2}-\Delta\mathfrak{\tilde{m}}_{z}\right)}}{\left(\lambda^{2}+\mathfrak{\tilde{m}}_{z}^{2}\right)^{3/2}}$\tabularnewline
\hline 
\end{tabular}

\caption{\label{tab:limits}Spectral regimes of the $C_{3v}$ model. To ease
the notation, all couplings are taken to be positive.}
\end{table*}
\begin{figure}
\includegraphics[width=0.9\columnwidth]{fig_02_supp.pdf}\caption{\label{fig:spin_text}Emergent pseudospin-spin texture of spin-majority
band for a 2D material|FM interface in the presence (thick lines)
and absence (thin lines) of orbital gap. The pseudospin-spin textures
are plotted along path $\mathbf{k}=(k_{x},0)$ (top panels) and $\mathbf{k}=(0,k_{y})$
(bottom panels). Band structure parameters: $\lambda=20$ meV, $\Delta=400$
meV, $\Delta_{\textrm{xc}}=15$meV (with $\phi=\pi/8$) and $\lambda_{\textrm{sv}}=0$.}
\end{figure}

\section*{Appendix B: Symmetry Analysis of Fermi-surface $K$-response tensor }

The Fermi surface density-current response functions are determined
from the zero temperature limit of the term I in the Kubo-Streda formula,
namely
\begin{equation}
K_{ai}=\frac{1}{2\pi}\int(d\mathbf{p})\textrm{Tr}\left\{ s_{a}\mathcal{G}_{\mathbf{p}}^{+}\tilde{J}_{i}\mathcal{G}_{\mathbf{p}}^{-}\right\} \,,\label{eq:B1}
\end{equation}
where $\tilde{J}_{i}$ is the renormalized current density vertex
and $(d\mathbf{p})\equiv\text{d}\mathbf{p}/\left(2\pi\right)^{2}$.
To determine the parity of $K_{ai}$ with respect to the field reversal
$\mathfrak{m}\rightarrow-\mathfrak{m}$, it suffices to consider the
``empty'' bubble \citep{Titov_2DEGRashba_SOT}. The disorder-averaged
Green's functions satisfy the following symmetry relations
\begin{align}
s_{x}\sigma_{y}\mathcal{G}^{a}(-p_{x},p_{y})\sigma_{y}s_{x} & =\left.\mathcal{G}^{a}(p_{x},p_{y})\right|_{S_{1}}\,,\label{eq:6}\\
s_{y}\sigma_{x}\mathcal{G}^{a}(p_{x},-p_{y})s_{y}\sigma_{x} & =\left.\mathcal{G}^{a}(p_{x},p_{y})\right|_{S_{2}}\,,\\
s_{z}\sigma_{z}\mathcal{G}^{a}(-p_{x},-p_{y})s_{z}\sigma_{z} & =\left.\mathcal{G}^{a}(p_{x},p_{y})\right|_{S_{3}}\,,\label{eq:7}
\end{align}
with $S_{1}\equiv\{\mathfrak{m}_{z}\rightarrow-\mathfrak{m}_{z},\lambda_{\textrm{sv}}\rightarrow-\lambda_{\textrm{sv}},\Delta\rightarrow-\Delta\}$,
$S_{2}\equiv\{\mathfrak{m}_{z}\rightarrow-\mathfrak{m}_{z},\mathfrak{m}_{x}\rightarrow-\mathfrak{m}_{x},\lambda_{\textrm{sv}}\rightarrow-\lambda_{\textrm{sv}},\Delta\rightarrow-\Delta\}$
and $S_{3}\equiv\{\mathfrak{m}_{x}\rightarrow-\mathfrak{m}_{x}\}$.
Using these symmetries, we find after some straightforward algebra
(sum over repeated indices $a=x,z$ is implied)\begin{widetext}

\begin{equation}
\{K_{ai}\}=\left(\begin{array}{cc}
\mathfrak{m}_{z}\kappa_{xx}+\mathfrak{m}_{a}^{2}\vec{\alpha}_{\xi}\cdot\vec{f_{xx}^{a}}+\vec{\alpha}_{\xi}\cdot\vec{g}\qquad & \kappa_{xy}+\mathfrak{m}_{z}\vec{\alpha}_{\xi}\cdot\vec{f}_{xy}+\mathfrak{m}_{a}^{2}h_{xy}^{a}\\
-\kappa_{xy}-\mathfrak{m}_{z}\vec{\alpha}_{\xi}\cdot\vec{f}_{xy}+\mathfrak{m}_{a}^{2}h_{yx}^{a} & \mathfrak{m}_{z}\kappa_{xx}+\vec{\alpha}_{\xi}\cdot\vec{g}+\mathfrak{m}_{a}^{2}\vec{\alpha}_{\xi}\cdot\vec{f_{yy}^{a}}\\
\mathfrak{m}_{x}\kappa_{zx}+\mathfrak{m}_{z}\mathfrak{m}_{x}\vec{\alpha}_{\xi}\cdot\vec{z}_{zx}\qquad & \mathfrak{m}_{z}\mathfrak{m}_{x}\kappa_{zy}+\mathfrak{m}_{x}\vec{\alpha}_{\xi}\cdot\vec{z}_{zy}
\end{array}\right),
\end{equation}
\end{widetext}where $\vec{\alpha}_{\xi}=\xi(\lambda_{\textrm{sv}},\Delta)$
and $\{\kappa_{ia},\mathbf{f}_{ij},\mathbf{f}_{ij}^{a},\mathbf{g},h_{ij}^{a},\mathbf{z}_{zi}\}_{i,j=x,y}$
are even functions of $\mathfrak{m}_{x}$ and $\mathfrak{m}_{z}$.
The terms linear in $\vec{\alpha}_{\xi}$ are activated by the breaking
of sublattice symmetry, vanishing upon the summation over the two
valleys. 

\section*{Appendix C: analytical results for weak scattering regime}

\begin{table*}
\begin{centering}
\begin{tabular}{|c|c|c|}
\hline 
 & $m_{x}=0$ & $\mathcal{O}\left(m_{x}\right)$\tabularnewline
\hline 
\hline 
$\,\,\,\tilde{J}_{x}\,\,\,$ & $\,\,\Sigma_{x}s_{0},\,\Sigma_{0}s_{y},\,\Sigma_{x}s_{z},\,\Sigma_{z}s_{y}\,\,$ & $\,\,\Sigma_{x}s_{x},\,\Sigma_{y}s_{y}\,\,$\tabularnewline
\hline 
$\tilde{J}_{y}$ & $\,\,\Sigma_{y}s_{0},\,\Sigma_{0}s_{x},\,\Sigma_{y}s_{z},\,\Sigma_{z}s_{x}\,\,$ & $\,\,\Sigma_{0}s_{0},\,\Sigma_{z}s_{0},\,\Sigma_{0}s_{z},\,\Sigma_{x}s_{y},\,\Sigma_{y}s_{x},\,\Sigma_{z}s_{z}\,\,$\tabularnewline
\hline 
\end{tabular}
\par\end{centering}
\caption{\label{tab:Matrixstructures}Gaussian matrix structures of the renormalized
current vertex. In-plane magnetic coupling ($m_{x}\protect\neq0$)
generates additional orbital-spin mixings, which are fundamental for
the accurate description of SOTs.}
\end{table*}
For 2D materials with $C_{6v}$ symmetry ($\Delta=\lambda_{\textrm{sv}}=0$),
the self energy at high carrier density reads as 
\begin{equation}
\Sigma^{\pm}\left(\epsilon\right)\simeq\mp\imath\eta\left(\boldsymbol{1}+\frac{\mathfrak{m}_{x}}{\epsilon}s_{x}+\frac{\mathfrak{m}_{z}}{\epsilon}s_{z}\right),
\end{equation}
where $\eta=nu_{0}^{2}\epsilon/(4v^{2})$ is the disorder-induced
quasiparticle broadening. There is already a significant difference
at this stage between magnetized Dirac fermions and 2DEGs since, in
the latter, the self-energy is a scalar. For 2D Dirac fermions, it
is fundamental to keep the full matrix structure of $\Sigma^{\pm}$
in order to obtain physically sensible results that comply with exact
symmetry relations of the four-point vertex function, known as Ward
identities \citep{manuel_covariant_laws}. In the main text, we presented
full nonperturbative results obtained with a numerical inversion of
the Bethe-Salpeter equation (Fig.\,\ref{fig:03}). Analytical expressions
for the \emph{weak scattering regime} can be obtained by evaluating
special subsets of diagrams \citep{MilletariFerreira16a}. Three responses,
$K_{xy}$, $K_{yx}$ and $K_{zy}$, are activated at the ``Gaussian
level'' with the standard correlator $\langle V(\mathbf{x})V(\mathbf{x}^{\prime})\rangle=nu_{0}^{2}\delta(\mathbf{x}-\mathbf{x}^{\prime})$
. To capture the antidamping responses $K_{xx},K_{yy}$ and $K_{zx}$,
one supplements the ladder series with the ``Y'' diagrams generated
by the high-order correlator $\langle V(\mathbf{x})V(\mathbf{x}^{\prime})V(\mathbf{x}^{\prime\prime})\rangle=nu_{0}^{3}\delta(\mathbf{x}-\mathbf{x}^{\prime})\delta(\mathbf{x}^{\prime}-\mathbf{x}^{\prime\prime})$.
We find for the Gaussian Fermi-surface responses:
\begin{equation}
K_{xy}=-K_{yx}=\frac{2}{v\eta}\frac{\lambda^{3}\epsilon^{2}\left(\epsilon^{2}+\mathfrak{m}_{z}^{2}\right)}{\epsilon^{4}\left(\lambda^{2}+\mathfrak{m}_{z}^{2}\right)-\mathfrak{m}_{z}^{4}\left(\epsilon^{2}-3\lambda^{2}\right)},
\end{equation}
\begin{equation}
K_{zy}=-\frac{2}{v\eta}\frac{\lambda\mathfrak{m}_{x}\mathfrak{m}_{z}}{\lambda^{2}+\mathfrak{m}_{z}^{2}}+\mathcal{O}\left(\epsilon^{-3}\right)\,,
\end{equation}
\begin{equation}
\sigma_{xx}=\sigma_{yy}=\frac{\epsilon}{\eta}\left(1-\frac{4\lambda^{2}\mathfrak{m}_{z}^{2}\left(\epsilon^{2}-2\lambda^{2}\right)}{\epsilon^{4}\left(\lambda^{2}+\mathfrak{m}_{z}^{2}\right)-\mathfrak{m}_{z}^{4}\left(\epsilon^{2}-3\lambda^{2}\right)}\right)\,,
\end{equation}
and
\begin{equation}
K_{xx}^{\text{Y}}=K_{yy}^{\text{Y}}=\frac{u_{0}\mathfrak{m}_{z}\lambda^{5}\epsilon^{2}\left(\epsilon^{2}-\mathfrak{m}_{z}^{2}\right)\left(\epsilon^{2}+\mathfrak{m}_{z}^{2}\right)^{2}}{v^{3}\eta\left(\epsilon^{4}\left(\mathfrak{m}_{z}^{2}+\lambda^{2}\right)-\mathfrak{m}_{z}^{4}\left(\epsilon^{2}-3\lambda^{2}\right)\right)^{2}},\label{eq:coll}
\end{equation}
\begin{equation}
K_{zx}^{\text{Y}}=-\frac{u_{0}\lambda^{5}\mathfrak{m}_{x}\mathfrak{m}_{z}^{2}}{2v^{3}\eta\left(\lambda^{2}+\mathfrak{m}_{z}^{2}\right)^{3}}+\mathcal{O}\left(\epsilon^{-3}\right),
\end{equation}
\begin{alignat}{1}
\sigma_{xy}^{\text{Y}} & =-\sigma_{yx}^{\text{Y}}=-\frac{2u_{0}\mathfrak{m}_{z}\lambda^{6}\epsilon\left(\epsilon^{2}+\mathfrak{m}_{z}^{2}\right)^{3}}{v^{2}\eta\left(\epsilon^{4}\left(\lambda^{2}+\mathfrak{m}_{z}^{2}\right)-\mathfrak{m}_{z}^{4}\left(\epsilon^{2}-3\lambda^{2}\right)\right)^{2}}\,,
\end{alignat}
for the skew-scattering Fermi-surface terms.

The linear response function at the Gaussian level is determined by
a \emph{single} component of the renormalized vertex \citep{manuel_MDPI},
which provides a transparent scheme to identify candidate nonzero
responses $K_{ia}$ based on a symmetry analysis (see table \ref{tab:Matrixstructures}).
The current vertex transverse to the in-plane magnetic coupling axis
($\tilde{J}_{y}$) displays a complex structure with a term proportional
to $s_{z}$, showing that $K_{zy}$ is finite already at the Gaussian
level, in agreement with the analysis of Boltzmann transport equations
outlined in the main text.


\begin{figure}[H]
\begin{centering}
\includegraphics[scale=0.4]{fig_03_supp.pdf}
\par\end{centering}
\caption{\label{fig:fig_4}\textit{ }Angular dependence of $\mathbf{m}$-odd
response functions (top) and associated SOT efficiencies (bottom).
System parameters: $\Delta_{\textrm{xc}}=0.1\,\text{eV}$, $\lambda=60\text{meV},$
$\lambda_{\textrm{sv}}=\Delta=0$, $n=10^{11}\text{cm}^{-2}$ and
$\epsilon=0.4$ eV.}
\end{figure}

Defining $\gamma_{ij}\equiv\Sigma_{i}s_{j}$, one obtains $\tilde{\gamma}_{ij}=\sum_{kl}c_{ijkl}\gamma_{kl}$
with the following nonzero coefficients at the Gaussian level
\begin{equation}
c_{1010}=c_{2020}=2-\frac{4\lambda^{2}\mathfrak{m}_{z}^{2}\left(\epsilon^{2}-2\lambda^{2}\right)}{\epsilon^{4}\left(\lambda^{2}+\mathfrak{m}_{z}^{2}\right)-\mathfrak{m}_{z}^{4}\left(\epsilon^{2}-3\lambda^{2}\right)},\label{eq:Long_cond}
\end{equation}
\begin{equation}
c_{1002}=-c_{2001}=-\frac{2\lambda^{3}\epsilon\left(\epsilon^{2}+\mathfrak{m}_{z}^{2}\right)}{\epsilon^{4}\left(\lambda^{2}+\mathfrak{m}_{z}^{2}\right)-\mathfrak{m}_{z}^{4}\left(\epsilon^{2}-3\lambda^{2}\right)},\label{eq:edelstein}
\end{equation}
\begin{equation}
c_{2003}=-\frac{2\lambda\mathfrak{m}_{x}\mathfrak{m}_{z}}{\epsilon\left(\lambda^{2}+\mathfrak{m}_{z}^{2}\right)}+\mathcal{O}\left(\epsilon^{-2}\right).\label{eq:zy}
\end{equation}
The longitudinal conductivity is determined by Eq. (\ref{eq:Long_cond}),
the Edelstein (inverse spin galvanic) effect is encoded in Eq. (\ref{eq:edelstein}),
and Eq. (\ref{eq:zy}) determines the generation of out-of-plane nonequilibrium
spin polarization. Replacing the renormalized vertex in the \textquotedblleft bubble\textquotedblright{}
{[}Eq. (\ref{eq:B1}){]} yields the Gaussian response functions presented
earlier. 

\section*{Appendix D: Full angular dependence of current-induced SOT}

From the knowledge of the spin-charge response tensor $\hat{K}^{J}(\phi)=\hat{K}(\phi)\cdot\hat{\sigma}(\phi)^{-1}$,
one can easily extract the SOT efficiencies or ``torkances'' \begin{widetext}
\begin{align}
t_{\textrm{o}1}(\phi)= & K_{xy}^{J}(\phi)-\tan(\phi)\,K_{zy}^{J}(\phi)\,,\label{eq:t_odd_12-1}\\
t_{\textrm{o}2}(\phi)= & \textrm{cosec}^{2}(\phi)\left[K_{xy}^{J}(\phi)+K_{yx}^{J}(\phi)\right]-2\textrm{cosec}(2\phi)K_{zy}^{J}(\phi)\,,\\
t_{\textrm{e}1}(\phi)= & \sec(\phi)K_{yy}^{J}(\phi)\,,\quad\quad\quad\label{eq:t_even_12-1}\\
t_{\textrm{e}2}(\phi)= & \textrm{cosec}^{2}(\phi)\sec(\phi)\left[K_{xx}^{J}(\phi)-K_{yy}^{J}(\phi)\right]-\textrm{cosec}(\phi)K_{zx}^{J}(\phi)-\sec(\phi)K_{xx}^{J}(\phi)\,,
\end{align}
\end{widetext} entering in the final expression for the current-induced
SOT\begin{widetext}
\begin{equation}
\mathbf{T}=t_{\textrm{o}1}(\phi)\,\mathbf{m}\times(\hat{z}\times\mathbf{J})+t_{\textrm{e}1}(\phi)\,\mathbf{m}\times(\mathbf{m}\times(\hat{z}\times\mathbf{J}))+t_{\textrm{o}2}(\phi)\,\mathbf{m}\times(\mathbf{m}\times\hat{z})(\mathbf{m}\cdot\mathbf{J})+t_{\textrm{e}2}(\phi)\,\mathbf{m}\times\hat{z}\,(\mathbf{m}\cdot\mathbf{J})\,.\label{eq:2}
\end{equation}
\end{widetext}

These expressions match those reported by I. A. Ado \emph{et al} \citep{Titov_2DEGRashba_SOT}
(apart from an overall minus sign in the $\mathbf{m}-$even torques).
The evaluation fo the angular dependence of the torkance funcions
$\{t_{\textrm{o}1}(\phi),t_{\textrm{o}2}(\phi),t_{\textrm{e}1}(\phi),t_{\textrm{e}2}(\phi)\}$
requires a full nonperturbative treatment beyond previous microscopic
formulations\textcolor{black}{{} \citep{Barnas_SOT_Graphene_PRB2015,Titov_2DEGRashba_SOT,TI_SOT,PhysRevLett.108.117201}}.
The angular dependence of the skew scattering-activated response function
and SOT efficiency parameters obtained by a numerically exact procedure
is depicted in Fig. \ref{fig:fig_4} for a magnetized graphene layer
with strong interfacial exchange and SOC effects. 

\bibliographystyle{unsrt}
