# Supplementary information to the paper “Quantum interference tuning of spin-
orbit coupling in twisted van der Waals trilayers”

Csaba G. Péterfalvi  Department of Physics, University of Konstanz, D-78464
Konstanz, Germany  Alessandro David  Peter Grünberg Institute - Quantum
Control (PGI-8), Forschungszentrum Jülich GmbH, Jülich, Germany  Péter Rakyta
Department of Physics of Complex Systems, Eötvös Loránd University, Budapest,
Hungary  Quantum Information National Laboratory, Hungary  Guido Burkard
Department of Physics, University of Konstanz, D-78464 Konstanz, Germany
Andor Kormányos  Department of Physics of Complex Systems, Eötvös Loránd
University, Budapest, Hungary

##  I Symmetry properties of the phase $\vartheta_{R}$

In order to be self-contained and to set the stage, in Sec. I.1 and I.2 we
carefully re-derive previous results of Ref. david_induced_2019  for the
Hamiltonian $H_{R}^{gr}$, which describes the proximity induced Rashba SOC in
graphene. We then present two short, but important discussions. First, in Sec.
I.3 we argue that the sign of the SOC coupling $\lambda_{R}$ depends on the
layer stacking order and this is important in trilayer stacks. Next, in Sec.
I.4 we discuss the effects of the interlayer twist on $H_{R}^{gr}$. Our most
important results in this section follow in Sec. I.5. We prove that the
results for the quantum phase $\vartheta_{R}$ simplify for interlayer twist
angles $\theta=0,\pi/6,\pi/3$, where the bilayer stack has $C_{3v}$ symmetry.
Finally, in Sec. I.6 we give a further discussion of the results obtained in
Sec. I.5.

###  I.1 Preliminaries

####  I.1.1 Lattice vectors

Both graphene and the TMDC layer have hexagonal lattices. We define the
lattice vectors $\mathbf{a}_{1}=\frac{a}{2}\left(1,\sqrt{3}\right)^{T}$,
$\mathbf{a}_{2}=\frac{a}{2}\left(1,-\sqrt{3}\right))^{T}$ (see Fig. 1(a)) and
the corresponding primitive reciprocal lattice vectors
$\mathbf{b}_{1}=\frac{2\pi}{a}\left(1,1/\sqrt{3}\right)^{T}$,
$\mathbf{b}_{2}=\frac{2\pi}{a}\left(1,-1/\sqrt{3}\right))^{T}$. Here
$a=a_{gr}$ ($a=a_{tmdc}$) is the lattice constant for graphene (TMDC). We will
distinguish the reciprocal lattice vectors $\mathbf{b}_{1,2}$ of graphene from
the TMDC ones by using the notation $\mathbf{b}^{\prime}_{1,2}$ for the
latter. The position of the $A$ and $B$ sublattice in the unit cell is given
by $\mathbf{t}_{A}=\frac{a_{gr}}{2}\left(1,1/\sqrt{3}\right)^{T}$,
$\mathbf{t}_{B}=\frac{a_{gr}}{2}\left(1,-1/\sqrt{3}\right)^{T}.$ The metal
(chalcogen) atoms occupy a position corresponding to the $A$ ($B$) sublattice
in the unit cell of the TMDC.

####  I.1.2 Formalism to calculate the induced SOC

As it was shown in Refdavid_induced_2019 , an important contribution to the
proximity induced Bychkov-Rashba type SOC in graphene comes from virtual
interlayer tunneling processes to the TMDC layer. In third order perturbation
theory it is given by

$(H^{gr}_{R})_{Xs,X^{\prime}s^{\prime}}=\sum_{j,b,b^{\prime},s^{\prime\prime},s^{\prime\prime\prime}}\frac{\left(T(\tau\mathbf{k}^{\prime}_{j})\right)_{Xs,bs^{\prime\prime}}(H_{\text{soc}})_{bs^{\prime\prime},b^{\prime}s^{\prime\prime\prime}}\left(T^{\dagger}(\tau\mathbf{k}^{\prime}_{j})\right)_{b^{\prime}s^{\prime\prime\prime},X^{\prime}s^{\prime}}}{[E_{D}^{\text{gr}}-E_{b}^{\text{tmdc}}(\tau\mathbf{k}^{\prime}_{j})][E_{D}^{\text{gr}}-E_{b^{\prime}}^{\text{tmdc}}(\tau\mathbf{k}^{\prime}_{j})]}.$
(1)

Here $X=\\{A,B\\}$ runs over sublattice indices of graphene, $s$, $s^{\prime}$
are spin indices and $b\neq b^{\prime}$ are band indices running over the
bands of the (isolated) TMDC layer. Spin-flip process can be facilitated by
certain off-diagonal elements of the intrinsic SOC matrix of the monolayer
TMDC denoted by $(H_{\text{soc}})_{bs,b^{\prime}s^{\prime}}$ The spin-flip
off-diagonal matrix elements are allowed between pairs of bands $b$,
$b^{\prime}$ if one of the bands is symmetric (even) and the other one is
antisymmetric (odd) with respect to reflection on the horizontal mirror plane
of the TMDC (see, e.g., Ref. kormanyos_spin-orbit_2014  for further discussion
of the SOC in monolayer TMDCs).

We assume that the graphene layer is rotated with respect to the TMDC layer by
a twist angle $\theta$. Therefore, as indicated in Fig. 1(b), the BZ of
graphene is also rotated. For each interlayer twist angle $\theta$ the right-
hand side of Eq. (1) needs to be evaluated for wavenumbers
$\tau\mathbf{k}^{\prime}_{j}$ in the TMDC Brillouin zone (see Fig. 1(b)) that
are determined by quasimomentum conservation condition
$\tau\mathbf{K}^{\theta}+\mathbf{G}^{\theta}=\mathbf{k}^{\prime}+\mathbf{G}^{\prime}$,
such that
$|\tau\mathbf{K}^{\theta}+\mathbf{G}^{\theta}|=|\mathbf{K}^{\theta}|$. Here
$\tau\mathbf{K}^{\theta}$, $\tau=\pm 1$ denote the wavevector to the Dirac
point of graphene and $\mathbf{G}^{\theta}$ and ($\mathbf{G}^{\prime}$) are
reciprocal lattice vectors of graphene (TMDC). As indicated in Fig. 1(b), for
each $\tau\mathbf{K}^{\theta}$ there are three wavevectors
$\tau\mathbf{k}^{\prime}_{j}$ and corresponding reciprocal lattice vectors
$\mathbf{G}^{\theta}_{j}$ and $\mathbf{G}^{\prime}_{j}$ that satisfy the
quasimomentum conservation condition. These $\tau\mathbf{k}^{\prime}_{j}$
vectors are related by $2\pi/3$ rotations. Explicitly, one finds
$\mathbf{G}^{\theta}_{1,2,3}=-\mathbf{b}_{2}^{\theta},-\mathbf{b}_{1}^{\theta},0$
and
$\mathbf{G}^{\prime}_{1,2,3}=-\mathbf{b}^{\prime}_{2},-(\mathbf{b}^{\prime}_{1}-\mathbf{b}^{\prime}_{2}),\mathbf{b}^{\prime}_{1}$.
$E_{b}^{\text{tmdc}}(\tau\mathbf{k}^{\prime}_{j})$ denote TMDC bands energies
and $E_{D}^{\text{gr}}$ is the energy of the Dirac point of graphene.
$(T_{\tau\mathbf{k}^{\prime}_{j}})_{Xs,bs^{\prime\prime}}$ is the tunneling
matrix between Bloch states at $\tau\mathbf{K}^{\theta}$ to states at
$\tau\mathbf{k}^{\prime}_{j}$. We assume spin-independent tunneling and
therefore $(T_{\tau\mathbf{k}^{\prime}_{j}})_{Xs,bs^{\prime\prime}}$ is
diagonal in the spin space.
$(T_{\tau\mathbf{k}^{\prime}_{j}})_{Xs,bs^{\prime\prime}}$ can be written as

$\displaystyle\left(T(\tau\bm{k}^{\prime}_{j})\right)_{Xb}=\langle\Psi^{(X)}_{gr}({\mathbf{K}}_{1}^{\theta},\mathbf{r})|H_{orb}|\Psi^{(b)}_{tmdc}(\mathbf{k}^{\prime}_{j},\mathbf{r})\rangle=e^{i\tau\bm{G}^{\theta}_{j}\cdot\bm{d}_{0}}t_{b}(\tau\mathbf{k}^{\prime}_{j})e^{i\tau\phi_{j}^{X}}$
(2)

Here $H_{orb}$ is the orbital part of the Hamiltonian of the system (without
the SOC), $\Psi^{(X)}_{gr}({\mathbf{K}}_{1}^{\theta},\mathbf{r})$ and
$\Psi^{(b)}_{tmdc}(\mathbf{k}^{\prime}_{j},\mathbf{r})$ are Bloch
wavefunctions of graphene and of band $b$ of the TMDC layer, respectively. A
possible lateral shift of the graphene lattice with respect to the TMDC
lattice is denoted by $\mathbf{d}_{0}$. The phase $\phi_{j}^{X}$ is defined as
$\phi_{j}^{X}=\bm{t}_{X}\cdot\mathbf{G}_{j}$. The tunneling amplitude
$t_{b}(\tau\mathbf{k}^{\prime}_{j})$ to band $b$ of the TMDC can be written as

$t_{b}(\tau\mathbf{k}^{\prime}_{j})=e^{-i\tau\bm{G}^{\prime}_{j}\cdot\bm{t}^{\prime}_{X}}t_{b}^{p}\left(\tau\mathbf{k}^{\prime}_{j}\right)+e^{-i\tau\bm{G}^{\prime}_{j}\cdot\bm{t}^{\prime}_{X^{\prime}}}t^{d}_{b}\left(\tau\mathbf{k}^{\prime}_{j}\right),$
(3)

where $\bm{t}^{\prime}_{X}$ ($\bm{t}^{\prime}_{X^{\prime}}$), $X\neq
X^{\prime}$ is the vector pointing to the chalcogen (metal) atom position in
the unit cell of the TMDC layer, and $t^{p}_{b}(\tau\mathbf{k}^{\prime}_{j})$
($t^{d}_{b}((\tau\mathbf{k}^{\prime}_{j})$) describes the interlayer tunneling
between graphene and the $p$ orbitals of the chalcogen ($d$ orbitals of the
metal) atoms. The calculation of the amplitude
$t_{b}^{p,d}(\tau\mathbf{k}^{\prime}_{1})$ in terms of materials dependent
parameters will be discussed in Sec. II.1.

Figure 1: a) Real space orientation of the hexagonal lattice and the primitive
lattice vectors $\mathbf{a}_{1}$, $\mathbf{a}_{2}$. b) The rotated BZ and
Dirac points $\mathbf{K}^{\theta}$ of graphene, along with the wavenumbers
$\mathbf{k}^{\prime}_{j}$ in the BZ of the TMDC. As $\theta$ changes in the
range $[0..\pi/3]$, the $\mathbf{k}^{\prime}_{j}$ trace out the arcs shown by
dashed lines. c) Reflection planes of the hexagonal BZ.
$\mathcal{R}_{v}^{(1,2,3)}$ is relevant for both graphene and the TMDC,
$\mathcal{R}_{d}^{(1,2,3)}$ only for graphene.

$(H_{\text{soc}})_{bs,b^{\prime}s^{\prime}}$ denote the matrix elements of the
spin-orbit coupling operator $\hat{H}_{soc}$ between Bloch wavefunctions
corresponding to bands $b$ and $b^{\prime}$ of the monolayer TMDC. For our
purposes the matrix elements between an even
$\Psi^{(e)}_{tmdc}(\tau\mathbf{k}^{\prime}_{j},\mathbf{r})$ and an odd
$\Psi^{(o)}_{tmdc}(\tau\mathbf{k}^{\prime}_{j},\mathbf{r})$ Bloch states are
important and they can be written as

$H_{e,o}(\tau\mathbf{k}^{\prime}_{j})=\langle\Psi^{(e)}_{tmdc}(\tau\mathbf{k}^{\prime}_{j},\mathbf{r})|\hat{H}_{soc}|\Psi^{(o)}_{tmdc}(\tau\mathbf{k}^{\prime}_{j},\mathbf{r})\rangle=i\lambda_{soc}\left[\alpha_{e,o}^{(x)}(\tau\mathbf{k}^{\prime}_{j}){s}_{x}+\alpha_{e,o}^{(y)}(\tau\mathbf{k}^{\prime}_{j}){s}_{y}\right].$
(4)

Here $\lambda_{soc}$ characterizes the strength of the intrinsic SOC in the
TMDC layer.
$\hat{H}_{soc}\sim\hat{L}_{x}\hat{S}_{x}+\hat{L}_{y}\hat{S}_{y}+\hat{L}_{z}\hat{S}_{z}=\frac{1}{2}(\hat{L}_{+}\hat{S}_{-}+\hat{L}_{-}\hat{S}_{+})+\hat{L}_{z}\hat{S}_{z}$,
$\hat{L}_{x,y,z}$ are angular momentum operators,
$\hat{L}_{\pm}=\hat{L}_{x}\pm i\hat{L}_{y}$,
$\hat{S}_{x,y,z}=\frac{\hbar}{2}s_{x,y,z}$ are spin operators and the Pauli
matrices $s_{x,y,z}$ act on the spin degree of freedom. The complex numbers
$\alpha_{e,o}^{(x,y)}(\tau\mathbf{k}^{\prime}_{j})$ are defined as

$i\,\alpha_{e,o}^{(x)}(\tau\mathbf{k}^{\prime}_{j})=\langle\Psi^{(e)}_{tmdc}(\tau\mathbf{k}^{\prime}_{j},\mathbf{r})\left|\hat{L}_{x}\right|\Psi^{(o)}_{tmdc}(\tau\mathbf{k}^{\prime}_{j},\mathbf{r})\rangle,$
(5a)
$i\,\alpha_{e,o}^{(x)}(\tau\mathbf{k}^{\prime}_{j})=\langle\Psi^{(e)}_{tmdc}(\tau\mathbf{k}^{\prime}_{j},\mathbf{r})\left|\hat{L}_{y}\right|\Psi^{(o)}_{tmdc}(\tau\mathbf{k}^{\prime}_{j},\mathbf{r})\rangle.$
(5b)

One can make use of the fact that the $\tau\mathbf{k}^{\prime}_{j}$ vectors
are related by $2\pi/3$ rotations, which is also a symmetry of the TMDC
lattice. Following Ref. dresselhaus_group_2010  one may write
$\Psi^{(e(o))}_{tmdc}(\mathcal{R}_{\pm
2\pi/3}\tau\mathbf{k}^{\prime}_{j},\mathbf{r})=\mathcal{R}_{\pm
2\pi/3}\Psi^{(e(o))}_{tmdc}(\tau\mathbf{k}^{\prime}_{j},\mathbf{r})$ where
$\mathcal{R}_{\pm 2\pi/3}$ denotes rotation by $\pm 2\pi/3$. Substituting this
into Eq. (4) and taking into account that
$\mathcal{R}_{2\pi/3}\hat{L}_{\pm}(\mathcal{R}_{2\pi/3})^{\dagger}=e^{\mp
i2\pi/3}L_{\pm}$, one finds

$H_{e,o}(\tau\mathbf{k}^{\prime}_{j})=i\lambda_{soc}\left(\begin{array}[]{cc}0&\alpha_{e,o}^{(-)}(\tau\mathbf{k}^{\prime}_{1})e^{-i(j-1)2\pi/3}\\\
\alpha_{e,o}^{(+)}(\tau\mathbf{k}^{\prime}_{1})e^{i(j-1)2\pi/3}&0\end{array}\right),$
(6)

where
$\alpha_{e,o}^{(\pm)}(\mathbf{k}^{\prime}_{1})=\alpha_{e,o}^{(x)}(\tau\mathbf{k}^{\prime}_{1})\pm
i\alpha_{e,o}^{(y)}(\tau\mathbf{k}^{\prime}_{1})$.

###  I.2 Derivation of $\vartheta_{R}$ and $\lambda_{R}$

It is instructive to consider the contribution of a single pair of $b=e$ and
$b^{\prime}=o$ bands in Eq. (1). Using Eq.(2), one finds

$(H^{gr}_{R})_{e,o}=\sum_{j}\left(\begin{array}[]{cc}1&e^{i\tau(\phi_{j}^{A}-\phi_{j}^{B})}\\\
e^{-i\tau(\phi_{j}^{A}-\phi_{j}^{B})}&1\end{array}\right)\otimes\frac{\left[{t}_{e,o}(\tau\mathbf{k}^{\prime}_{j})H_{e,o}(\tau\mathbf{k}^{\prime}_{j})+{t}_{o,e}(\tau\mathbf{k}^{\prime}_{j})H_{o,e}(\tau\mathbf{k}^{\prime}_{j})\right]}{[E_{D}^{\text{gr}}-E_{e}^{\text{tmdc}}(\tau\mathbf{k}^{\prime}_{j})][E_{D}^{\text{gr}}-E_{o}^{\text{tmdc}}(\tau\mathbf{k}^{\prime}_{j})]}\\\
=\sum_{j}\mathbf{M}_{j}(\tau\phi_{j}^{A},\tau\phi_{j}^{B})\otimes\frac{\left[{t}_{e,o}(\tau\mathbf{k}^{\prime}_{j})H_{e,o}(\tau\mathbf{k}^{\prime}_{j})+{t}_{o,e}(\tau\mathbf{k}^{\prime}_{j})H_{o,e}(\tau\mathbf{k}^{\prime}_{j})\right]}{[E_{D}^{\text{gr}}-E_{e}^{\text{tmdc}}(\tau\mathbf{k}^{\prime}_{j})][E_{D}^{\text{gr}}-E_{o}^{\text{tmdc}}(\tau\mathbf{k}^{\prime}_{j})]}$
(7)

As explained in Sec. I.1, the phases $\phi_{j}^{A,B}$ in
$\mathbf{M}_{j}(\tau\phi_{j}^{A},\tau\phi_{j}^{B})$ are defined as
$\phi_{j}^{A,B}=\bm{t}_{A,B}\cdot\mathbf{G}_{j}$. There are a number of
possible choices for the vectors $\bm{t}_{A,B}$ and for the primitive
reciprocal lattice vectors $\mathbf{b}_{1,2}$ of graphene, and the latter
determine $\mathbf{G}_{j}$. For a given $j$ the phase $\phi_{j}^{A,B}$ depends
on these choices, but in general $\phi_{j}^{A,B}\in\\{0,2\pi/3,-2\pi/3\\}$.
The different choices lead to unitary equivalent
$\mathbf{M}^{\prime}_{j}(\tau\phi_{j}^{A},\tau\phi_{j}^{B})=\mathbf{U}(\tau)\mathbf{M}_{j}(\tau\phi_{j}^{A},\tau\phi_{j}^{B})\mathbf{U}^{\dagger}(\tau)$
matrices and hence unitary equivalent $(H^{gr}_{R})_{e,o}$ Hamiltonians.
Importantly, $\mathbf{U}(\tau)$ does not depend on the geometric twist angle
$\theta$. Using Sec. I.1.1, the $\mathbf{M}_{j}$ matrices read

$\mathbf{M}_{j}(\tau\phi_{j}^{A},\tau\phi_{j}^{B})=\left(\begin{array}[]{cc}1&e^{i\tau\frac{2\pi}{3}j}\\\
e^{-i\tau\frac{2\pi}{3}j}&1\end{array}\right).$ (8)

Furthermore,
$H_{o,e}(\tau\mathbf{k}^{\prime}_{j})=[H_{e,o}(\tau\mathbf{k}^{\prime}_{j})]^{\dagger}$,
${t}_{e,o}(\tau\mathbf{k}^{\prime}_{j})={t}_{e}(\tau\mathbf{k}^{\prime}_{j}){t}_{o}^{*}(\tau\mathbf{k}^{\prime}_{j})=(-1)^{\mu}|{t}_{e,o}(\tau\mathbf{k}^{\prime}_{j})|e^{i\eta(\tau\mathbf{k}^{\prime}_{j})}$
and
${t}_{o,e}(\tau\mathbf{k}^{\prime}_{j})={t}_{e,o}^{*}(\tau\mathbf{k}^{\prime}_{j})$.
We define $\eta(\tau\mathbf{k}^{\prime}_{j})={\rm
Arg}[{t}_{e,o}(\tau\mathbf{k}^{\prime}_{j})]\,{\rm mod}\,\pi$ and, as it will
be explained in more detail in Sec. I.3, the index $\mu=0,1$ gives the sign of
the induced Rashba SOC which depends on the layer stacking. Note that the
phase factors $e^{i\tau\bm{G}^{\theta}_{j}\cdot\bm{d}_{0}}$, which encode the
effect of a lateral shift between the layers in Eq. (2), have dropped out of
the calculations. This means that our results for the induced Rashba SOC do
not depend on such shift. The same conclusion holds for the valley Zeeman SOC
as well, which depends on $|t_{b}(\tau\mathbf{k}^{\prime}_{j})|^{2}$, see Ref.
david_induced_2019, .

Because of the $C_{3}$ symmetry of the system,
${t}_{e,o}(\tau\mathbf{k}^{\prime}_{j})$ does not depend on the index $j$ and
from time reversal symmetry follows that
$t_{e,o}(-\mathbf{k}^{\prime}_{1})=t_{e,o}^{*}(\mathbf{k}^{\prime}_{1})$.
Therefore, by taking into account Eq. (6),

$\displaystyle{t}_{e,o}(\tau\mathbf{k}^{\prime}_{j})H_{e,o}(\tau\mathbf{k}^{\prime}_{j})+{t}_{e,o}^{*}(\tau\mathbf{k}^{\prime}_{j})[H_{e,o}(\tau\mathbf{k}^{\prime}_{j})]^{\dagger}=2i\lambda_{soc}|{t}_{e,o}(\tau\mathbf{k}^{\prime}_{1})|\left(\begin{array}[]{cc}0&\Lambda_{e,o}(\tau\mathbf{k}^{\prime}_{1})e^{-i(j-1)2\pi/3}\\\
-\Lambda_{e,o}^{*}(\tau\mathbf{k}^{\prime}_{1})e^{i(j-1)2\pi/3}&0\end{array}\right),$
(11)

where
$\Lambda_{e,o}(\tau\mathbf{k}^{\prime}_{1})=Im[e^{i\eta(\tau\mathbf{k}^{\prime}_{1})}\alpha_{e,o}^{(y)}(\tau\mathbf{k}^{\prime}_{1})]+iIm[e^{i\eta(\tau\mathbf{k}^{\prime}_{1})}\alpha_{e,o}^{(x)}(\tau\mathbf{k}^{\prime}_{1})]$.
One can show that
$\Lambda_{e,o}(\tau\mathbf{k}^{\prime}_{1})=\tau\Lambda_{e,o}(\mathbf{k}^{\prime}_{1})$
and one may write
$\Lambda_{e,o}(\mathbf{k}^{\prime}_{1})=|\Lambda_{e,o}(\mathbf{k}^{\prime}_{1})|e^{i\tilde{\vartheta}_{e,o}(\mathbf{k}^{\prime}_{1})}$,
where

$\tilde{\vartheta}_{e,o}(\mathbf{k}^{\prime}_{1})={\rm
Arg}[\Lambda_{e,o}(\mathbf{k}^{\prime}_{1})]={\rm Arctan}\left[\frac{{\rm
Im}[e^{i\eta(\mathbf{k}^{\prime}_{1})}\alpha_{e,o}^{(x)}(\mathbf{k}^{\prime}_{1})]}{{\rm
Im}[e^{i\eta(\mathbf{k}^{\prime}_{1})}\alpha_{e,o}^{(y)}(\mathbf{k}^{\prime}_{1})]}\right]+m\pi,\hskip
14.22636ptm=0,\pm 1,\pm 2\dots$ (12)

Here ${\rm Im}[\dots]$ denotes the imaginary part.
$\tilde{\vartheta}_{e,o}(\mathbf{k}^{\prime}_{1})$ is a quantum phase which
depends on i) the phase $\eta(\mathbf{k}^{\prime}_{1})$ of the interlayer
tunneling amplitude between the Bloch states of graphene and the TMDC layer
and ii) on off-diagonal matrix elements of the intrinsic SOC of the TMDC layer
defined in Eq.(4).

Due to $C_{3}$ symmetry, the denominator of Eq. (7) is the same for all $j$
and because of time reversal symmetry it does not depend on $\tau$. Therefore
the denominator can be pulled out of the summation. Substituting Eqs. (8) and
(11) into Eq. (7) and performing the summation over $j$ one finds

$(H^{gr}_{R})_{e,o}=\lambda_{R,(e,o)}(\mathbf{k}^{\prime}_{1})e^{i\frac{s_{z}}{2}\vartheta_{e,o}(\mathbf{k}^{\prime}_{1})}\left(\begin{array}[]{cccc}0&0&0&i\frac{\tau+1}{2}\\\
0&0&-i\frac{\tau-1}{2}&0\\\ 0&i\frac{\tau-1}{2}&0&0\\\
-i\frac{\tau+1}{2}&0&0&0\end{array}\right)e^{-i\frac{s_{z}}{2}\vartheta_{e,o}(\mathbf{k}^{\prime}_{1})}.$
(13)

Here we defined

${\vartheta}_{e,o}(\mathbf{k}^{\prime}_{1})=\tilde{\vartheta}_{e,o}(\mathbf{k}^{\prime}_{1})+2\pi/3,$
(14)

and $\lambda_{R,(e,o)}(\mathbf{k}^{\prime}_{1})$ is given bydavid_induced_2019

$\lambda_{R,(e,o)}(\mathbf{k}^{\prime}_{1})=\frac{(-1)^{\mu}6\lambda_{soc}|{t}_{e,o}(\mathbf{k}^{\prime}_{1})||\Lambda_{e,o}(\mathbf{k}^{\prime}_{1})|}{\left[E^{\text{gr}}_{D}-E_{e}^{\text{tmdc}}(\mathbf{k}^{\prime}_{1})\right]\left[E^{\text{gr}}_{D}-E_{o}^{\text{tmdc}}(\mathbf{k}^{\prime}_{1})\right]}.$
(15)

Using the $\tau_{z}$ Pauli matrix, one may rewrite Eq. (13) in a more compact
form:

$(H^{gr}_{R})_{e,o}=-\frac{\lambda_{R,(e,o)}(\mathbf{k}^{\prime}_{1})}{2}e^{i\frac{s_{z}}{2}{\vartheta}_{e,o}(\mathbf{k}^{\prime}_{1})}\left(\tau_{z}\sigma_{x}s_{y}+\sigma_{y}s_{x}\right)e^{-i\frac{s_{z}}{2}{{\vartheta}_{e,o}(\mathbf{k}^{\prime}_{1})}}.$
(16)

Note, that in the derivation of Eqs. (14)-(16) we only made use of time
reversal symmetry and that the $\mathbf{k}^{\prime}_{j}$ vectors are related
by $2\pi/3$ rotations. Therefore these results are valid for arbitrary
interlayer twist angle $\theta$. On the other hand, we considered the
contribution of a single pair of $e$ and $o$ bands of the TMDC to the induced
Rashba type SOC. As Eq. (1) shows, one needs to sum over all $q$ pairs of
$(e,o)$ bands and one may define the complex Rashba SOC coefficient by

$\bar{\lambda}_{R}=\sum_{q}\lambda_{\text{R},(e,o)_{q}}e^{i\vartheta_{(e,o)_{q}}}.$
(17)

In terms of its magnitude
$\lambda_{R}(\mathbf{k}^{\prime}_{1})=|\bar{\lambda}_{R}(\mathbf{k}^{\prime}_{1})|$
and phase
$\vartheta_{R}(\mathbf{k}^{\prime}_{1})=\mathrm{Arg}[\bar{\lambda}_{R}(\mathbf{k}^{\prime}_{1})]$,
the Hamiltonian of the induced Rashba SOC can be written as

$H^{gr}_{R}=-\frac{\lambda_{R}(\mathbf{k}^{\prime}_{1})}{2}e^{i\frac{s_{z}}{2}{\vartheta}_{R}(\mathbf{k}^{\prime}_{1})}\left(\tau_{z}\sigma_{x}s_{y}+\sigma_{y}s_{x}\right)e^{-i\frac{s_{z}}{2}{{\vartheta}_{R}(\mathbf{k}^{\prime}_{1})}}.$
(18)

Note, that both $\lambda_{R}(\mathbf{k}^{\prime}_{1})$ and
${\vartheta}_{R}(\mathbf{k}^{\prime}_{1})$ depend implicitly on $\theta$
through $\mathbf{k}^{\prime}_{1}$, as shown in Fig. 1. We will discuss the
$\theta$ dependence of $\vartheta_{R}$ in more detail in Sec. I.5.

###  I.3 Staking order dependent sign of $\lambda_{R}$

In our theory the sign of $\lambda_{R}$ depends on the stacking order of the
graphene and the TMDC layers, i.e., whether is graphene below or above the
TMDC layer. Such a sign has no physical significance in TMDC/graphene bilayer
stacks, but it is important in TMDC/graphene/TMDC trilayers, because it means
that the contributions of the two TMDC layers to the induced Rashba SOC add up
with an opposite sign.

As explained below Eq. (8), the induced Rashba SOC is a second order process
in interlayer tunneling. The tunneling amplitude is given by
${t}_{e,o}(\tau\mathbf{k}^{\prime}_{j})={t}_{e}(\tau\mathbf{k}^{\prime}_{j}){t}_{o}^{*}(\tau\mathbf{k}^{\prime}_{j})$,
i.e., it depends on the product of the tunneling amplitudes to an $e$ and to
an $o$ band of the TMDC. $\Psi^{(e)}_{tmdc}(\mathbf{r})$ of an $e$ band of the
TMDC is even with respect to the mirror plane that contains the metal atoms of
the TMDC layer. However, the graphene $p_{z}$ orbitals have a different sign
above and below the graphene layer. Therefore one can expect that the matrix
element given in Eq. (2) between the graphene Bloch wavefunction
$\Psi_{gr}(\mathbf{r})$ and $\Psi^{(e)}_{tmdc}(\mathbf{r})$ should have a
global sign difference depending on the stacking order of the two layers,
i.e., whether the graphene layer is above or below the TMDC layer. On the
other hand, similar consideration suggests that the sign of the amplitude
${t}_{o}(\tau\mathbf{k}^{\prime}_{j})$ between $\Psi_{gr}(\mathbf{r})$ and
$\Psi^{(o)}_{tmdc}(\mathbf{r})$ does not depend on the stacking order.
Therefore, there will be a stacking order dependent sign factor in
${t}_{e,o}(\tau\mathbf{k}^{\prime}_{j})={t}_{e}(\tau\mathbf{k}^{\prime}_{j}){t}_{o}^{*}(\tau\mathbf{k}^{\prime}_{j})$,
which is described by the index $\mu=0,1$ introduced below Eq. (8). The
stacking order dependence of ${t}_{e}(\tau\mathbf{k}^{\prime}_{j})$ and
${t}_{o}(\tau\mathbf{k}^{\prime}_{j})$ can be explicitly shown if they are
calculated e.g., using the two-center Slater-Koster parametrization of the
corresponding transfer integrals, see Refs.david_induced_2019 ; li_twist-
angle_2019  for details.

Because of this stacking order dependent sign, the induced Rashba type SOC in
TMDC/graphene/TMDC trilayers, which is given by the sum of the contributions
from the two TMDC layer, reads $H_{R}^{(tls)}=H_{R}^{(b)}-H_{R}^{(t)}$. Taking
into account the interlayer twist angle dependence of $\lambda_{R}$ (see Sec.
II.3) and of $\theta+\vartheta_{R}$ (shown in Fig.1 of the main text) one can
then easily see that $H_{R}^{(tls)}$ vanishes whenever the TMDC/graphene/TMDC
stack is either inversion symmetric or has a horizontal mirror symmetry.

###  I.4 Effects of the interlayer twist

$H^{gr}_{R}$ in Eq. (18) was obtained in the $x^{\prime}-y^{\prime}$
coordinate system of the TMDC layer (see Fig. 2), where the matrix elements of
$\hat{H}_{soc}$ were evaluated.

Figure 2: Schematics of a twisted graphene/TMDC bilayer with twist angle
$\theta$ and the coordinate systems $x-y$ and $x^{\prime}-y^{\prime}$.

Our choice of reciprocal lattice vectors leads to
$H_{orb}^{gr}=v_{F}(\tau_{z}\sigma_{x}\delta k_{x}-\sigma_{y}\delta k_{y})$
for the orbital Hamiltonian of graphene in the $x-y$ coordinate system, where
$\delta\mathbf{k}=(\delta k_{x},\delta k_{y})$ is a small wavevector measured
from the $\tau\mathbf{K}$ point. In the $x^{\prime}-y^{\prime}$ system
$H_{orb}^{gr}$ is given by
$H_{orb}^{gr,\theta}=e^{-i\tau_{z}\frac{\sigma_{z}}{2}{\theta}}H_{orb}^{gr}e^{i\tau_{z}\frac{\sigma_{z}}{2}{\theta}}$.
The total effective Hamiltonian of graphene is therefore

$H_{\rm
eff}^{gr,\theta}=e^{-i\tau_{z}\frac{\sigma_{z}}{2}{\theta}}H_{orb}^{gr}e^{i\tau_{z}\frac{\sigma_{z}}{2}{\theta}}+H^{gr}_{R}+H^{gr}_{vZ},$
(19)

where we have also taken into account the valley-Zeeman induced SOC, which is
described by

$H_{vZ}^{gr}=\lambda_{vZ}\tau_{z}s_{z}.$ (20)

$H_{vZ}^{gr}$ is invariant under rotations around the $\hat{z}$ axis and has
the same form in both the $x-y$ and $x^{\prime}-y^{\prime}$ coordinate
systems. One may perform a unitary transformation
$e^{i\tau_{z}\frac{\sigma_{z}}{2}{\theta}}H_{\rm
eff}^{gr,\theta}e^{-i\tau_{z}\frac{\sigma_{z}}{2}{\theta}}$. This
transformation leaves $H_{vZ}$ unchanged and one finds

$H_{\rm
eff}^{gr}=H_{orb}^{gr}+e^{i\tau_{z}\frac{\sigma_{z}}{2}{\theta}}H_{R}^{gr}e^{-i\tau_{z}\frac{\sigma_{z}}{2}{\theta}}+H_{vZ}^{gr}.$
(21)

It is convenient to perform another unitary transformation $H_{\rm
eff}^{gr}\rightarrow\sigma_{x}H_{\rm eff}^{gr}\sigma_{x}$. This changes the
orbital Hamiltonian $v_{F}(\tau_{z}\sigma_{x}\delta k_{x}-\sigma_{y}\delta
k_{y})\rightarrow v_{F}(\tau_{z}\sigma_{x}\delta k_{x}+\sigma_{y}\delta
k_{y}),$ leaves $H_{vZ}^{gr}$ unchanged, and transforms the Rashba Hamiltonian
to

$\sigma_{x}e^{i\tau_{z}\frac{\sigma_{z}}{2}{\theta}}H_{R}^{gr}e^{-i\tau_{z}\frac{\sigma_{z}}{2}{\theta}}\sigma_{x}=-\frac{\lambda_{R}(\mathbf{k}^{\prime}_{1})}{2}e^{-i\tau_{z}\frac{\sigma_{z}}{2}{\theta}}e^{i\frac{s_{z}}{2}{\vartheta}_{R}(\mathbf{k}^{\prime}_{1})}\left(\tau_{z}\sigma_{x}s_{y}-\sigma_{y}s_{x}\right)e^{-i\frac{s_{z}}{2}{\vartheta}_{R}(\mathbf{k}^{\prime}_{1})}e^{i\tau_{z}\frac{\sigma_{z}}{2}{\theta}}.$
(22)

Eq. (22) shows more explicitly the effects of interlayer rotation on the
graphene Rashba Hamiltonian
$\tilde{H}^{gr}_{R}=\frac{\lambda_{R}}{2}\left(\tau_{z}\sigma_{x}s_{y}-\sigma_{y}s_{x}\right)$
derived previously in Refs. kane_quantum-spin-Hall_2005 ; min_intrinsic_2006 :
i) $\tilde{H}^{gr}_{R}$ is rotated in real space by $\theta$, an ii) if
$\vartheta_{R}\neq 2m\pi$, $m$ integer, then a rotation in spin-space appears.

###  I.5 Calculation of $\vartheta_{R}$ for $\theta=0$, $\theta=\pi/6$ and
$\theta=\pi/3$ twist angles

The calculations of Sec. I.2 are valid for arbitrary interlayer twist angle
$\theta$. For $\theta_{l}=l\pi/6$, $l=0,1,2\dots$ the graphene/TMDC bilayer
has a higher, $C_{3v}$ symmetry, which implies further constraints on the form
of ${H}^{gr}_{R}$. By calculating $\vartheta_{e,o}(\mathbf{k}^{\prime}_{1})$
defined in Eq. (14) for $\theta=0,\pi/6,\pi/3$ explicitly, we show how the
results of Sec. I.2 are simplified in this case.

####  I.5.1 $\theta=0$ twist angle

In Fig. 3 we show the wavevectors $\mathbf{k}^{\prime}_{1,2,3}$ that satisfy
the quasimomentum conservation as well as the three reflection planes
$\mathcal{R}_{v}^{(1,2,3)}$ that the system possess in addition to the $C_{3}$
symmetry. First, we will derive a useful relation for the tunneling amplitude
$t_{e,o}(\mathbf{k}^{\prime}_{1})=t_{e}(\mathbf{k}^{\prime}_{1})t_{o}^{*}(\mathbf{k}^{\prime}_{1})$,
defined below Eq. (8).

Figure 3: The BZ of graphene (purple) and of the TMDC (orange), along with
$\mathbf{k}^{\prime}_{1,2,3}$ vectors and the reflection planes
$\mathcal{R}_{v}^{1,2,3}$ for $\theta=0$.

Using the definition $t_{e(o)}(\tau\mathbf{k}^{\prime}_{j})$ given in Eq.(2)
and denoting $\mathbf{K}^{\theta=0}=\mathbf{K}^{(0)}$, consider the following:

$\displaystyle t_{e(o)}(\mathbf{k}^{\prime}_{1})e^{-i\phi^{X}_{1}}$
$\displaystyle=\langle\Psi^{(X)}_{gr}({\mathbf{K}}_{1}^{(0)},\mathbf{r})|H_{orb}|\Psi^{(e(o))}_{tmdc}(\mathbf{k}^{\prime}_{1},\mathbf{r})\rangle=\langle(H_{orb}\Psi^{(e(o))}_{tmdc}(\mathbf{k}^{\prime}_{1},\mathbf{r}))^{*}|(\Psi^{(X)}_{gr}({\mathbf{K}}_{1}^{(0)},\mathbf{r}))^{*}\rangle$
$\displaystyle=\langle
H_{orb}(\Psi^{(e(o))}_{tmdc}(\mathbf{k}^{\prime}_{1},\mathbf{r}))^{*}|(\Psi^{(X)}_{gr}({\mathbf{K}}_{1}^{(0)},\mathbf{r}))^{*}\rangle=\langle
H_{orb}\Psi^{(e(o))}_{tmdc}(-\mathbf{k}^{\prime}_{1},\mathbf{r})|\Psi^{(X)}_{gr}({-\mathbf{K}}_{1}^{(0)},\mathbf{r})\rangle$
$\displaystyle=\langle
H_{orb}\Psi^{(e(o))}_{tmdc}(\mathcal{R}_{v}^{(2)}(\mathbf{k}^{\prime}_{1}-\mathbf{b}^{\prime}_{1}),\mathbf{r})|\Psi^{(X)}_{gr}(\mathcal{R}_{v}^{(1)}\mathbf{K}_{1}^{(0)},\mathbf{r})\rangle$
$\displaystyle=\langle
H_{orb}\mathcal{R}_{v}^{(2)}\Psi^{(e(o))}_{tmdc}(\mathbf{k}^{\prime}_{1}-{\mathbf{b}^{\prime}_{1}},\mathbf{r})|\Psi^{(X)}_{gr}(\mathcal{R}_{v}^{(1)}\mathbf{K}_{1}^{(0)},\mathbf{r})\rangle$
$\displaystyle=\langle\mathcal{R}_{v}^{(2)}H_{orb}\Psi^{(e(o))}_{tmdc}(\mathbf{k}^{\prime}_{1}-{\mathbf{b}^{\prime}_{1}},\mathbf{r})|\Psi^{(X)}_{gr}(\mathcal{R}_{v}^{(1)}\mathbf{K}_{1}^{(0)},\mathbf{r})\rangle$
$\displaystyle=\langle
H_{orb}\Psi^{(e(o))}_{tmdc}(\mathbf{k}^{\prime}_{1}-{\mathbf{b}^{\prime}_{1}},\mathbf{r})|\mathcal{R}_{v}^{(2)}\Psi^{(X)}_{gr}(\mathcal{R}_{v}^{(1)}\mathbf{K}_{1}^{(0)},\mathbf{r})\rangle$
$\displaystyle=\langle
H_{orb}\Psi^{(e(o))}_{tmdc}(\mathbf{k}^{\prime}_{1}-{\mathbf{b}^{\prime}_{1}},\mathbf{r})|\Psi^{(X)}_{gr}(\mathcal{R}_{v}^{(2)}\mathcal{R}_{v}^{(1)}\mathbf{K}_{1}^{(0)},\mathbf{r})\rangle$
$\displaystyle=\langle
H_{orb}\Psi^{(e(o))}_{tmdc}(\mathbf{k}^{\prime}_{1}-{\mathbf{b}^{\prime}_{1}},\mathbf{r})|\Psi^{(X)}_{gr}(\mathbf{K}_{1}^{(0)}-{\mathbf{b}_{1}^{(0)}},\mathbf{r})\rangle$
$\displaystyle=\left(\langle\Psi^{(X)}_{gr}(\mathbf{K}_{1}^{(0)}-{\mathbf{b}_{1}^{(0)}},\mathbf{r})|H_{orb}|\Psi^{(e(o))}_{tmdc}(\mathbf{k}^{\prime}_{1}-{\mathbf{b}^{\prime}_{1}},\mathbf{r})\rangle\right)^{*}$
$\displaystyle=\left(e^{-i\omega^{X}_{1}}\langle\Psi^{(X)}_{gr}(\mathbf{K}_{1}^{(0)},\mathbf{r})|H_{orb}|\Psi^{(e(o))}_{tmdc}(\mathbf{k}^{\prime}_{1},\mathbf{r})\rangle
e^{i\omega^{(e(o))}_{1}}\right)^{*}=t_{e(o)}^{*}(\mathbf{k}^{\prime}_{1})e^{i\phi^{X}_{1}}e^{i\omega^{X}_{1}}e^{-i\omega^{(e(o))}_{1}}.$
(23)

Here we have made use of i) time reversal symmetry when assuming
$H^{*}_{orb}=H_{orb}$, ii)
$\mathcal{R}_{v}^{(2)}=(\mathcal{R}_{v}^{(2)})^{-1}$, and iii)
$\mathcal{R}_{v}^{(2)}$ is a symmetry of the system therefore
$\mathcal{R}_{v}^{(2)}H_{orb}=H_{orb}\mathcal{R}_{v}^{(2)}$. Moreover, Bloch
wavefunctions at equivalent wavenumbers can differ at most by a phase factor,
i.e.,
$|\Psi^{(X)}_{gr}(\mathbf{K}_{1}^{(0)}-{\mathbf{b}_{1}^{(0)}},\mathbf{r})\rangle=|\Psi^{(X)}_{gr}(\mathbf{K}_{1}^{(0)},\mathbf{r})\rangle
e^{i\omega^{X}_{1}}$ and
$|\Psi^{(e(o))}_{tmdc}(\mathbf{k}^{\prime}_{1}-{\mathbf{b}^{\prime}_{1}},\mathbf{r})\rangle=|\Psi^{(e(o))}_{tmdc}(\mathbf{k}^{\prime}_{1},\mathbf{r})\rangle
e^{i\omega^{(e(o))}_{1}}$. From Eq. (23) one finds that

$t_{e,o}(\mathbf{k}^{\prime}_{1})=t_{e}(\mathbf{k}^{\prime}_{1})t_{o}^{*}(\mathbf{k}^{\prime}_{1})=t_{e,o}^{*}(\mathbf{k}^{\prime}_{1})e^{-i\omega^{(e)}_{1}}e^{i\omega^{(o)}_{1}}.$
(24)

A similar relation can be derived for the matrix elements of $\hat{H}_{soc}$,
namely, for
$\langle\Psi^{(e)}(\mathbf{k}^{\prime}_{1},\mathbf{r})|L_{\pm}|\Psi^{(o)}(\mathbf{k}^{\prime}_{1},\mathbf{r})\rangle$:

$\displaystyle\langle\Psi^{(e)}(\mathbf{k}^{\prime}_{1},\mathbf{r})|L_{\pm}|\Psi^{(o)}(\mathbf{k}^{\prime}_{1},\mathbf{r})\rangle$
$\displaystyle=\langle(L_{\pm}\Psi^{(o)}(\mathbf{k}^{\prime}_{1},\mathbf{r}))^{*}|(\Psi^{(e)}(\mathbf{k}^{\prime}_{1},\mathbf{r}))^{*}\rangle=-\langle
L_{\mp}\Psi^{(o)}(-\mathbf{k}^{\prime}_{1},\mathbf{r})|\Psi^{(e)}(-\mathbf{k}^{\prime}_{1},\mathbf{r})\rangle$
$\displaystyle=-\langle
L_{\mp}\Psi^{(o)}(\mathcal{R}_{v}^{(2)}(\mathbf{k}^{\prime}_{1}-\mathbf{b}^{\prime}_{1}),\mathbf{r})|\Psi^{(e)}(\mathcal{R}_{v}^{(2)}(\mathbf{k}^{\prime}_{1}-\mathbf{b}^{\prime}_{1}),\mathbf{r})\rangle$
$\displaystyle=-\langle
L_{\mp}\mathcal{R}_{v}^{(2)}\Psi^{(o)}(\mathbf{k}^{\prime}_{1}-\mathbf{b}^{\prime}_{1},\mathbf{r})|\mathcal{R}_{v}^{(2)}\Psi^{(e)}(\mathbf{k}^{\prime}_{1}-\mathbf{b}^{\prime}_{1},\mathbf{r})\rangle$
$\displaystyle=-\langle e^{\pm
i\frac{2\pi}{3}}\mathcal{R}_{v}^{(2)}[L_{\pm}\Psi^{(o)}(\mathbf{k}^{\prime}_{1}-\mathbf{b}^{\prime}_{1},\mathbf{r})]|\mathcal{R}_{v}^{(2)}\Psi^{(e)}(\mathbf{k}^{\prime}_{1}-\mathbf{b}^{\prime}_{1},\mathbf{r})\rangle$
$\displaystyle=-e^{\mp
i\frac{2\pi}{3}}\langle[L_{\pm}\Psi^{(o)}(\mathbf{k}^{\prime}_{1}-\mathbf{b}^{\prime}_{1},\mathbf{r})]|\Psi^{(e)}(\mathbf{k}^{\prime}_{1}-\mathbf{b}^{\prime}_{1},\mathbf{r})\rangle$
$\displaystyle=-e^{\mp
i\frac{2\pi}{3}}\left(\langle\Psi^{(e)}(\mathbf{k}^{\prime}_{1}-\mathbf{b}^{\prime}_{1},\mathbf{r})|L_{\pm}|\Psi^{(o)}(\mathbf{k}^{\prime}_{1}-\mathbf{b}^{\prime}_{1},\mathbf{r})\rangle\right)^{*}$
$\displaystyle=-e^{\mp
i\frac{2\pi}{3}}\left(\langle\Psi^{(e)}(\mathbf{k}^{\prime}_{1},\mathbf{r})|L_{\pm}|\Psi^{(o)}(\mathbf{k}^{\prime}_{1},\mathbf{r})\rangle\right)^{*}e^{i\omega^{(e)}_{1}}e^{-i\omega^{(o)}_{1}}.$
(25)

In terms of $\alpha_{e,o}^{(\pm)}(\mathbf{k}^{\prime}_{1})$, introduced after
Eq. (6), Eq. (25) means that
$\alpha_{e,o}^{(\pm)}(\mathbf{k}^{\prime}_{1})=-e^{\mp
i\frac{2\pi}{3}}[\alpha_{e,o}^{(\pm)}(\mathbf{k}^{\prime}_{1})]^{*}e^{i\omega^{(e)}_{1}}e^{-i\omega^{(o)}_{1}}$.
Therefore,

$\displaystyle\alpha_{e,o}^{(x)}({\mathbf{k}}_{1})$
$\displaystyle=\frac{1}{2}\left(\alpha^{(+)}_{e,o}(\mathbf{k}^{\prime}_{1})+\alpha^{(-)}_{e,o}(\mathbf{k}^{\prime}_{1})\right)=-\frac{1}{2}\left(e^{-i\frac{2\pi}{3}}[\alpha^{(+)}_{e,o}(\mathbf{k}^{\prime}_{1})]^{*}+e^{i\frac{2\pi}{3}}[\alpha^{(-)}_{e,o}(\mathbf{k}^{\prime}_{1})]^{*}\right)e^{i\omega^{(e)}_{1}}e^{-i\omega^{(o)}_{1}}$
$\displaystyle=-\frac{1}{2}\left(-\frac{1}{2}\left([\alpha^{(+)}_{e,o}(\mathbf{k}^{\prime}_{1})]^{*}+[\alpha^{(-)}_{e,o}(\mathbf{k}^{\prime}_{1})]^{*}\right)-i\frac{\sqrt{3}}{2}\left([\alpha^{(+)}_{e,o}(\mathbf{k}^{\prime}_{1})]^{*}-[\alpha^{(-)}_{e,o}(\mathbf{k}^{\prime}_{1})]^{*}\right)\right)e^{i\omega^{(e)}_{1}}e^{-i\omega^{(o)}_{1}},$
(26a) $\displaystyle\alpha_{e,o}^{(y)}({\mathbf{k}}_{1})$
$\displaystyle=\frac{1}{2i}\left(\alpha^{(+)}_{e,o}(\mathbf{k}^{\prime}_{1})-\alpha^{(-)}_{e,o}(\mathbf{k}^{\prime}_{1})\right)=-\frac{1}{2i}\left(e^{-i\frac{2\pi}{3}}[\alpha^{(+)}_{e,o}(\mathbf{k}^{\prime}_{1})]^{*}-e^{i\frac{2\pi}{3}}[\alpha^{(-)}_{e,o}(\mathbf{k}^{\prime}_{1})]^{*}\right)e^{i\omega^{(e)}_{1}}e^{-i\omega^{(o)}_{1}}$
$\displaystyle=-\frac{1}{2i}\left(-\frac{1}{2}\left([\alpha^{(+)}_{e,o}(\mathbf{k}^{\prime}_{1})]^{*}-[\alpha^{(-)}_{e,o}(\mathbf{k}^{\prime}_{1})]^{*}\right)-i\frac{\sqrt{3}}{2}\left([\alpha^{(+)}_{e,o}(\mathbf{k}^{\prime}_{1})]^{*}+[\alpha^{(-)}_{e,o}(\mathbf{k}^{\prime}_{1})]^{*}\right)\right)e^{i\omega^{(e)}_{1}}e^{-i\omega^{(o)}_{1}}.$
(26b)

Eqs. (26) can be written as

$\left(\begin{array}[]{c}\alpha_{e,o}^{(x)}(\mathbf{k}^{\prime}_{1})\\\
\alpha_{e,o}^{(y)}(\mathbf{k}^{\prime}_{1})\end{array}\right)=-\frac{1}{2}\left(\begin{array}[]{cc}1&\sqrt{3}\\\
\sqrt{3}&-1\end{array}\right)\left(\begin{array}[]{c}\left(\alpha_{e,o}^{(x)}(\mathbf{k}^{\prime}_{1})\right)^{*}\\\
\left(\alpha_{e,o}^{(y)}(\mathbf{k}^{\prime}_{1})\right)^{*}\end{array}\right)e^{i\omega^{(e)}_{1}}e^{-i\omega^{(o)}_{1}},$
(27)

One now calculate ${\vartheta}_{e,o}(\mathbf{k}^{\prime}_{1})$ for $\theta=0$.
First, combining Eq. (24) and Eq. (27) one arrives at

$t_{e,o}(\mathbf{k}^{\prime}_{1})\left(\begin{array}[]{c}\alpha_{e,o}^{(x)}(\mathbf{k}^{\prime}_{1})\\\
\alpha_{e,o}^{(y)}(\mathbf{k}^{\prime}_{1})\end{array}\right)=-\frac{1}{2}\left(\begin{array}[]{cc}1&\sqrt{3}\\\
\sqrt{3}&-1\end{array}\right)\left(\begin{array}[]{c}\left(\alpha_{e,o}^{(x)}(\mathbf{k}^{\prime}_{1})\right)^{*}\\\
\left(\alpha_{e,o}^{(y)}(\mathbf{k}^{\prime}_{1})\right)^{*}\end{array}\right)t_{e,o}^{*}(\mathbf{k}^{\prime}_{1}).$
(28)

This means that ${\rm
Im}\left[t_{e,o}(\mathbf{k}^{\prime}_{1})\alpha_{e,o}^{(x)}(\mathbf{k}^{\prime}_{1})\right]=-\frac{1}{2}{\rm
Im}\left[[t_{e,o}(\mathbf{k}^{\prime}_{1})\alpha_{e,o}^{(x)}(\mathbf{k}^{\prime}_{1})]^{*}\right]-\frac{\sqrt{3}}{2}{\rm
Im}\left[[t_{e,o}(\mathbf{k}^{\prime}_{1})\alpha_{e,o}^{(y)}(\mathbf{k}^{\prime}_{1})]^{*}\right],$
i.e., ${\rm
Im}\left[t_{e,o}(\mathbf{k}^{\prime}_{1})\alpha_{e,o}^{(x)}(\mathbf{k}^{\prime}_{1})\right]=\sqrt{3}\,{\rm
Im}\left[t_{e,o}(\mathbf{k}^{\prime}_{1})\alpha_{e,o}^{(y)}(\mathbf{k}^{\prime}_{1})\right].$
Regarding Eq. (12), it follows that

${\rm Arctan}\left[\frac{{\rm
Im}[e^{i\eta(\mathbf{k}^{\prime}_{1})}\alpha_{e,o}^{(x)}(\mathbf{k}^{\prime}_{1})]}{{\rm
Im}[e^{i\eta(\mathbf{k}^{\prime}_{1})}\alpha_{e,o}^{(y)}(\mathbf{k}^{\prime}_{1})]}\right]={\rm
Arctan}\left[\frac{{\rm
Im}[t_{e,o}(\mathbf{k}^{\prime}_{1})\alpha_{e,o}^{(x)}(\mathbf{k}^{\prime}_{1})]}{{\rm
Im}[t_{e,o}(\mathbf{k}^{\prime}_{1})\alpha_{e,o}^{(y)}(\mathbf{k}^{\prime}_{1})]}\right]=\frac{\pi}{3},$
(29)

and using Eq. (14) one finds
${\vartheta}_{e,o}(\mathbf{k}^{\prime}_{1})=(m+1)\pi$. Therefore, one finds
from Eq. (22) that for $\theta=0$
$H_{R}^{gr}=(-1)^{m}\frac{\lambda_{R}(\mathbf{k}^{\prime}_{1})}{2}\left(\tau_{z}\sigma_{x}s_{y}-\sigma_{y}s_{x}\right).$

####  I.5.2 $\theta=\pi/6$ twist angle

Graphene has two, inequivalent sets of vertical reflection planes, see Fig. 1.
For $\theta=\pi/6$, the $\mathcal{R}_{d}^{(1,2,3)}$ set aligns with the TMDC’s
$\mathcal{R}_{v}^{(1,2,3)}$ reflection planes, therefore the stack again has
$C_{3v}$ symmetry.

We denote the three wavevectors in the BZ of the TMDC that satisfy the
quasimomentum conservation by $\tilde{\mathbf{k}}^{\prime}_{1,2,3}$, see Fig.
4. As one can see, the $\tilde{\mathbf{k}}^{\prime}_{j}$ vectors now lie on
the $\Gamma$-$M$ lines of the BZ of the TMDC.

Figure 4: The BZ of graphene (purple) and of the TMDC (orange), along with
$\tilde{\mathbf{k}}^{\prime}_{1,2,3}$ vectors and the reflection planes
$\mathcal{R}_{v}^{1,2,3}$ for $\theta=\pi/6$.

Let us consider the matrix element
$\langle\Psi^{(e)}(\tilde{\mathbf{k}}_{1},\mathbf{r})|L_{\pm}|\Psi^{(o)}(\tilde{\mathbf{k}}_{1},\mathbf{r})\rangle$.

$\displaystyle\langle\Psi^{(e)}(\tilde{\mathbf{k}}^{\prime}_{1},\mathbf{r})|L_{\pm}|\Psi^{(o)}(\tilde{\mathbf{k}}^{\prime}_{1},\mathbf{r})\rangle$
$\displaystyle=\langle\Psi^{(e)}(\tilde{\mathbf{k}}^{\prime}_{1},\mathbf{r})|(\mathcal{R}_{v}^{(3)})^{\dagger}\mathcal{R}_{v}^{(3)}L_{\pm}(\mathcal{R}_{v}^{(3)})^{\dagger}\mathcal{R}_{v}^{(3)}|\Psi^{(o)}(\tilde{\mathbf{k}}^{\prime}_{1},\mathbf{r})\rangle$
$\displaystyle=\langle\Psi^{(e)}(\mathcal{R}_{v}^{(3)}\tilde{\mathbf{k}}^{\prime}_{1},\mathbf{r})|\mathcal{R}_{v}^{(3)}L_{\pm}(\mathcal{R}_{v}^{(3)})^{\dagger}|\Psi^{(o)}(\mathcal{R}_{v}^{(3)}\tilde{\mathbf{k}}^{\prime}_{1},\mathbf{r})\rangle$
$\displaystyle=\langle\Psi^{(e)}(\tilde{\mathbf{k}}^{\prime}_{1},\mathbf{r})|\mathcal{R}_{v}^{(3)}L_{\pm}(\mathcal{R}_{v}^{(3)})^{\dagger}|\Psi^{(o)}(\tilde{\mathbf{k}}^{\prime}_{1},\mathbf{r})\rangle$
$\displaystyle=e^{\pm
i\frac{2\pi}{3}}\langle\Psi^{(e)}(\tilde{\mathbf{k}}^{\prime}_{1},\mathbf{r})|L_{\mp}|\Psi^{(o)}(\tilde{\mathbf{k}}^{\prime}_{1},\mathbf{r})\rangle,$
(30)

where we made use of
$\mathcal{R}_{v}^{(3)}L_{\pm}(\mathcal{R}_{v}^{(3)})^{\dagger}=e^{\pm
i\frac{2\pi}{3}}L_{\mp}$ and that
$\mathcal{R}_{v}^{(3)}\tilde{\mathbf{k}}^{\prime}_{1}=\tilde{\mathbf{k}}^{\prime}_{1}$.
In terms of $\alpha_{e,o}^{(\pm)}(\mathbf{k}^{\prime}_{1})$, Eq.(30) can be
rewritten as
$\alpha_{e,o}^{(-)}(\tilde{\mathbf{k}}^{\prime}_{1})=e^{-i\frac{2\pi}{3}}\alpha_{e,o}^{(+)}(\tilde{\mathbf{k}}^{\prime}_{1})$.
Therefore,

$\displaystyle\alpha_{e,o}^{(x)}(\tilde{\mathbf{k}}^{\prime}_{1})$
$\displaystyle=\frac{1}{2}\left(\alpha_{e,o}^{(+)}(\tilde{\mathbf{k}}^{\prime}_{1})+\alpha_{e,o}^{(-)}(\tilde{\mathbf{k}}^{\prime}_{1})\right)=\frac{1}{2}\left(1+e^{-i\frac{2\pi}{3}}\right)\alpha_{e,o}^{(+)}(\tilde{\mathbf{k}}^{\prime}_{1}),$
(31a) $\displaystyle\alpha_{e,o}^{(y)}(\tilde{\mathbf{k}}^{\prime}_{1})$
$\displaystyle=\frac{1}{2i}\left(\alpha_{e,o}^{(+)}(\tilde{\mathbf{k}}^{\prime}_{1})-\alpha_{e,o}^{(-)}(\tilde{\mathbf{k}}^{\prime}_{1})\right)=\frac{1}{2i}\left(1-e^{-i\frac{2\pi}{3}}\right)\alpha_{e,o}^{(+)}(\tilde{\mathbf{k}}^{\prime}_{1}).$
(31b)

It follows that
$\alpha_{e,o}^{(y)}(\tilde{\mathbf{k}}^{\prime}_{1})=\sqrt{3}\alpha_{e,o}^{(x)}(\tilde{\mathbf{k}}^{\prime}_{1})$.
Regarding Eq. (12), one finds ${\rm Arctan}\left[\frac{{\rm
Im}[e^{i\eta(\tilde{\mathbf{k}}^{\prime}_{1})}\alpha_{e,o}^{(x)}(\tilde{\mathbf{k}}^{\prime}_{1})]}{{\rm
Im}[e^{i\eta(\tilde{\mathbf{k}}^{\prime}_{1})}\alpha_{e,o}^{(y)}(\tilde{\mathbf{k}}^{\prime}_{1})]}\right]=\frac{\pi}{6},$
and from Eq. (14)
$\vartheta_{e,o}(\tilde{\mathbf{k}}^{\prime}_{1})=\frac{5\pi}{6}+m^{\prime}\pi=(m^{\prime}+1)\pi-\pi/6$.

####  I.5.3 $\theta=\pi/3$ twist angle

We again denote the three wavevectors in the BZ of the TMDC that satisfy the
quasimomentum conservation by $\tilde{\mathbf{k}}^{\prime}_{1,2,3}$, see Fig.
5. Note, that $\tilde{\mathbf{k}}^{\prime}_{1}$ is related to
$\mathbf{k}^{\prime}_{1}$ in the $\theta=0$ case by
$\tilde{\mathbf{k}}^{\prime}_{1}=\mathcal{R}_{v}^{(3)}\mathbf{k}^{\prime}_{1}$
and the $\mathbf{K}^{(\pi/3)}$ point of graphene is given by
$-\mathbf{K}^{(0)}+\mathbf{b}_{1}^{0}$.

Figure 5: The BZ of graphene (purple) and of the TMDC (orange), along with
$\tilde{\mathbf{k}}^{\prime}_{1,2,3}$ vectors and the reflection planes
$\mathcal{R}_{v}^{1,2,3}$ for $\theta=\pi/3$.

This can be used to calculate the amplitude
$t_{e(o)}(\tilde{\mathbf{k}}^{\prime}_{1})$ and
$\alpha_{e,o}^{(\pm)}(\tilde{\mathbf{k}}^{\prime}_{1})$. Similarly to Eq.
(23), one may write

$\displaystyle t_{e(o)}(\tilde{\mathbf{k}}^{\prime}_{1})e^{-i\phi_{1}^{X}}$
$\displaystyle=\langle\Psi^{(X)}_{gr}({\mathbf{K}}_{1}^{(\pi/3)},\mathbf{r})|H_{orb}|\Psi^{(e(o))}_{tmdc}(\tilde{\mathbf{k}}^{\prime}_{1},\mathbf{r})\rangle$
$\displaystyle=\langle\Psi^{(X)}_{gr}({-\mathbf{K}}_{1}^{(0)}+\mathbf{b}_{1}^{(0)},\mathbf{r})|H_{orb}|\Psi^{(e(o))}_{tmdc}(\mathcal{R}_{v}^{(3)}\mathbf{k}^{\prime}_{1},\mathbf{r})\rangle$
$\displaystyle=\langle\Psi^{(X)}_{gr}({-\mathbf{K}}_{1}^{(0)}+\mathbf{b}_{1}^{(0)},\mathbf{r})|H_{orb}|\Psi^{(e(o))}_{tmdc}(\mathcal{R}_{v}^{(3)}\mathbf{k}^{\prime}_{1},\mathbf{r})\rangle$
$\displaystyle=\langle\Psi^{(X)}_{gr}({-\mathbf{K}}_{1}^{(0)}+\mathbf{b}_{1}^{(0)},\mathbf{r})|H_{orb}|\mathcal{R}_{v}^{(3)}\Psi^{(e(o))}_{tmdc}(\mathbf{k}^{\prime}_{1},\mathbf{r})\rangle$
$\displaystyle=\langle\Psi^{(X)}_{gr}({-\mathbf{K}}_{1}^{(0)}+\mathbf{b}_{1}^{(0)},\mathbf{r})|\mathcal{R}_{v}^{(3)}H_{orb}|\Psi^{(e(o))}_{tmdc}(\mathbf{k}^{\prime}_{1},\mathbf{r})\rangle$
$\displaystyle=\langle\Psi^{(X)}_{gr}(\mathcal{R}_{v}^{(3)}({-\mathbf{K}}_{1}^{(0)}+\mathbf{b}_{1}^{(0)}),\mathbf{r})|H_{orb}|\Psi^{(e(o))}_{tmdc}(\mathbf{k}^{\prime}_{1},\mathbf{r})\rangle$
$\displaystyle=\langle\Psi^{(X)}_{gr}({\mathbf{K}}_{1}^{(0)}-\mathbf{b}_{1}^{(0)}),\mathbf{r})|H_{orb}|\Psi^{(e(o))}_{tmdc}(\mathbf{k}^{\prime}_{1},\mathbf{r})\rangle=t_{e(o)}(\mathbf{k}^{\prime}_{1})e^{-i\phi_{1}^{X}}e^{-i\omega_{1}^{X}}.$
(32)

This meas that
$t_{e,o}(\tilde{\mathbf{k}}^{\prime}_{1})=t_{e}(\tilde{\mathbf{k}}^{\prime}_{1})t_{o}^{*}(\tilde{\mathbf{k}}^{\prime}_{1})=t_{e}(\mathbf{k}^{\prime}_{1})t_{o}^{*}(\mathbf{k}^{\prime}_{1})=t_{e,o}(\mathbf{k}^{\prime}_{1})$.
Regarding the matrix elements of $\hat{H}_{soc}$, one finds

$\displaystyle\langle\Psi^{(e)}(\tilde{\mathbf{k}}^{\prime}_{1},\mathbf{r})|L_{\pm}|\Psi^{(o)}(\tilde{\mathbf{k}}^{\prime}_{1},\mathbf{r})\rangle$
$\displaystyle=\langle\Psi^{(e)}(\mathcal{R}_{v}^{(3)}\mathbf{k}^{\prime}_{1},\mathbf{r})|L_{\pm}|\Psi^{(o)}(\mathcal{R}_{v}^{(3)}\mathbf{k}^{\prime}_{1},\mathbf{r})\rangle$
$\displaystyle=\langle\Psi^{(e)}(\mathbf{k}^{\prime}_{1},\mathbf{r})|[\mathcal{R}_{v}^{(3)}]^{\dagger}L_{\pm}\mathcal{R}_{v}^{(3)}|\Psi^{(o)}(\mathbf{k}^{\prime}_{1},\mathbf{r})\rangle$
$\displaystyle=e^{\pm
i\frac{2\pi}{3}}\langle\Psi^{(e)}(\mathbf{k}^{\prime}_{1},\mathbf{r})|L_{\mp}|\Psi^{(o)}(\mathbf{k}^{\prime}_{1},\mathbf{r})\rangle\rangle$
(33)

Here we have used
$[\mathcal{R}_{v}^{(3)}]^{\dagger}L_{\pm}\mathcal{R}_{v}^{(3)}=e^{\pm
i\frac{2\pi}{3}}L_{\mp}$. Performing analogous calculations as in Eqs. (26),
one finds

$t_{e,o}(\tilde{\mathbf{k}}^{\prime}_{1})\left(\begin{array}[]{c}\alpha_{e,o}^{(x)}(\tilde{\mathbf{k}}^{\prime}_{1})\\\
\alpha_{e,o}^{(y)}(\tilde{\mathbf{k}}^{\prime}_{1})\end{array}\right)=\frac{1}{2}\left(\begin{array}[]{cc}-1&\sqrt{3}\\\
\sqrt{3}&1\end{array}\right)\left(\begin{array}[]{c}\alpha_{e,o}^{(x)}(\mathbf{k}^{\prime}_{1})\\\
\alpha_{e,o}^{(y)}(\mathbf{k}^{\prime}_{1})\end{array}\right)t_{e,o}(\mathbf{k}^{\prime}_{1}).$
(34)

One can now easily calculate
${\vartheta}_{e,o}(\tilde{\mathbf{k}}^{\prime}_{1})$. Using Eq.(34)

${\rm
Im}[t_{e,o}(\tilde{\mathbf{k}}^{\prime}_{1})\alpha_{e,o}^{(x)}(\tilde{\mathbf{k}}^{\prime}_{1})]=-\frac{1}{2}{\rm
Im}[t_{e,o}(\mathbf{k}^{\prime}_{1})\alpha_{e,o}^{(x)}(\mathbf{k}^{\prime}_{1})]+\frac{\sqrt{3}}{2}{\rm
Im}[t_{e,o}(\mathbf{k}^{\prime}_{1})\alpha_{e,o}^{(y)}(\mathbf{k}^{\prime}_{1})].$
(35)

However, in Sec. I.5.1 we found that ${\rm
Im}\left[t_{e,o}(\mathbf{k}^{\prime}_{1})\alpha_{e,o}^{(x)}(\mathbf{k}^{\prime}_{1})\right]=\sqrt{3}\,{\rm
Im}\left[t_{e,o}(\mathbf{k}^{\prime}_{1})\alpha_{e,o}^{(y)}(\mathbf{k}^{\prime}_{1})\right]$,
which means that ${\rm
Im}[t_{e,o}(\tilde{\mathbf{k}}^{\prime}_{1})\alpha_{e,o}^{(x)}(\tilde{\mathbf{k}}^{\prime}_{1})]=0$
and from Eq.(14)
${\vartheta}_{e,o}(\tilde{\mathbf{k}}^{\prime}_{1})=2\pi/3+m^{\prime\prime}\pi=(m^{\prime\prime}+1)\pi-\pi/3$.

###  I.6 Discussion of the results of Sec. I.5

Similar calculations to those in Sec. I.5.1-I.5.3 can be performed for all
interlayer twist angles $\theta_{l}=l\pi/6$, $l=0,1,2,\dots$. We now summarize
the most important findings.

  * •

One finds that $\vartheta_{e,o}(\theta_{l})$ can be expressed as
$\vartheta_{e,o}(\theta_{l})=m(\theta_{l})\pi-\theta_{l}$, where
$m(\theta_{l})$ is an integer.

  * •

The results of Sec. I.5.1-I.5.3 were obtained for an arbitrary pair of $e$ and
$o$ bands. This means, that using Eq. (17) one finds
$\vartheta_{R}(\theta_{l})=n(\theta_{l})\pi-\theta_{l}$, where $n(\theta_{l})$
is an integer.

  * •

According to Eq. (22), the non-zero matrix elements of $H_{R}^{gr}$ are
$\propto\lambda_{R}(\theta)e^{\pm i(\theta+\vartheta_{R})}$. For
$\theta=\theta_{l}$ one finds that these matrix elements are
$\lambda_{R}(\theta_{l})(-1)^{n(\theta_{l})}$. This means that for interlayer
twist angles where the $C_{3v}$ symmetry of the stack is restored, the
Hamiltonian of the induced Rashba SOC simplifies to

$H_{R}^{gr}=(-1)^{[n(\theta_{l})+1]}\frac{{\lambda}_{R}(\theta_{l})}{2}\left(\tau_{z}\sigma_{x}s_{y}-\sigma_{y}s_{x}\right),$
(36)

which is, apart from the sign $(-1)^{[n(\theta_{l})+1]}$, the well-known
result of Refs. kane_quantum-spin-Hall_2005 ; min_intrinsic_2006 . This shows
our results are in agreement with general expectations based on the symmetry
of the system.

  * •

For twist angles $\theta\in(\theta_{l},\theta_{l+1})$, i.e., when the stack
has only $C_{3}$ symmetry, $\vartheta_{R}$ is a continuous function of
$\theta$ (through the wavenumbers $\tilde{\mathbf{k}}^{\prime}_{j}$).
Therefore the matrix elements $\lambda_{R}(\theta)e^{\pm
i(\theta+\vartheta_{R})}$ of $H_{R}^{gr}$ are complex numbers.

  * •

One would expect that ${\lambda}_{R}(\theta)={\lambda}_{R}(\theta+2\pi/3)$. As
shown in Figs. 7(b) and 8(b) below, we find that
${\lambda}_{R}(\theta)={\lambda}_{R}(\theta+\pi/3)$. Note, that Eq.(1) gives
the lowest order non-vanishing contribution. We expect that higher order
contributions in the perturbation series, albeit small in magnitude, would
lead to $2\pi/3$ periodicity. Note, that
$\vartheta_{R}(\theta)\neq\vartheta_{R}(\theta+\pi/3)\,\,{\rm mod}\,\,2\pi$.
This can be seen, e.g., comparing the results of Sec. I.5.1 and Sec. I.5.3.
However, the results of our calculations show that
$\vartheta_{R}(\theta)+\theta$, which gives the phase of the induced Rashba
SOC, does change by multiples of $2\pi$ upon interlayer rotation by
$\theta=\pi/3$.

##  II Details of the numerical calculations for $\vartheta_{R}$,
$\lambda_{R}$ and $\lambda_{vZ}$

In this section we discuss questions related to the choice of material and
other parameters in our numerical calculations and compare some of our results
to corresponding results in Ref. david_induced_2019 .

###  II.1 Interlayer tunneling parameters

In Ref. david_induced_2019  it was argued that for the description of the
tunneling between the monolayer graphene (MLG) sheet and the MoS2 layer it is
sufficient to consider only the $p$ orbitals of the chalcogen atoms and one
can neglect the $d$ orbitals of the metal atoms. Namely, the the $d$ orbitals
are farther away from the carbon atoms of MLG and therefore their overlap with
the carbon $p_{z}$ orbitals should be small. We have extended the calculations
of Ref. david_induced_2019  by taking into account contributions from
tunneling processes involving the $d$ orbitals of the metal atoms as well. The
results shown in the main text have been obtained in this way.

As already mentioned in Sec. I.1, the tunneling matrix
$(T({\tau\mathbf{k}^{\prime}_{j}}))_{b}$ can be written as

$\left(T(\tau\bm{k}^{\prime}_{j})\right)_{b}=t_{b}(\tau\mathbf{k}^{\prime}_{j})e^{i\tau\bm{G}^{\theta}_{j}\cdot\bm{r}_{0}}\left(\begin{array}[]{c}e^{i\tau\phi_{j}^{A}}\\\
e^{i\tau\phi_{j}^{B}}\end{array}\right).$

When only the $p$ orbitals are considered in the interlayer tunneling, then
the interlayer tunneling amplitude
$t_{b}^{p}(\tau\mathbf{K}^{\theta},\tau\mathbf{k}^{\prime}_{j})\equiv
t_{b}^{p}(\tau\mathbf{k}^{\prime}_{j})$ between the $\mathbf{K}^{\theta}$
point of graphene and electronic states at wavenumber
$\mathbf{k}^{\prime}_{j}$ in band $b$ of the TMDC is given by

$t_{b}^{p}(\tau\mathbf{k}^{\prime}_{j})=i\tau[c_{bx}(\tau\mathbf{k}^{\prime}_{j})\cos\theta+c_{by}(\tau\mathbf{k}^{\prime}_{j})\sin\theta]\;t_{\parallel}^{p}+c_{bz}(\tau\mathbf{k}^{\prime}_{j})\;t_{\perp}^{p}.$
(37)

Here $c_{bx}(\tau\mathbf{k}^{\prime}_{j})$,
$c_{by}(\tau\mathbf{k}^{\prime}_{j})$ and
$c_{bz}(\tau\mathbf{k}^{\prime}_{j})$ are the complex amplitudes of the
$p_{x}$, $p_{y}$ and $p_{z}$ orbitals of the chalcogen atoms at the BZ point
$\mathbf{k}^{\prime}_{j}$ in band $b$ of the TMDC. As it was shown in Ref.
david_induced_2019 ,
$t_{b}^{p}(\tau\mathbf{k}^{\prime}_{2,3})=t_{b}^{p}(\tau\mathbf{k}^{\prime}_{1})$.
According to Eq. (37), two independent tunneling coefficients are needed:
$t_{\parallel}^{p}$ and $t_{\perp}^{p}$, and Ref. david_induced_2019  obtained
them from a fitting procedure using previous results Gmitra_TMDC-
proximity_2016  of DFT band structures calculations for MLG/TMDC
heterostructures. The values of these parameters were estimated to be
$t_{\parallel}^{p}\approx t_{\perp}^{p}\approx 100\,$meV.

Considering now the $d$ orbitals of the TMDC metal atoms one can derive a
similar expression for the tunneling amplitude
$t_{b}^{d}(\tau\mathbf{k}^{\prime}_{j})$ between MLG and the $d$ orbitals:

$\displaystyle
t_{b}^{d}(\tau\mathbf{k}^{\prime}_{2,3})=t_{b}^{d}(\tau\mathbf{k}^{\prime}_{1})$
$\displaystyle=$
$\displaystyle\left(c_{b,x^{2}-y^{2}}(\tau{\bm{k}^{\prime}_{1}})\cos(2\theta)+c_{b,xy}(\tau{\bm{k}^{\prime}_{1}})\sin(2\theta)\right)t^{d}_{\square}+$
(38) $\displaystyle
i\tau\left(c_{b,xz}(\tau{\bm{k}^{\prime}_{1}})\cos(\theta)+c_{b,yz}(\tau{\bm{k}^{\prime}_{1}})\sin(\theta)\right)t^{d}_{||}+c_{b,3z^{2}-r^{2}}(\tau{\bm{k}^{\prime}_{1}})t^{d}_{\perp}.$

Here $c_{b,d}(\tau{\bm{k}^{\prime}_{1}})$
($d=\\{x^{2}-y^{2},xy,xz,yz,3z^{2}-r^{2}\\}$) are the complex amplitudes of
the $d$ orbitals of the TMDC metal atoms at the BZ point
$\mathbf{k}^{\prime}_{1}$ in band $b$. The derivation of Eq(38) involves the
same steps as that of Eq. (37) and therefore we do not give further details
here. As one can see, $t_{b}^{d}(\tau\mathbf{k}^{\prime}_{1})$ involves three
more overlap parameters: $t^{d}_{\square}$, $t^{d}_{||}$, $t^{d}_{\perp}$.
Therefore there are altogether five $t^{p,d}$ parameters that describe the
overlap between graphene’s $p_{z}$ orbitals and the $p$, $d$ orbitals of the
TMDC layer. The approach used in Ref. david_induced_2019  to estimate
$t_{\parallel}^{p}$ and $t_{\perp}^{p}$ is not applicable to estimate all five
$t^{p,d}$ parameters.

Therefore we used the Slater–Koster method to re-calculate the tunneling
coefficients both for the chalcogen $p$ orbitals and the $d$ orbitals of the
metal atoms. The steps needed to calculate $t_{\parallel}^{p}$,
$t_{\perp}^{p}$, $t^{d}_{\square}$, $t^{d}_{||}$, and $t^{d}_{\perp}$ in terms
of the set of parameters of the Slater–Koster method are explained in Refs.
david_induced_2019 ; li_twist-angle_2019 . We only note that the
parametrization of the transfer integrals is based on the Harrison’s
modelharrison_elementary_1999  and the numerical values of the necessary
parameters, with the exception of $\eta_{\rm C-Se}$, are listed in
Ref.li_twist-angle_2019 , see Eqs(4), (5) therein. All our calculations for
MLG/MoS2 are based on the parameters given in Ref.li_twist-angle_2019 .

However, for the calculations for MLG/WSe2 we also needed the Harrison’s model
parameter $\eta_{\rm C-Se}$. We have fixed its value in the following way. We
took the values $\lambda_{vZ}$ and $\lambda_{R}$ of the induced valley Zeeman
and Rashba SOC from Ref. Gmitra_TMDC-proximity_2016 , where they were obtained
by fitting band structure calculations of graphene/WSe2 heterostructures at
$\theta=0$ twist angle. Using the methodology of Ref. david_induced_2019  to
calculate $\lambda_{vZ}$ and $\lambda_{R}$, we adjusted the value of
$\eta_{\rm C-Se}$ such that we obtain a reasonable agreement with Ref.
Gmitra_TMDC-proximity_2016 , for $\theta=0$. We found that $\eta_{\rm
C-Se}=0.748$ minimizes the deviation to $\approx\pm 20\%$ between our results
for $\lambda_{vZ}$, $\lambda_{R}$ and those of Ref. Gmitra_TMDC-proximity_2016
.

|  $t_{\parallel}^{p}$ [meV] |  $t_{\perp}^{p}$ [meV] |  $t^{d}_{\square}$ [meV] |  $t^{d}_{||}$ [meV] |  $t^{d}_{\perp}$ [meV]
---|---|---|---|---|---
MoS2 | $-143.6$ | $144.7$ | $9.2$ | $17.6$ | $-14.9$
WSe2 | $-76.2$ | $76.2$ | $2.2$ | $4.1$ | $-3.5$
Table 1: Tunneling coefficients

To summarize, we used tunneling coefficients shown in Table 1. Regarding the
$t^{p}$ coefficients for MoS2, they are indeed not very different from those
estimated in Ref. david_induced_2019  ($\approx 100\,$meV) using a different
method, but the sign of $t_{||}^{p}$ is opposite to the one given in Ref.
david_induced_2019 . This can be traced back to the fact that the TB model of
Ref.fang_ab_2015  uses a different orientation of the TMDC lattice with
respect to Ref. david_induced_2019 . This gave rise to an overlooked
$\mathbf{k}\to-\mathbf{k}$ transformation in reciprocal space, which resulted
in a sign change in the estimation of the tunneling parameter. One can also
see that the values of the $t^{d}$ coefficients are indeed much smaller than
$t_{\parallel}^{p}$ and $t_{\perp}^{p}$. On the other hand, the complex
amplitude factors $c_{b,x^{2}-y^{2}}$, $c_{b,xy}$, $c_{b,xz}$, $c_{b,yz}$, and
$c_{b,3z^{2}-r^{2}}$ turn out to be significantly larger than $c_{bx}$,
$c_{by}$ and $c_{bz}$ in large regions of the BZ. Therefore
$t_{b}^{d}(\tau\mathbf{k}^{\prime}_{1})$ is found to be small, but non-
negligible with respect to $t_{b}^{p}(\tau\mathbf{k}^{\prime}_{1})$.

###  II.2 Further parameters

In addition to $t_{b}(\tau\mathbf{k}^{\prime}_{j})$, there are a couple of
further parameters that enter the numerical calculations. We will briefly
discuss them here.

  * •

The value of the band gap $E_{g}$ of the TMDC.

As briefly discussed in the main text, one can take its value either from
experiments (when available) or from previous theoretical works. Regarding the
latter, since we use the TMDC TB model of Ref. fang_ab_2015  to calculate the
band structure, the spin-orbit coupling Hamiltonian matrix elements and the
interlayer tunneling amplitude, we also used the $E_{g}$ of this model. (We
implemented the TB model in the Kwant codegroth_kwant_2014 ). Although the
model of Ref. fang_ab_2015  itself is based on DFT calculations, the value of
$E_{g}$ is different from what one can extract from the calculations of Ref.
Gmitra_TMDC-proximity_2016  that were performed for MLG/TMDC supercells at
$\vartheta=0$. For example, in the case of MoS2 the model of Ref. fang_ab_2015
has a band gap that is $17\%$ larger than the corresponding $E_{g}$ given in
Ref. Gmitra_TMDC-proximity_2016 . See Table 2 for the $E_{g}$ values used in
this work.

  * •

The position of the Dirac point of graphene within the the band gap of the
TMDC.

We describe the energy of the Dirac point of graphene in the band gap of the
TMDC by a number $f_{G}\in[0,1]$. Its value is a linear function of the
position of the Dirac point in the TMDC band gap. When $f_{G}=0$, the Dirac
point is aligned with the TMDC valence band edge, for $f_{G}=1$ the Dirac
point has the same energy as the TMDC conduction band edge. When available,
the value of $f_{G}$ can also be taken from experimentspierucci_graphene-MoS2
; nakamura_graphene-WSe2 . See Table 2 for the $f_{g}$ values used in this
work.

  * •

The number of bands in the model for the TMDC layer.

The explicit expression to calculate $\lambda_{vZ}$ involve a sum over
contributions of individual bands. Similarly, in order to calculate
$\vartheta_{R}$ and $\lambda_{R}$, one needs to sum over pairs of even (e) and
odd (o) bands. Ref. david_induced_2019  used the approximation that the
conduction band (CB) and the valence band (VB) were taken into account for
$\lambda_{vZ}$, and three pairs of e-o bands for $\lambda_{R}$. As it will be
shown below, this already leads to qualitatively good results in most cases,
see Sec. II.3 below. An exception is the value of $\lambda_{R}$ for
$\theta\approx\pi/6$, where, due to band crossings and near-degeneracies of
certain $e$ bands, more pairs of $e$ and $o$ bands need to be taken into
account. The TB model that we use involves $11$ bands of the TMDC layer.
Unless otherwise indicated, we use the contributions of all $11$ bands to
calculate $\lambda_{vZ}$, and $30$ pairs of $e$ and $o$ bands for the
$\vartheta_{R}$ and $\lambda_{R}$ calculations.

|  $E_{g}$ (DFT) |  $E_{g}$ (exp) |  $f_{G}$ (DFT) |  $f_{G}$ (exp)
---|---|---|---|---
MoS2 |  $1.807\,$eVfang_ab_2015  |  $2.0\,$eVpierucci_graphene-MoS2  |  $0.974\,$ Gmitra_TMDC-proximity_2016  |  $0.55$pierucci_graphene-MoS2
WSe2 |  $1.638\,$eVfang_ab_2015  |  $1.95\,$eVnakamura_graphene-WSe2  |  $0.161$Gmitra_TMDC-proximity_2016  |  $0.426$nakamura_graphene-WSe2
Table 2: Band gap and MLG Dirac point energy position parameters from DFT
calculation and experiments.

  * •

The TB model of the TMDC

We also note that the TB model of the TMDC appears to be less accurate for the
bands above the band gapfang_ab_2015 . This can add further uncertainty to the
results, especially in the case graphene/MoS2, where $f_{G}=0.974$, i.e., the
relative contributions of $e$-$o$ band pairs above $E_{g}$ is larger than the
contribution coming from the valence bands.

These details, in addition to those already discussed in Sec. II.1, also need
to be considered when comparing our results to related previous
worksGmitra_TMDC-proximity_2016 ; li_twist-angle_2019 ; david_induced_2019 ;
Pezo2020twist .

###  II.3 Numerical calculations of the induced SOC in MLG/TMDC twisted
bilayers

We now show our results for the twist angle dependence of $\lambda_{vZ}$ and
$\lambda_{R}$ for graphene/MoS2 and graphene/WSe2. In the former case we also
compare our results to those of Ref. david_induced_2019 . For the numerical
calculations of $\lambda_{R}$ we evaluated Eqs. (15) and (17). To calculate
$\lambda_{vZ}$ we used the result of Ref. david_induced_2019 :

$\lambda_{\rm
VZ}=3\sum_{b}\frac{\left|t_{b}(\bm{k}^{\prime}_{1})\right|^{2}\Delta_{0,b}(\bm{k}^{\prime}_{1})}{E_{b}^{2}(\bm{k}^{\prime}_{1})-\Delta_{0,b}^{2}(\bm{k}^{\prime}_{1})},$
(39)

where $t_{b}(\bm{k}^{\prime}_{1})$ is defined in Eq. (3) and the summation
runs over all bands $b$ of the tight-binding model.

####  II.3.1 Comparison to the results of Ref. david_induced_2019

First, we briefly discuss how the tunneling parameters shown in Table 1 modify
the results for tunneling amplitude $t_{b}$ compared to the corresponding
results in Ref. david_induced_2019 . In Fig. 6 we plot the tunneling strength
$|t_{c}|^{2}$ and $|t_{v}|^{2}$ to the CB and VB of MoS2 as a function of the
twist angle $\theta$ using the parameters from Table 1 (solid lines). We also
show the corresponding results of Ref. david_induced_2019 , where only
$t_{\parallel}^{p}\approx t_{\perp}^{p}$ were non-zero.

Figure 6: Tunneling strength $|t_{c}|^{2}$ to the conduction band (light
green) and to the valence band $|t_{v}|^{2}$ (blue) of monolayer MoS2. The
solid curves are the results of this work, the dashed curves are reproduced
from Ref. david_induced_2019 .

As one can see in Fig. 6, the present calculations tend to yield larger
tunneling strengths. However, the twist angle dependence of $|t_{c}|^{2}$ and
$|t_{v}|^{2}$ remain qualitatively the same.

In Fig. 7(a) we show calculations for $\lambda_{vZ}$ vs $\theta$ for
graphene/MoS2 using the DFT parameters for $E_{g}$ and $f_{G}$, see Table 2.
The dashed curve in Fig. 7(a) again corresponds to the result of Ref.
david_induced_2019 , which uses the contributions of the CB and the VB only.
The thin solid line is obtained using the tunneling amplitudes from Table 1
and only tunneling to the VB and the CB is taken into account. One can see a
pronounced increase of $\lambda_{\rm VZ}$ around $\theta=18^{\circ}$ with
respect to Ref. david_induced_2019 , which is mainly due to the larger
tunneling strengths, as shown in Fig. 6. Nevertheless, the curves remain
qualitatively very similar. The thick solid line shows the result that we
obtain by using the tunneling amplitudes from Table 1 and taking into account
tunneling to all $11$ bands of the TB model. It is again qualitatively similar
to the other two results, the most important change is the two new zero
crossings close to $\theta=30^{\circ}$.

Figure 7: a) $\lambda_{vZ}$ and b) $\lambda_{R}$ for MLG/MoS2 as a function of
interlayer twist angle $\theta$ using the DFT parameters given in Table 2.
Different lines correspond to different tunneling amplitudes and number of
bands in the calculations, see text for details. The red marks at $\theta=0$
indicate the results of Ref. Gmitra_TMDC-proximity_2016 .

In Fig. 7(b) we show similar calculations for $\lambda_{R}$ vs $\theta$. The
larger tunneling strength $t_{b}$ used in this work leads again to an
enhancement of $\lambda_{R}$ (thin solid line) compared to Ref.
david_induced_2019  (dashed line), but the $\theta$ dependence of the two
results otherwise agree, including the deep minimum at $\theta=30^{\circ}$.
However, if all possible $e$-$o$ pairs of bands are included in the
calculation (thick solid line), then this minimum becomes a small dip and the
value of $\lambda_{R}$ is significantly enhanced. As already mentioned, this
is because of band crossings and near-degeneracies of certain $e$ bands which
means that more pairs of $e$ and $o$ bands need to be taken into account than
it was done in Ref. david_induced_2019 .

For comparison, in Fig. 7 we also indicate by red marks the results of the DFT
calculations of Ref. Gmitra_TMDC-proximity_2016  for $\theta=0$. As one can
see, all our curves take values close to these reference values. As we
explained in Sec. II.2, the exact results depend on a number of details,
including potentially fine-tuned Slater–Koster parameters. This leaves us with
much freedom to adjust our parameters and fit previously published results,
but this goes beyond the scope of this paper.

We found a qualitatively similar $\lambda_{vZ}$ and $\lambda_{R}$ vs $\theta$
dependence using the experimental parameters of $E_{g}$ and $f_{G}$ from Table
2, but the maximum values of the induced SOC are significantly smaller.

####  II.3.2 Calculations for MLG/WSe2

Our calculations for the twist angle dependence of $\lambda_{vZ}$ and
$\lambda_{R}$ are shown in Fig. 8. We used the tunneling parameters given in
Tables 1 and the DFT parameters in Table 2. Contributions from all bands or
pairs of bands of the TB model are taken into account. Red marks denote the
results of Ref. Gmitra_TMDC-proximity_2016  for $\theta=0$, they differ
$\approx 20\%$ from our results.

Figure 8: a) $\lambda_{vZ}$ and b) $\lambda_{R}$ for MLG/WSe2 as a function of
interlayer twist angle $\theta$ using the DFT parameters given in Table 2. The
red marks at $\theta=0$ indicate the results of Ref. Gmitra_TMDC-
proximity_2016 .

Importantly, one can see that both $\lambda_{vZ}(\theta)$ and
$\lambda_{R}(\theta)$ are different from the MLG/MoS2 heterostructures in Fig.
7. This is mainly due to the fact that the Dirac point of graphene is closer
to the VB (we used $f_{G}=0.161$ in these calculations versus $f_{G}=0.976$ in
Fig. 7, i.e., here the valence bands give larger contributions. Similarly to
MLG/MoS2, $\lambda_{vZ}$ and $\lambda_{R}$ vs $\theta$ dependence using the
experimental parameters for $E_{g}$ and $f_{G}$ from Table 2 are qualitatively
similar, but the maximum values of the induced SOC are significantly smaller.

_Comparison of MLG/MoS 2 and MLG/WSe2_

The different $\theta$ dependence of the induced SOC for MLG/MoS2 (Fig. 7) and
MLG/WSe2 (Fig. 8) explain the differences in the results for the corresponding
twisted trilayer systems shown in Fig. 3 of the main text. We point out, in
particular, that $\lambda_{R}$ has a single maximum for MLG/MoS2, while it has
two maxima for MLG/WSe2. Together with the $\theta$ dependence of
$\vartheta_{R}$ given in Fig. 1 of the main text, this can explain the finding
that $\lambda_{R}^{(tls)}$ can be enhanced for WSe2/MLG/WSe2, but not for
MoS2/MLG/MoS2.

####  II.3.3 Comparison to the results of Refs. Naimer2021twistangle  and
Pezo2020twist

During the preparation of this manuscript, a publication by Naimer _et.
al._Naimer2021twistangle  has appeared discussing the twist angle dependence
of the proximity-induced SOC in graphene/TMDC bilayers. Their methodology is
based on DFT calculations. The incommensurate graphene/TMDC heterostructures
were approximated by periodic supercells, which required several
approximations, e.g., straining the graphene sheet.

In agreement with our results, the calculations of Ref. Naimer2021twistangle
indicate that the strength of the proximity induced SOC in graphene/TMDC
heterostructures is tunable by interlayer twist. The twist angle dependence of
$\lambda_{vZ}$ for graphene/MoS2 and graphene/WSe2 are qualitatively similar:
in the former case there is a maximum at $\theta\approx 20^{\circ}$, in the
latter case $\lambda_{vZ}$ monotonically decreases as $\theta$ is changed from
$0^{\circ}$ to $30^{\circ}$. They also find that the phase factor
$\vartheta_{R}$ in $H_{R}$ is in general non-zero and depends on the twist
angle $\theta$. Regarding the $\theta$ dependence of $\lambda_{R}$, for
graphene/MoS2 their result seems to be qualitatively similar to ours, while
for graphene/WSe2 the results show differences. On a more quantitative level,
however, our results for $\lambda_{vZ}$ and $\lambda_{R}$ display a much
stronger enhancement as a function of $\theta$ than in Ref.
Naimer2021twistangle . The reasons for this discrepancy remains to be
investigated in future works.

Another relevant recent work is by Pezo _et. al._Pezo2020twist , who performed
supercell based DFT calculations for graphene/MoTe2 bilayers for three
different interlayer twist angles ($\theta=0,\pi/12,\pi/6$). It was noted in
Ref. Pezo2020twist  that the position of graphene’s Dirac point is closer to
the valence band of MoTe2. Therefore one would expect that the $\theta$
dependence of the induced SOC should be similar to the one in graphene/WSe2,
meaning, e.g., that magnitude of both $\lambda_{vZ}$ and $\lambda_{R}$ should
decrease for $\theta=\pi/6$ compared to their $\theta=0$ value. This seems to
be in agreement with results of Ref. Pezo2020twist .

Finally, we note that both Refs. Naimer2021twistangle  and Pezo2020twist
emphasize the effect of the strain that is always present in their supercell
based DFT calculations.

##  III Supermoiré effects

Let us first consider the moiré effects in a graphene/TMDC system. One can
define the primitive moiré wavelength $\lambda_{m}$ for a hexagonal
heterostructureYankowitz2012  consisting of graphene and a substrate as

$\lambda_{m}=\frac{q}{\sqrt{1+q^{2}-2q\cos\theta}}\,a_{gr},$ (40)

where $q=a_{s}/a_{gr}$ is the ratio of the lattice constants of the substrate
($a_{s}$) and of graphene ($a_{gr}$) and $\theta$ is the interlayer twist
angle. In the case of graphene and e.g., monolayer MoS2 this ratio is
$q\approx 1.283$. The maximum of $\lambda_{m}$ can be found for $\theta=0$
where $\lambda_{m}\approx 1.1$ nm (see Fig.9). As a comparison, the moiré
length scale is $\lambda_{m}=13.9$ nm for a graphene/hBN bilayer at
$\theta=0$.

The moiré potential due to this periodic perturbation leads to, e.g., gap
openings at the moiré Brillouin zone boundariesYankowitz2012  at energies
$E_{M}=\pm\frac{\hbar v_{F}G_{M}}{2}$, where
$G_{M}=\frac{4\pi}{\sqrt{3}\lambda_{m}}$. Using $v_{F}=10^{6}$m/s, for
graphene/MoS2 this energy scale is $E_{M}\approx\pm 2.2$eV for $\theta=0$,
which is far from the Dirac point. Indeed, in the ARPES experimental results
of Ref. pierucci_graphene-MoS2  on graphene/MoS2, superlattice effects were
observed at binding energies $E_{b}=-3.55$eV, from which the authors concluded
that the interlayer twist angle was $\theta\approx 6^{\circ}$ in their sample.
In short, because of the relatively large difference between the graphene and
TMDC lattice constants, no moiré effects are to be expected close to the Dirac
point of graphene for any twist angle $\theta$ in a graphene/TMDC system.

The situation might be different in TMDC/graphene/TMDC trilayers. Namely,
studying hBN/graphene/hBN trilayers, Refs. Peeters-supermoire, ; Koshino-
supermoire,  found that band gaps can appear in the spectrum of these
heterostructures at energies closer to the Dirac point of graphene than what
was found in hBN/graphene bilayers. This can be understood as a consequence of
an interference of the moiré patterns coming from the two hBN layers. Ref.
Peeters-supermoire  found that if one of the interlayer twist angles, e.g.,
$\theta^{(b)}$ is kept fixed at $\theta^{(b)}=0$ and $\theta^{(t)}$ is
changed, then the four longest wavelength component $\lambda_{sm}^{(i)}$ of
such a _supermoiré_ potential in graphene are given by

$\displaystyle\lambda_{sm}^{(1)}$
$\displaystyle=\frac{q}{\sqrt{2-2\cos\theta^{(t)}}}\,a_{gr},$ (41a)
$\displaystyle\lambda_{sm}^{(2)}$
$\displaystyle=\frac{q}{\sqrt{(2-\delta)(1-\cos\theta^{(t)})+\delta^{2}-\sqrt{3}\delta\sin\theta^{(t)}}}\,a_{gr},$
(41b) $\displaystyle\lambda_{sm}^{(3)}$
$\displaystyle=\frac{q}{\sqrt{2-2\cos\theta^{(t)}+3\delta^{2}-2\sqrt{3}\delta\sin\theta^{(t)}}}\,a_{gr},$
(41c) $\displaystyle\lambda_{sm}^{(4)}$
$\displaystyle=\frac{q}{\sqrt{(2-\delta)(1-\cos\theta^{(t)})+\delta^{2}+\sqrt{3}\delta\sin\theta^{(t)}}}\,a_{gr},$
(41d)

where we used the notation $q-1=\delta$. Eqs. (41) are valid also for
TMDC/graphene/TMDC trilayers because the constituent layers have hexagonal
lattices. We plot $\lambda_{m}$ as well as $\lambda_{sm}^{(i)}$ as a function
of $\theta^{(t)}$ using the lattice constants of graphene and MoS2 in Fig. 9.

Figure 9: Moiré wavelength $\lambda_{m}$ calculated from Eq. (40) and
supermoiré wavelengths calculated from Eqs. (41) as a function of the
interlayer twist angle.

As one can see, $\lambda_{sm}^{(2)}$, $\lambda_{sm}^{(3)}$,
$\lambda_{sm}^{(4)}$ are always quite short. However, $\lambda_{sm}^{(1)}$
does increase when the two TMDC layers are nearly aligned. For
$\theta^{(t)}=1^{\circ}$ one finds $\lambda_{sm}^{(1)}=18$nm corresponding to
$E_{M}=0.13$eV, which is comparable to the ones found in hBN/Graphene/hBN.
Based on these results one may expect an interesting interplay of induced SOC
and supermoiré effects in TMDC/graphene/TMDC for
$|\theta^{(b)}-\theta^{(t)}|\lesssim 1^{\circ}$.

##  IV Magnetotransport calculations for graphene $npn$ junctions

In this section we briefly explain the theoretical model that was used in the
magnetotransport calculations through a $npn$ graphene junction with proximity
induced SOC shown in Fig. 4 of the main text. We assume that the graphene
flake is perfectly ballistic and due to external gates it is doped such that
it hosts an $npn$ junction, see Fig. 10. We numerically calculate the
transmission through the junction as a function of perpendicular magnetic
field magnetic field $B_{z}$ and Fermi energy $E_{F}$ using the in-house code
EQuUsequus-code .

Figure 10: Schematics of a $npn$ junction in graphene. The Dirac point of
graphene is shifted in energy due to the doping potential $U(x)$ in the $n$
and $p$ doped regions. The length of the $n$ ($p$) doped region is $L_{n}$(
$L_{p}$). $E_{F}$ denotes the Fermi energy.

We use the nearest-neighbor tight-binding (TB) model of graphene:

$\hat{H}_{orb}=\sum_{i,s}U_{i}\left(a_{i,s}^{\dagger}a_{i,s}+b_{i,s}^{\dagger}b_{i,s}\right)-\gamma\sum_{\langle
ij\rangle}e^{\phi_{ij}}a_{i,s}^{\dagger}b_{j,s}+h.c.$ (42)

where $a_{i,s}^{\dagger}$ ($a_{i,s}$) and $b_{i,s}^{\dagger}$ ($b_{i,s}$) are
creation (annihilation) operators for electrons of spin $s$ on the $A$ and $B$
sublattice site, respectively. $\gamma=2.97\,$eV is the hopping amplitude
between the nearest-neighbor atomic sites $\langle ij\rangle$ in the graphene
lattice, $\phi_{ij}$ encodes the effect of perpendicular magnetic field
through Peierls substitution, and $U_{i}$ is the on-site energy on the atomic
site $i$. Following Ref. CPR-graphene , we used

$U(x)=U_{n}+\frac{U_{p}-U_{n}}{2}\left(\tanh\left[\frac{x-L_{n}}{l_{1}}\right]-\tanh\left[\frac{x-L_{n}-L_{p}}{l_{2}}\right]\right).$
(43)

Here $L_{n}=50$nm gives the length of the $n$-doped regions and we used
$U_{n}=-240\,$meV. The $p$-doped middle region was $L_{p}=150\,$nm long and
the doping was $U_{p}=40\,$meV. The parameters $l_{1}=l_{2}=25\,$nm set the
smoothness of the transition between the $n$ and $p$-doped regions.

In the TB formalism the Bychkov-Rashba SOC can be written askane_quantum-spin-
Hall_2005 ; frank_protected_2018 :

$H_{R}=i\,\frac{\lambda_{R}(\theta)}{3}e^{i\frac{s_{z}}{2}{\vartheta_{R}}}\sum\limits_{\left\langle
i,j\right\rangle,s,s^{\prime}}\left[{a}^{\dagger}_{i,s}\left(\bm{s}\times\mathbf{\widehat{d}}_{\left\langle
i,j\right\rangle}\right)_{z}{b}_{j,s^{\prime}}-h.c.\right]e^{-i\frac{s_{z}}{2}{\vartheta_{R}}}\;.$
(44)

where $\bm{s}=(s_{x},s_{y},s_{z})$ are the Pauli matrices representing the
electron spin operator. Moreover, $\mathbf{\widehat{d}}_{\left\langle
i,j\right\rangle}=\mathbf{d}_{\left\langle i,j\right\rangle}/d$ are unit
vectors, where $\mathbf{d}_{\left\langle i,j\right\rangle}$ points from atom
$j$ to its nearest neighbors $i$ and $d=|\mathbf{d}_{\left\langle
i,j\right\rangle}|$. The corresponding continuum SOC Hamiltonian, that can be
obtained by Fourier transforming Eq(44) and expanding it at the $\pm K$ points
of the Brillouin zone, reads
$H_{R}=\frac{\lambda_{R}(\theta)}{2}e^{i\frac{s_{z}}{2}{\vartheta_{R}}}\left(\tau_{z}\sigma_{x}s_{y}-\sigma_{y}s_{x}\right)e^{-i\frac{s_{z}}{2}{\vartheta_{R}}}.$
As long as only a perpendicular magnetic field is applied, the phase
$\vartheta_{R}$ will not affect the results, because the corresponding terms
in Eq. (44) can be removed by a unitary transformation.

The valley-Zeeman SOC can be written asfrank_protected_2018 :

$H_{vZ}=\frac{i}{3\sqrt{3}}\,\sum\limits_{\left\langle\langle
i,j\right\rangle\rangle
s,s^{\prime}}\left[\lambda_{vZ}^{(A)}\nu_{ij}[s_{z}]_{s,s^{\prime}}{a}^{\dagger}_{i,s}{a}_{j,s^{\prime}}+\lambda_{vZ}^{(B)}\nu_{ij}[s_{z}]_{s,s^{\prime}}{b}^{\dagger}_{i,s}{b}_{j,s^{\prime}}\right]\;.$
(45)

It couples same spins and depends on clockwise ($\nu_{ij}=1$) or
counterclockwise ($\nu_{ij}=-1$) paths along a hexagonal ring from site $j$ to
$i$. In the case of valley-Zeeman SOC
$\lambda_{vZ}^{(A)}=-\lambda_{vZ}^{(B)}$. The term Eq.(45) can also describe
the intrinsic SOC if the SOC strength $\lambda_{I}$ is the same on the two
sublattices, i.e., $\lambda_{I}^{(A)}=\lambda_{I}^{(B)}$. The corresponding
continuum SOC Hamiltonian reads $H_{vZ}=\tau_{z}s_{z}\lambda_{vZ}$.

We used periodic boundary conditions perpendicular to the transport direction
($\hat{x}$ in Fig.10). Thus, the transverse momentum $q_{n}$ is a good quantum
number and the total transmission can be calculated as a sum over all $q_{n}$:

$T(E_{F},B_{z})=\sum\limits_{n}t(q_{n},E_{F},B_{z})\;,$ (46)

where $t(q_{n},E_{F},B_{z})$ is the transmission coefficient for mode $n$.
When there are many transverse momenta the exact form of the boundary
conditions is not important and therefore we used the infinite mass boundary
conditiongraphene-boundary-cond  to obtain $q_{n}$:
$q_{n}=\left(n+\frac{1}{2}\right)\frac{\pi}{W}\;,$ where $n=0,1,2,\dots$ and
$W$ is the width of the junction.

##  V Electron scattering at a planar junction in a graphene-WSe2 system

We will discuss a further setup to investigate twist-angle dependent transport
in a graphene-TMDC system. It is motivated by recent electron optics
experiment of Ref.barnard_absorptive_2017 ; bhandari_imaging_2018  where
collimated electron beams were created in monolayer graphene. By making use of
such collimated electron beams one can test the twist-angle dependence of the
induced SOC through spin-dependent scattering.

To demonstrate this, we consider the TMDC/MLG/TMDC stack shown in Fig. 11. It
consists of a MLG/WSe2 bilayer of a fixed twist angle $\theta^{(b)}$ (purple
region). Part of this bilayer is covered by a top WSe2 layer (orange region)
of twist angle $\theta^{(t)}$. Thus the induced SOC in graphene changes at the
bilayer-trilayer boundary and in the trilayer region it is assumed to be
tunable by $\theta^{(t)}$. The left edge of the top layer defines a junction
along the $y$ axis, where ballistic electrons scatter into various forward-
propagating and back-propagating modes, as indicated in Fig. 11. We calculate
below the reflection coefficient of this junction for a collimated beam of
electrons.

Figure 11: Ballistic electron scattering in MLG/WSe2-WSe2/MLG/WSe2 junction.
The top layer WSe2 is present only on the right hand side (orange). The
electron source produces a collimated beam of $20$ meV electrons (green),
which scatters into different modes at the junction. Transmitted electrons
(red) get absorbed at the far side of the system, while reflected electrons
(blue) are absorbed by the drain electrode shown in blue. The reflection
coefficient of the junction depends on the rotation angle $\theta^{(t)}$ of
the top layer on the right-hand side

The general form of the Hamiltonian is
$H^{gr}_{eff}=H_{orb}^{gr}+H_{R}^{gr}+H_{vZ}^{gr}$, see Sec. I.4 for the
definition of each term. On the left-hand side of the junction $H_{R}^{gr}$ is
given in Eq. (22), with $\theta=\theta^{(b)}$ and
$\lambda_{vZ}=\lambda_{vZ}^{(b)}(\theta^{(b)})$ in $H_{vZ}^{gr}$. On the right
hand side, as shown in the main text, the Rashba Hamiltonian is
$H_{R}^{(tls)}=\frac{\lambda^{(tls)}_{R}}{2}e^{i\frac{s_{z}}{2}{\vartheta^{(tls)}}}\left(\tau_{z}\sigma_{x}s_{y}-\sigma_{y}s_{x}\right)e^{-i\frac{s_{z}}{2}{\vartheta^{(tls)}}},$
while
$\lambda_{vZ}=\lambda_{vZ}^{(b)}(\theta^{(b)})+\lambda_{vZ}^{(t)}(\theta^{(t)})$
in $H_{vZ}^{gr}$. One can easily show that the reflection coefficient is the
same for electrons in both valleys, for concreteness, we perform the
calculation for the $K$ valley. We express the electronic states as plane
waves multiplied by spinors in the basis
$\left\\{|A\uparrow\rangle,|B\uparrow\rangle,|A\downarrow\rangle,|B\downarrow\rangle\right\\}$
and we assume that their energy is $E=20$ meV. Due to translational symmetry
along $\hat{y}$ axis, the $\delta k_{y}$ component of the wavevector of the
electronic states is conserved. The SOC leads to the splitting of the bands,
therefore (at the $E_{F}$ given above) for each direction of propagation one
can find two modes for the incoming electrons with slightly different wave
vectors and opposite pointing spins. We assume a net zero spin polarization
with an equal mixture of the two incoming modes, and for each incoming mode,
one can find two reflected and two refracted solutions. In total this gives
four back-scattered and four forward-scattered modes within one valley.
Depending on the exact parameters and the angle of incidence, some of these
modes can be decaying waves carrying no current. The particle current operator
can be calculated as $\mathcal{J}_{x}=\frac{i}{\hbar}[\mathcal{H},x]=\tau
v_{F}s_{0}\otimes\sigma_{x},$ and in a similar fashion,
$\mathcal{J}_{y}=v_{F}s_{0}\otimes\sigma_{y}.$ The ratio of the current
carried by all the back-scattered waves and that of all the incoming modes
give us the reflection coefficient:

$R=-\sum_{i=1}^{4}\langle\mathcal{J}_{x,i}^{\text{(back)}}\rangle\left/\sum_{i=1}^{2}\langle\mathcal{J}_{x,i}^{\text{(in)}}\rangle\right..$
(47)

In an experimental scenario one can never create a beam of electrons with an
exact direction of propagation. Nevertheless, as shown in
Refs.barnard_absorptive_2017 ; bhandari_imaging_2018  there are experimental
techniques to collimate electrons in graphene from a source. In Ref
barnard_absorptive_2017  absorptive pinhole collimators were devised to
collimate ballistic electrons into a triangular shaped angular distribution
with a half width at half maximum (HWHM) of just below $10^{\circ}$; while in
Ref bhandari_imaging_2018  a cosine-shaped distribution was achieved with a
similar HWHM of $9^{\circ}$ using zigzag contacts to absorb stray electrons.
Therefore we consider the effective reflection index $\overline{R}$ that is
obtained by averaging $R$ with a certain normalized distribution $d$ centered
around an angle of incidence $\alpha$:

$\overline{R}(\alpha)=\int
d(\alpha^{\prime}-\alpha)R(\alpha^{\prime})\text{d}\alpha^{\prime}.$ (48)

In Fig 12, $\overline{R}$ is plotted as function of the rotation angle
$\theta^{(t)}$ of the top TMDC layer for fixed twist angle
$\theta^{(b)}=\pi/4$ of the bottom layer. As shown in Fig.2 of the main text,
by changing $\theta^{(t)}$ all three SOC parameters $\lambda_{R}^{(tls)}$,
$\vartheta^{(tls)}$ and $\lambda_{vZ}^{(tls)}$ discussed for the trilayer case
are changing and they all affect $\overline{R}$. As one can see, the main
features in $\overline{R}$ do not seem to depend much on the specific angular
distribution of the incoming plane waves or on the exact value of the central
angle of incidence.

Figure 12: The effective reflection $\overline{R}$ as a function of
$\theta^{(t)}$. Angular distribution of incoming electrons taken from Ref
barnard_absorptive_2017  (blue) and from Ref bhandari_imaging_2018  (red)
centered at $\alpha=70^{\circ}$. The green curve shows the result for
$\alpha=75^{\circ}$ using the angular distribution from Ref
bhandari_imaging_2018 .

Experimentally, in order to check the predictions in Fig. 12, a
straightforward approach would be to prepare a few different samples with
different $\theta^{(t)}$. Alternatively, a setup similar to Ref. ribeiro-
palau_twistable_2018  could also be feasible, whereby the top layer is rotated
in situ.

## References

  * [1] Alessandro David, Péter Rakyta, Andor Kormányos, and Guido Burkard.  Induced spin-orbit coupling in twisted graphene–transition metal dichalcogenide heterobilayers: Twistronics meets spintronics.  Phys. Rev. B, 100(8):085412, August 2019.
  * [2] Andor Kormányos, Viktor Zólyomi, Neil D. Drummond, and Guido Burkard.  Spin-Orbit Coupling, Quantum Dots, and Qubits in Monolayer Transition Metal Dichalcogenides.  Phys. Rev. X, 4(1):011034, March 2014.
  * [3] M. S Dresselhaus, G Dresselhaus, and A Jorio.  Group theory: application to the physics of condensed matter.  Springer-Verlag, Berlin, 2010.  OCLC: 692760083.
  * [4] Yang Li and Mikito Koshino.  Twist-angle dependence of the proximity spin-orbit coupling in graphene on transition-metal dichalcogenides.  Phys. Rev. B, 99(7):075438, February 2019.  Publisher: American Physical Society.
  * [5] C. L. Kane and E. J. Mele.  Quantum Spin Hall Effect in Graphene.  Phys. Rev. Lett., 95(22):226801, November 2005.
  * [6] Hongki Min, J. E. Hill, N. A. Sinitsyn, B. R. Sahu, Leonard Kleinman, and A. H. MacDonald.  Intrinsic and Rashba spin-orbit interactions in graphene sheets.  Phys. Rev. B, 74(16):165310, October 2006.
  * [7] Martin Gmitra, Denis Kochan, Petra Högl, and Jaroslav Fabian.  Trivial and inverted Dirac bands and the emergence of quantum spin Hall states in graphene on transition-metal dichalcogenides.  Phys. Rev. B, 93(15):155104, April 2016.  Publisher: American Physical Society.
  * [8] Walter Ashley Harrison.  Elementary Electronic Structure.  World Scientific, 1999.  Google-Books-ID: eAfi3ajSTcUC.
  * [9] Shiang Fang, Rodrick Kuate Defo, Sharmila N. Shirodkar, Simon Lieu, Georgios A. Tritsaris, and Efthimios Kaxiras.  Ab initio tight-binding Hamiltonian for transition metal dichalcogenides.  Phys. Rev. B, 92(20):205108, November 2015.  Publisher: American Physical Society.
  * [10] Christoph W. Groth, Michael Wimmer, Anton R. Akhmerov, and Xavier Waintal.  Kwant: a software package for quantum transport.  New J. Phys., 16(6):063065, June 2014.  Publisher: IOP Publishing.
  * [11] Debora Pierucci, Hugo Henck, Jose Avila, Adrian Balan, Carl H. Naylor, Gilles Patriarche, Yannick J. Dappe, Mathieu G. Silly, Fausto Sirotti, A. T. Charlie Johnson, Maria C. Asensio, and Abdelkarim Ouerghi.  Band Alignment and Minigaps in Monolayer MoS ${}_{\textrm{2}}$ -Graphene van der Waals Heterostructures.  Nano Letters, 16(7):4054–4061, July 2016.
  * [12] H. Nakamura, A. Mohammed, Ph. Rosenzweig, K. Matsuda, K. Nowakowski, K. Küster, P. Wochner, S. Ibrahimkutty, U. Wedig, H. Hussain, J. Rawle, C. Nicklin, B. Stuhlhofer, G. Cristiani, G. Logvenov, H. Takagi, and U. Starke.  Spin splitting and strain in epitaxial monolayer WSe2 on graphene.  Phys. Rev. B, 101:165103, Apr 2020.
  * [13] Armando Pezo, Zeila Zanolli, Nils Wittemeier, Adalberto Fazzio, Stephan Roche, and Jose H. Garcia.  A twist for tuning the spin-orbit coupling in graphene/transition metal dichalcogenide heterobilayers.  arXiv:2011.06714, 2020.
  * [14] Thomas Naimer, Klaus Zollner, Martin Gmitra, and Jaroslav Fabian.  Twist-angle dependent proximity induced spin-orbit coupling in graphene/transition metal dichalcogenide heterostructures.  Phys. Rev. B, 104:195156, Nov 2021.
  * [15] Matthew Yankowitz, Jiamin Xue, Daniel Cormode, Javier D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, Pablo Jarillo-Herrero, Philippe Jacquod, and Brian J. LeRoy.  Emergence of superlattice dirac points in graphene on hexagonal boron nitride.  Nature Physics, 8(5):382–386, May 2012.
  * [16] Miša Andjelković, Slaviša P. Milovanović, Lucian Covaci, and François M. Peeters.  Double moiré with a twist: Supermoiré in encapsulated graphene.  Nano Letters, 20(2):979–988, 2020.  PMID: 31961161.
  * [17] Hiroki Oka and Mikito Koshino.  Fractal energy gaps and topological invariants in hbn/graphene/hbn double moiré systems.  Phys. Rev. B, 104:035306, Jul 2021.
  * [18] Eötvös Quantum Utilities.
  * [19] G. Nanda, J. L. Aguilera-Servin, P. Rakyta, A. Kormányos, R. Kleiner, D. Koelle, K. Watanabe, T. Taniguchi, L. M. K. Vandersypen, and S. Goswami.  Current-phase relation of ballistic graphene josephson junctions.  Nano Letters, 17(6):3396–3401, 2017.
  * [20] Tobias Frank, Petra Högl, Martin Gmitra, Denis Kochan, and Jaroslav Fabian.  Protected Pseudohelical Edge States in ${Z}_{2}$-Trivial Proximitized Graphene.  Phys. Rev. Lett., 120(15):156402, April 2018.
  * [21] A. R. Akhmerov and C. W. J. Beenakker.  Boundary conditions for dirac fermions on a terminated honeycomb lattice.  Phys. Rev. B, 77:085423, Feb 2008.
  * [22] Arthur W. Barnard, Alex Hughes, Aaron L. Sharpe, Kenji Watanabe, Takashi Taniguchi, and David Goldhaber-Gordon.  Absorptive pinhole collimators for ballistic Dirac fermions in graphene.  Nature Communications, 8(1):15418, May 2017.  Number: 1 Publisher: Nature Publishing Group.
  * [23] S. Bhandari, G. H. Lee, K. Watanabe, T. Taniguchi, P. Kim, and R. M. Westervelt.  Imaging electron flow from collimating contacts in graphene.  2D Mater., 5(2):021003, March 2018.  Publisher: IOP Publishing.
  * [24] Rebeca Ribeiro-Palau, Changjian Zhang, Kenji Watanabe, Takashi Taniguchi, James Hone, and Cory R. Dean.  Twistable electronics with dynamically rotatable heterostructures.  Science, 361(6403):690–693, August 2018.