# Twist-angle dependent proximity induced spin-orbit coupling
in graphene/transition-metal dichalcogenide heterostructures

Thomas Naimer<EMAIL_ADDRESS>Institute for
Theoretical Physics, University of Regensburg, 93040 Regensburg, Germany
Klaus Zollner  Institute for Theoretical Physics, University of Regensburg,
93040 Regensburg, Germany  Martin Gmitra  Institute of Physics, Pavol Jozef
Šafárik University in Košice, 04001 Košice, Slovakia  Jaroslav Fabian
Institute for Theoretical Physics, University of Regensburg, 93040 Regensburg,
Germany

###### Abstract

We investigate the proximity-induced spin-orbit coupling in heterostructures
of twisted graphene and monolayers of transition-metal dichalcogenides (TMDCs)
MoS2, WS2, MoSe2, and WSe2 from first principles. We identify strain, which is
necessary to define commensurate supercells, as the key factor affecting the
band offsets and thus magnitudes of the proximity couplings. We establish that
for biaxially strained graphene the band offsets between the Dirac point and
conduction (valence) TMDC bands vary linearly with strain, regardless of the
twist angle. This relation allows to identify the apparent zero-strain band
offsets and find a compensating transverse electric field correcting for the
strain. The resulting corrected band structure is then fitted around the Dirac
point to an established spin-orbit Hamiltonian. This procedure yields the
dominant, valley-Zeeman and Rashba spin-orbit couplings. The magnitudes of
these couplings do not vary much with the twist angle, although the valley-
Zeeman coupling vanishes for 30$\degree$ and Mo-based heterostructures exhibit
a maximum of the coupling at around 20$\degree$. The maximum for W-based
stacks is at 0$\degree$. The Rashba coupling is in general weaker than the
valley-Zeeman coupling, except at angles close to 30$\degree$. We also
identify the Rashba phase angle which measures the deviation of the in-plane
spin texture from tangential, and find that this angle is very sensitive to
the applied transverse electric field. We further discuss the reliability of
the supercell approach with respect to atomic relaxation (rippling of
graphene), relative lateral shifts of the atomic layers, and transverse
electric field.

spintronics, graphene, TMDC, heterostructures, proximity spin-orbit coupling

##  I Introduction

Graphene’s weak (tens of $\mu$eV) spin-orbit coupling (SOC) Gmitra  _et al._
(2009); Sichau  _et al._ (2019) enables long spin diffusion paths but
precludes it from forming a platform for spin manipulation (Han  _et al._ ,
2014). Fortunately, proximity effects enable strong SOC—on the meV scale—in
graphene (Sierra  _et al._ , 2021), making it a suitable candidate for
spintronics devices Žutić  _et al._ (2004); Avsar  _et al._ (2020) and for
observing topological states (Kane and Mele, 2005a, b; Qiao  _et al._ , 2010;
Ren  _et al._ , 2016; Frank  _et al._ , 2018; Högl  _et al._ , 2020).

Unlike pristine graphene, whose spin-orbit coupling is of the Kane-Mele type
Kane and Mele (2005a); Sichau  _et al._ (2019); Gmitra  _et al._ (2009),
graphene proximitized by transition-metal dichalcogenides (TMDCs) exhibits
valley-Zeeman SOC which acts as an effective magnetic field, opposite at $K$
and $K^{\prime}$, but otherwise momentum independent close to these points
(Gmitra and Fabian, 2015; Wang  _et al._ , 2015a; Gmitra  _et al._ , 2016;
Frank  _et al._ , 2018); similarly for graphene on topological insulators Song
_et al._ (2018); Zollner and Fabian (2019, 2021a); Khokhriakov  _et al._
(2018). There is by now plenty of experimental evidence for the large
proximity SOC in graphene (Avsar  _et al._ , 2014; Wang  _et al._ , 2015a;
Yang  _et al._ , 2016; Wang  _et al._ , 2016; Völkl  _et al._ , 2017; Zihlmann
_et al._ , 2018; Khokhriakov  _et al._ , 2018; Garcia  _et al._ , 2018; Safeer
_et al._ , 2019; Herling  _et al._ , 2020; Khoo  _et al._ , 2017; Wang  _et
al._ , 2015a; Omar and van Wees, 2018, 2017; Fülöp  _et al._ , 2021).
Ramifications of the proximity effect include a spin-orbit valve (Gmitra and
Fabian, 2017; Island  _et al._ , 2019; Tiwari  _et al._ , 2021; Amann  _et
al._ , 2021; Khoo  _et al._ , 2017) and giant spin-relaxation anisotropies
(Cummings  _et al._ , 2017; Ghiasi  _et al._ , 2017; Benítez  _et al._ , 2018;
Zihlmann  _et al._ , 2018).

The discovery of superconductivity in magic-angle twisted bilayer graphene
(Cao  _et al._ , 2018; Arora  _et al._ , 2020; Stepanov  _et al._ , 2020;
Balents  _et al._ , 2020) gave the impetus for twistronics Carr  _et al._
(2017); Hennighausen and Kar (2021); Ribeiro-Palau  _et al._ (2018); Carr  _et
al._ (2020), recognizing the twist angle as a critical new control parameter
for tailoring electronic properties of van der Waals heterostructures Zollner
and Fabian (2021b); Yan  _et al._ (2021); David  _et al._ (2019); Li and
Koshino (2019); Pezo  _et al._ (2021). It is then natural to ask how does the
twist angle influence proximity SOC in graphene. Two recent pioneering studies
(Li and Koshino, 2019; David  _et al._ , 2019) based on tight-binding modeling
(Koshino, 2015) suggest that the twist angle between the TMDC and graphene has
a major influence on both the magnitude and the type of proximity induced SOC.
In those works, the interlayer interaction in incommensurate systems is
described by Umklapp processes which connect selected points of the TMDC
Brillouin zone with the $K$ and $K^{\prime}$ points of the graphene Brillouin
zone. As these momentum points vary with the twist angle, the proximity SOC
gets modified. The reliability of such calculations depends on how well both
the energy dispersion of the TMDC layer and the interlayer hybridization
between graphene and TMDC orbitals are described by the effective microscopic
models.

Perhaps the most direct theoretical approach for extracting proximity SOC in
graphene/TMDC bilayers is to perform ab initio calculations on supercell
geometries incorporating the twisted monolayers. Since density functional
theory (DFT) calculations in this regard have been performed almost
exclusively on aligned heterostructures (0° twist angle), there is yet no
clear ab initio perspective on this topic. There already exist studies (Wang
_et al._ , 2015a; Hou  _et al._ , 2017; Felice  _et al._ , 2017) discussing
electronic properties of twisted graphene/TMDC stacks, but not examining the
proximity SOC. Recently though, Pezo et al. (Pezo  _et al._ , 2021) reported
SOC parameters from DFT calculations on three twist angles (0°, 15°, and 30°)
for graphene/MoTe2 and graphene/WSe2, and examined ramifications for spin
relaxation anisotropy.

Here, we aim to provide a comprehensive picture, applying DFT to large
graphene/TMDC heterostructures employing four different semiconducting TMDCs
MoS2, MoSe2, WSe2, and WS2. For several twist angles we build commensurate
supercells with built-in strains in graphene (TMDCs are left unstrained) of up
to 17%, although we deem as excessive strains of more than 10%. From the band
structures of the supercells we find that the apparent dependence of the band
offsets between the Dirac point and the TMDC bands is due to the strain. In
other words, heterostructures of different twist angles but similar strain
exhibit similar band offsets. By plotting the band offsets as a function of
the strain we find a scattered linear relation which allows us to extract the
“zero-strain” offset. To correct for the strain we then apply a transverse
electric field, different for different heterostructures, so that at the end
we compare different extracted parameters for the same band offsets. We
believe that this currently is the best way to extract quantitatively reliable
twist-angle dependence of the proximity spin-orbit parameters for
computationally feasible supercell sizes. We compare our results with previous
tight-binding studies (Li and Koshino, 2019; David  _et al._ , 2019).

Overall, the magnitude of the SOC parameters are on the meV scale, while the
Rashba coupling is typically weaker than the dominant valley-Zeeman coupling.
Only close to 30$\degree$, at which the valley-Zeeman coupling vanishes due to
symmetry, does Rashba coupling prevail. There is also a marked difference for
Mo- and W-based heterostructures. For Mo-based stacks the valley-Zeeman
coupling exhibits a global maximum at a peak around 20$\degree$, reaching a
value twice that at 0$\degree$. In contrast, while there is a local maximum
for W-based stacks close in the region of 15$\degree-20\degree$, the largest
value of the valley-Zeeman coupling is found at 0$\degree$. The largest
parameters for W-based structures are roughly twice as large as the largest
ones for Mo-structures.

Unlike at 0$\degree$ and 30$\degree$, at a general twist angle between the two
extremes (outside of this interval the parameters can be obtained from
symmetry arguments as presented below) the Rashba coupling exhibits a radial
component which results in a spin texture which deviates from the typical
tangential Rashba pattern. The deviation is measured by the Rashba phase angle
Li and Koshino (2019); David  _et al._ (2019). We also report these angles for
selected supercells and demonstrate their rather high sensitivity to the
electric field, but not to lateral shifts of the layers.

There certainly are still many pitfalls of the methodology we use. Finite
supercells are susceptible to the atomic registry and atomic relaxation
leading to graphene rippling. We also investigate the dependence of the
proximity spin-orbit coupling parameters on such effects, although in a
limited way due to computational complexity, focusing on selected
heterostructures to provide realistic expectations rather than sweeping
proofs. For example, we find that lateral shifts of the atomic layers do not
change significantly the spin-orbit parameters, while the rippling of graphene
strongly enhances the staggered potential and increases the Kane-Mele
coupling, though valley-Zeeman coupling still dominates. Also, the
compensating electric field to correct the band offsets is not the panacea.
Especially for excessive strain above 10% the electric field can modify the
band structure itself and induce Rashba spin-orbit couplings that are not
native to the heterostructures. We provide a discussion of the electric field
effects as well.

The paper is organized as follows. In Sec. II we introduce supercell
geometries, the computational methods are detailed in App. A. The influence of
strain on the band offsets and our method of correcting them is presented in
Sec. III. In Sec. IV we introduce the effective low energy Dirac Hamiltonian,
which is used to extract orbital and proximity SOC parameters. Finally, in
Sec. V we discuss the extracted parameters for different angles. We
additionally comment on the effect of lateral shifts and structural relaxation
in App. 6 and App. C respectively. In App. 7 we quantify the Rashba phase
angle for different TMDCs and twist angles. Using a specific example, we
demonstrate in App. E how a transverse electric field tunes the extracted
spin-orbit parameters.

##  II Supercell geometries

Starting from the primitive hexagonal unit cells of graphene and four TMDCs
(MoS2, MoSe2, WSe2 and WS2; geometry parameters listed in Tab. 2), we
construct the supercells by implementing the coincidence lattice method (Koda
_et al._ , 2016; Carr  _et al._ , 2020; Wang  _et al._ , 2015b). Using the
integers $n$ and $m$, we define the lattice vectors $\mathbf{a}^{S}_{(n,m)}$
and $\mathbf{b}^{S}_{(n,m)}$ of a new supercell as a linear combination of the
primitive lattice vectors $\mathbf{a}$ and $\mathbf{b}$:

$\displaystyle\mathbf{a}^{S}_{(n,m)}$
$\displaystyle=n\cdot\mathbf{a}+m\cdot\mathbf{b}$ (1)
$\displaystyle\mathbf{b}^{S}_{(n,m)}$
$\displaystyle=-m\cdot\mathbf{a}+(n+m)\cdot\mathbf{b}.$ (2)

We give attributes $(n,m)$ to such a supercell. Its lattice constant is

$\displaystyle a^{S}_{(n,m)}=a\cdot\sqrt{(n^{2}+m^{2})+n\cdot m},$ (3)

where $a=|\mathbf{a}|=|\mathbf{b}|$ is the lattice constant of the primitive
unit cell. The relative twist angle with respect to the primitive unit cell is
given by

$\displaystyle\Theta_{(n,m)}=\arctan(\frac{\sqrt{3}m}{2n+m}).$ (4)

A graphene/TMDC heterostructure supercell contains a ($n$,$m$) graphene
supercell beneath a ($n^{\prime}$,$m^{\prime}$) TMDC supercell resulting in a
relative twist angle

$\displaystyle\Theta=$
$\displaystyle\Theta_{(n,m)}-\Theta_{(n^{\prime},m^{\prime})}$
$\displaystyle=$
$\displaystyle\arctan\Big(\frac{\sqrt{3}m}{2n+m}\Big{missing})-\arctan\Big(\frac{\sqrt{3}m^{\prime}}{2n^{\prime}+m^{\prime}}\Big{missing}).$

Table 1:  Structural information of the investigated graphene/TMDC heterostructures. We list the supercell attributes $(n,m)$ of graphene and $(n^{\prime},m^{\prime})$ of the TMDC. We also list the absolute value of the twist angle $|\Theta|$ between the two monolayers (the sign of $\Theta$ is not relevant here), the strain $\epsilon^{\text{TMDC}}$ in graphene (which depends on the specific choice of TMDC) and the number of atoms ($N_{\text{at}}$) in the heterostructure. For completeness, we also list strains corresponding to supercells, which were not investigated, in grey - for a full list of the supercells used in Fig. 5 and Fig. 2, see Tab. 4. These supercells, which were not investigated, either had too much built-in strain, too many atoms or entailed computational difficulties (e.g., convergence problems). $|\Theta|$ | $(n,m)$ | $(n^{\prime},m^{\prime})$ | $\epsilon^{\text{MoS}_{2}}$ | $\epsilon^{\text{WS}_{2}}$ | $\epsilon^{\text{MoSe}_{2}}$ | $\epsilon^{\text{WSe}_{2}}$ | $N_{\text{at}}$
---|---|---|---|---|---|---|---
$[\degree]$ |  |  | [%] | [%] | [%] | [%] |
0.0 | (4,0) | (3,0) | -2.9 | -3.05 | 1.19 | 1.19 | 59
0.0 | (0,5) | (0,4) | 3.58 | 3.41 | 7.93 | 7.93 | 98
5.2 | (3,1) | (2,1) | -4.99 | -5.14 | -1.0 | -1.0 | 47
6.6 | (3,2) | (2,2) | 2.89 | 2.73 | 7.22 | 7.22 | 74
9.5 | (3,2) | (3,1) | 7.1 | 6.93 | 11.6 | 11.6 | 77
10.9 | (2,1) | (1,1) | -15.24 | -15.37 | -11.67 | -11.67 | 23
13.9 | (3,1) | (3,0) | 7.73 | 7.56 | 12.26 | 12.26 | 53
13.9 | (0,4) | (1,3) | 16.7 | 16.52 | 21.61 | 21.61 | 71
13.9 | (5,0) | (3,1) | -6.64 | -6.78 | -2.71 | -2.71 | 89
19.1 | (2,1) | (2,0) | -2.13 | -2.28 | 1.99 | 1.99 | 26
22.7 | (3,2) | (1,3) | 7.1 | 6.93 | 11.6 | 11.6 | 77
23.4 | (3,2) | (3,0) | -10.89 | -11.03 | -7.14 | -7.14 | 65
27.0 | (3,1) | (1,2) | -4.99 | -5.14 | -1.0 | -1.0 | 47
30.0 | (2,0) | (1,1) | 12.13 | 11.95 | 16.84 | 16.84 | 17
30.0 | (5,0) | (2,2) | -10.3 | -10.44 | -6.53 | -6.53 | 86
Table 2: Unstrained geometries of the primitive unit cells of graphene and four selected TMDCs (Zollner  _et al._ , 2019). The structure of the TMDCs stay unchanged in the supercells, while the graphene layers are strained by the factors listed in Tab. 1 to ensure commensurability. | Graphene | MoS2 | WS2 | MoSe2 | WSe2
---|---|---|---|---|---
$a$[Å] | 2.46 | 3.185 | 3.18 | 3.319 | 3.319
$d_{\text{XX}}$[Å] | - | 3.138 | 3.145 | 3.357 | 3.364
Figure 1: (a)-(c) Bottom view of the graphene/TMDC heterostructure supercells
with twist angles $|\Theta|=0\degree,19.1\degree,30\degree$ (d) Side view of
the 0° unit cell with indicated interlayer distance $d$ and chalcogen-
chalcogen distance $d_{\text{XX}}$.

In our supercells there is at a corner a carbon atom directly beneath a metal
atom of the TMDC, see Fig. 1. Considering different configurations (see App.
6) we find that proximity SOC is rather insensitive to the changes of the
atomic registry, as was already shown for $0\degree$ in Ref. Gmitra  _et al._
(2016). In order to obtain commensurate supercells for periodic DFT
calculations, one of the layers (or both) need to be strained. We thus
introduce the strain factor $\epsilon$ which depends on the twist angle and
the lattice constant of the TMDC. TMDCs are very sensitive to strain (Zollner
_et al._ , 2019), while the low energy Dirac spectrum of graphene is—apart
from a renormalization of the Fermi velocity—rather robust against biaxial
strain smaller than 20% (Si  _et al._ , 2016; Choi  _et al._ , 2010).
Therefore, we choose to leave the TMDC unstrained and instead strain graphene.
Also, we use the same interlayer distance $d=3.3$ Å separating the monolayers
(see App. C) for all studied supercells, to focus on twist-angle effects.
However, in App. C we discuss the effects of structural relaxation and
rippling on the proximity band structure at the Dirac cone. Finally, we add a
vacuum of 20 Å to avoid interactions between periodic images in our slab
geometry. The graphene/TMDC heterostructures are set up using the atomic
simulation environment (ASE) (Bahn and Jacobsen, 2002) code. The structural
parameters of the heterostructures are collected in Tab. 1 and some prominent
examples are visualized in Fig. 1.

When a hexagonal system is described by a $(n,m)$ supercell, the K point of
the Brillouin zone can fold back to either K, K’, or $\Gamma$ of the reduced
Brillouin zone. The following rule, which can be derived from geometrical
considerations, determines, which of the three options is the case:

$\displaystyle\text{backfolding to $\Gamma$ for: }n-m=0+3\cdot l$ (5)
$\displaystyle\text{backfolding to K for: }n-m=1+3\cdot l$ (6)
$\displaystyle\text{backfolding to K' for: }n-m=2+3\cdot l$ (7)

with $l\in\mathbb{Z}$. Because of these backfolding effects, the TMDC band gap
resides at the K- or the $\Gamma$-point depending on the exact supercell. The
Dirac cone can in principle also fold back to $\Gamma$, e.g., for a
$(n,m)=(1,1)$ supercell (i.e., a $\sqrt{3}\crossproduct\sqrt{3}$ supercell).
In the following we only consider supercells for which the Dirac cone folds
back to the K/K’-point to be able to compare different structures and extract
the spin-orbit parameters using our model Hamiltonian presented below.

Computational methodology for obtaining DFT electronic band structures of the
supercells is detailed in App. A.

##  III Correcting band offsets for strain

For most investigated twisted heterostructures we find that the Dirac cone of
graphene lies within the TMDC band gap for most supercells. Exceptions are
stacks with heavily strained graphene.

Figure 2: Correcting for the strain. For all the investigated supercells of
graphene with MoS2, WS2, MoSe2, and WSe2 monolayers we plot the band offsets
$E_{C}$ of the Dirac cone with respect to conduction band (see inset) against
the strain on graphene $\epsilon$; $\epsilon>0$ indicates tensile strain while
$\epsilon<0$ indicates compressive strain. Each of the data points (red solid
circles) is annotated with the twist angle of the corresponding supercell.
From the linear fit (red line) we extract the (apparent) zero-strain band
offsets (empty red circles) which are collected in Tab. 3. The green circles
show the band offsets after the correction by the transverse electric field
employed to compensate the influence of strain. Strains above 10$\%$ and below
-10$\%$ and negative band offsets $E_{C}<0$ were not included in the linear
fit. Figure 3: Calculated band structures along high symmetry lines for
selected graphene/TMDC supercells. For each of the four TMDCs, we show three
different angles: $|\Theta|=0\degree,19.1\degree,30\degree$, as indicated. The
grey circles originate from TMDC orbitals, while the solid (dashed) lines come
from graphene states with (without) the electric field corrections. The energy
scale is set to zero for the Fermi energy of the electric-field corrected band
structures. Because of band folding effects, the TMDC bands look different,
while the global TMDC band gap stays intact for different twist angles.

Crucially, the relative twist angle $\Theta$ barely influences the band offset
between the Dirac point and the TMDC conduction (valence) band minimum
(maximum). Instead, we trace the considerable band offset variations among the
investigated heterostructures to the applied strain which is otherwise
necessary to build commensurate supercells of manageable sizes. In Fig. 2 we
collect the data for the investigated structures and plot the band offset
$E_{C}$ to the conduction band as a function of the strain $\epsilon$ applied
to graphene. The band offset varies, with some scatter, linearly with
$\epsilon$. The above thesis, that it is the strain and not twist that causes
the band offset variation, is most strikingly seen in the data of disparate
twist angles and very similar strain, such as the three angles 9.5$\degree$,
22.7$\degree$, and 13.9$\degree$ for graphene/MoS2, see Fig. 2. These three
heterostructures feature essentially the same offset $E_{C}$.

The absolute deformation potential of the Dirac cone can be defined as the
change of its energy when applying strain. Using the TMDCs band edges (which
remain at the same energy, because we leave the TMDC unstrained) as a fix
point, we can identify $|\alpha_{D}|$ as the absolute value of the slope of
the linear fits in Fig. 2:

$\displaystyle\alpha_{D}=\frac{\partial
E_{D}}{\partial\epsilon}=\frac{\partial
E_{V}}{\partial\epsilon}=-\frac{\partial E_{C}}{\partial\epsilon}$ (8)

Here $E_{D}$ is the energy level of the Dirac cone and $E_{V}$ and $E_{C}$ are
the band offsets to the valence and conduction band respectively. We find
$\alpha_{D}$ to be roughly -80 meV/%, i.e. the Dirac cone is lowered in energy
towards the valence band by 80 meV for each percent of tensile strain on the
graphene. Using this relation we were able to extrapolate an estimate for the
zero strain band offset for all four TMDCs (see Tab. 3).

Table 3: Extrapolated zero-strain valence and conduction band offsets (with respect to the graphene Dirac cone) $E_{V}^{0}$ and $E_{C}^{0}$, respectively, and the absolute deformation potential $\alpha_{D}$, as obtained by the fit in Fig. 2. Data points with excessive strain ($|\epsilon|>10\%$) or negative band offset $E_{V/C}<0$ were excluded from the fit. | MoS2 | WS2 | MoSe2 | WSe2
---|---|---|---|---
$E_{V}^{0}$[meV] | 1365 | 1027 | 680 | 290
$E_{C}^{0}$[meV] | 238 | 546 | 673 | 978
$\alpha_{D}$ [meV/%] | -74 | -79 | -78 | -89

The large magnitude of the extracted deformation potential indicates that even
small strains can cause large band offset changes. From the perspective of
perturbation theory, it is clear that the energy distance between the Dirac
cone and the TMDC bands, i.e., the bands offsets, influences the proximity SOC
(David  _et al._ , 2019) (see App. E). Therefore, to obtain reliable proximity
orbital and spin-orbit parameters the offsets need to be corrected. In order
to bring the Dirac cone to its apparent zero-strain level, as specified by the
obtained offsets in Tab. 3, we apply a transverse (perpendicular to the
layers) electric field to compensate the effects of strain $\epsilon$. This is
done in Fig. 2. The values of the compensating electric fields are listed in
Tab. 4. A field is defined as positive if it points from the TMDC layer to the
graphene layer. Figure 3 shows the global band structure along high-symmetry
lines for three selected angles, before and after the offset correction. The
Dirac points are located well inside the TMDC band gaps and the proximity
effects can be well discerned, see Fig. 4.

##  IV Effective Hamiltonian

Figure 4: Proximity Dirac band structures (after correcting band offsets for
strain) at K. We show the results for $|\Theta|=0\degree$, $19.1\degree$, and
$30\degree$. The dots represent the DFT data while the solid line represents
the fit from the model Hamiltonian in Eq. (9). The spin-z expectation value is
color coded from $\langle s_{z}\rangle=-0.5$ (spin down; blue) over $\langle
s_{z}\rangle=0$ (unpolarized; grey) to $\langle s_{z}\rangle=0.5$ (spin up;
red). The extracted parameters are listed in Tab. 4. The last row shows the
spin-x and spin-y expectation values of the spin-up (red arrows) and spin-down
(blue arrows) valence band of the Dirac cone of graphene proximitized by MoSe2
for the three twist angles. The k-path (dotted line) goes along a circular
path around the K-point. We also present the extracted value for the Rashba
angle $\Phi$ which measures the uniform deviation of the spin texture from
tangential.

To find the twist-angle dependence of the proximity induced SOC in our
heterostructures, we fit the DFT band structures at the Dirac points to a
model Hamiltonian (Gmitra and Fabian, 2015). The Hamiltonian $H$ comprises the
orbital part $H_{\text{orb}}$ and the spin-orbit part $H_{\text{so}}$, which
is composed of the intrinsic spin-orbit coupling $H_{\text{so,I}}$ and the
Rashba coupling $H_{\text{so,R}}$:

$H(\mathbf{k})=H_{\text{orb}}(\mathbf{k})+H_{\text{so}}=H_{\text{orb}}(\mathbf{k})+H_{\text{so,I}}+H_{\text{so,R}}.$
(9)

The orbital part describes the dispersion of the graphene Dirac cone
linearized around the K/K’-point; accordingly, $\mathbf{k}$ is the electron
wave vector measured from K/K’. It also includes a staggered potential
$\Delta$, caused by the substrate’s asymmetrical influence on the graphene A-
and B-sublattice:

$H_{\text{orb}}(\mathbf{k})=\hbar
v_{F}(\kappa\sigma_{x}k_{x}+\sigma_{y}k_{y})+\Delta\sigma_{z}.$ (10)

Here, $v_{F}$ is the Fermi velocity of the Dirac electrons and $\sigma$ are
the Pauli matrices operating on the sublattice (A/B) space. The parameter
$\kappa=1$ for K and $\kappa=-1$ for K’.

The intrinsic spin-orbit Hamiltonian

$H_{\text{so,I}}=\frac{1}{2}\Big{[}\lambda_{\text{I}}^{A}(\sigma_{z}+\sigma_{0})+\lambda_{\text{I}}^{B}(\sigma_{z}-\sigma_{0})\Big{]}\kappa
s_{z},$ (11)

and the Rashba spin-orbit Hamiltonian

$H_{\text{so,R}}=-\lambda_{\text{R}}\exp(-i\Phi\frac{s_{z}}{2})\Big{[}\kappa\sigma_{x}s_{y}-\sigma_{y}s_{x}\Big{]}\exp(i\Phi\frac{s_{z}}{2}),$
(12)

both additionally act on the spin space, described by the spin Pauli matrices
$s$; $\lambda_{\text{I}}^{A}$ and $\lambda_{\text{I}}^{B}$ are the intrinsic
spin-orbit parameters for the sublattices A and B respectively, while
$\lambda_{\text{R}}$ is the Rashba SOC strength. The Rashba phase angle
$\Phi$—present in $C_{3}$ symmetric structures Li and Koshino (2019); David
_et al._ (2019)—rotates the spin texture about the $z$-axis, adding a radial
component to the Rashba field. For $\Theta=0\degree$ and $\Theta=30\degree$,
the reflection symmetry (Li and Koshino, 2019) along the $x$\- and $y$-axis
respectively forces either $\Phi=0$ or $\Phi=180\degree$, therefore this
parameter does not appear in the Hamiltonian of Ref. Gmitra and Fabian (2015)
valid for aligned structures.

We choose to limit the Rashba parameter to positive values
$\lambda_{\text{R}}>0$. A sign change of $\lambda_{\text{R}}$ then corresponds
to an additional phase shift of $\Phi$ by a half rotation, i.e.
$\Phi\rightarrow\Phi+180\degree$. To make this clear we always list
$|\lambda_{\text{R}}|$ below. We also evaluated the angles $\Phi$ for a few
selected supercells (see Fig. 4 and App. 7). While the energy dispersion is
not affected by $\Phi$, the radial component of the Rashba field can affect
spin physics and need to be considered when interpreting experimental results
on spin transport and spin-orbit torque. For example, the in-plane spin
accumulation in the presence of electrical current can have a component along
the current, unlike for the usual Rashba effect.

We also define

$\displaystyle\lambda_{\text{VZ}}=\frac{\lambda_{\text{I}}^{A}-\lambda_{\text{I}}^{B}}{2}\text{\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ and\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
}\lambda_{\text{KM}}=\frac{\lambda_{\text{I}}^{A}+\lambda_{\text{I}}^{B}}{2}$
(13)

as the valley-Zeeman Gmitra and Fabian (2015); Wang  _et al._ (2015a) SOC
(sublattice-odd) and the Kane-Mele Kane and Mele (2005a) SOC (sublattice-even)
respectively. It turns out that $\lambda_{\text{KM}}$ is negligible for
graphene/TMDC heterostructures Gmitra and Fabian (2015); Wang  _et al._
(2015a). This is also true for the twisted heterostructures presented below,
and already predicted by tight-binding modeling David  _et al._ (2019); Li and
Koshino (2019).

Our DFT structures are for angles between 0° and 30°. The results for all
other twist angles can be obtained by symmetry as follows. Twisting clockwise
or counterclockwise from 0° does not influence the parameters:

$\displaystyle\lambda_{\text{VZ}}(-\theta)$
$\displaystyle=\lambda_{\text{VZ}}(\theta)$ (14)
$\displaystyle\lambda_{\text{R}}(-\theta)$
$\displaystyle=\lambda_{\text{R}}(\theta)$ (15) $\displaystyle\Delta(-\theta)$
$\displaystyle=\Delta(\theta).$ (16)

Additionally a twist by 60° corresponds to switching the sublattices of
graphene and therefore changes the sign of the sublattice-sensitive
parameters:

$\displaystyle\lambda_{\text{VZ}}(\theta+60\degree)$
$\displaystyle=-\lambda_{\text{VZ}}(\theta)$ (17)
$\displaystyle\lambda_{\text{R}}(\theta+60\degree)$
$\displaystyle=\lambda_{\text{R}}(\theta)$ (18)
$\displaystyle\Delta(\theta+60\degree)$ $\displaystyle=-\Delta(\theta).$ (19)

Using the first set of rules, one can infer the parameters’ values for
$\Theta\in[-30\degree,0\degree]$ from the values for
$\Theta\in[0\degree,30\degree]$. Then, using the second set of rules, one can
infer the values for $\Theta\in[30\degree,60\degree]$. For example,
$\lambda_{\text{VZ}}(32\degree)=-\lambda_{\text{VZ}}(-28\degree)=-\lambda_{\text{VZ}}(28\degree)$.
This leads to the conclusion that $\lambda_{\text{VZ}}$ changes sign at
$\Theta=30\degree$ and $\lambda_{\text{VZ}}(30\degree)=0$.

##  V Results

The fitted valley-Zeeman $\lambda_{\text{VZ}}$ and Rashba $\lambda_{\text{R}}$
spin-orbit couplings are collected in Tab. 4 for both uncorrected and
corrected electronic band structures. The latter are plotted in Fig. 5. In
agreement with earlier DFT studies of aligned graphene/TMDC heterostructures
Gmitra and Fabian (2015); Gmitra  _et al._ (2016); Wang  _et al._ (2015a), the
proximity spin-orbit coupling of the Dirac electrons is on the meV scale. In
all cases valley-Zeeman SOC $\lambda_{\text{VZ}}$ vanishes at
$\Theta=30\degree$, while the Rashba SOC $\lambda_{\text{R}}$ reaches a
minimum there. We also find that $\lambda_{\text{VZ}}$ for the Mo-based TMDC
heterostructures exhibit a peak at about $\theta=20\degree$, where the
parameters jump to about twice the value at 0$\degree$. This pronounced peak
is preceded by a slight dip at around $\Theta$ between 10$\degree$ and
15$\degree$. Similarly, for the heterostructures based on WSe2 and WS2 there
appears to be a peak in the valley-Zeeman coupling below 20$\degree$, although
its magnitude is less that the value at 0$\degree$. Except for twist angles
close to 30$\degree$, the Rashba parameters are generally smaller in magnitude
than $\lambda_{\text{VZ}}$.

Tight-binding models Li and Koshino (2019); David  _et al._ (2019) predicted a
peak structure with a global maximum in the region of 15$\degree$-20$\degree$.
We can confirm this picture for Mo-based but not for W-based heterostructures.
The two tight-binding models differ in their predictions for the sign of
$\lambda_{\text{VZ}}$: Ref. Li and Koshino (2019) predicts positive
$\lambda_{\text{VZ}}$, while Ref. David  _et al._ (2019) predicts a sign
change at around 10$\degree$. As seen in Fig. 5, the extracted valley-Zeeman
coupling parameters do not change sign, in agreement with Ref. Li and Koshino
(2019).

To check our approach, for 0° —and for MoS2 and MoSe2 additionally for
13.9°—we plot in Fig. 5 two data points corresponding to two different
supercells with the same twist angle but different strain. Such data pairs
indicate how reliable the DFT approach based on strained supercells is.
Ideally, the extracted parameters would agree for the given twist angle after
correcting for strain by transverse electric field. This is indeed the case at
0° where data pairs for TMDCs agree rather well. However, for the 13.9° data
pairs the parameters differ significantly. At this angle one of the
heterostructures in the pairs have strain above 10%. We believe that such a
strain is already too large for obtaining quantitatively reliable proximity
band structures even after correcting by the electric field. We exclude such
data from considering the trends seen in Fig. 5. The reason for the
sensitivity to strain is that the size of the Brillouin zone of strained
graphene varies with $\epsilon$. Following the interlayer coupling picture of
Refs. (Koshino, 2015; David  _et al._ , 2019), the graphene K-point thus
couples to somewhat different regions of the TMDC Brillouin zone. This
momentum-specific coupling (hybridization of the corresponding graphene and
TMDC orbitals) is responsible for the proximity effect.

Figure 5: The E-field corrected spin-orbit parameters $\lambda_{\text{R}}$
(red squares) and $\lambda_{\text{VZ}}$ (blue dots) of the graphene Dirac cone
as function of the twist angle $\Theta$ for the heterostructures with TMDCs
MoS2, MoSe2, WSe2, and WS2, as indicated. We de-emphasize (transparent, light
red) data points from supercells with excessive built-in strain of
$\epsilon>10\%$, deeming them less reliable. The dotted line is a guide for
the eyes.

We note that the differences in our reported data in Tab. 4 and previous ones
Gmitra  _et al._ (2016) at 0$\degree$ are due to the unrelaxed atomic
structures used here. The main effect of atomic relaxation is to introduce
rippling of the graphene layer which enhances especially the staggered
potential $\Delta$ and the Kane-Mele SOC $\lambda_{\text{KM}}$ parameter. In
our calculations both $\Delta$ and $\lambda_{\text{KM}}$ are rather
suppressed. We have studied the atomic relaxation effects for selected
heterostructures with 0 and 19.1$\degree$ twist angles; the results are
presented in App. C).

Finally, let us compare the proximity SOC for uncorrected and corrected band
structures. The former is plotted in Fig. 9, the latter in Fig. 5 which we
discussed above. We see that the while the overall magnitudes of
$\lambda_{\text{VZ}}$ and $\lambda_{\text{R}}$ are similar for both cases, the
correction by electric field has an influence on the overall twist-angle
dependence. The main effect of the correction is shifting the band offsets
which modify the interlayer coupling and charge transfer between the layers.
Both have an effect on the valley-Zeeman and Rashba couplings. In App. E we
further investigate the influence of the electric field on the proximity spin-
orbit parameters for a 19.1$\degree$ graphene/MoS2 heterostructure,
demonstrating a rather strong decrease of the magnitude of
$\lambda_{\text{VZ}}$ upon increasing of E-field and accompanying increase of
the band offset $E_{C}$, see Fig. 8.

##  VI Summary

We have performed systematic investigations of the twist-angle dependence of
the proximity spin-orbit coupling induced by four TMDC monolayers MoS2, MoSe2,
WS2, and WSe2 in graphene. For each graphene/TMDC bilayer we have studied
several twist angles, being computationally limited by the supercell sizes. We
have identified strain, which is necessary to build commensurate supercells of
reasonable sizes, as the key factor influencing the reliability of the
extracted SOC parameters. To correct for strain, we apply a transverse (to the
layers) electric field which compensates the changes of the band offsets
caused by strain.

Extracting the SOC parameters by fitting the corrected DFT band structures to
an already established model Hamiltonian, we present as the main result the
magnitudes of the valley-Zeeman and Rashba couplings of the Dirac electrons as
functions of the twist angles from 0$\degree$ to 30$\degree$; symmetry
relations can be used to deduce the parameters for other angles. A
qualitatively new effect of twisting is the emergence of the Rashba phase
angle $\Phi$ which rotates the spin texture away from tangential; $\Phi$
vanishes at 0$\degree$ and 30$\degree$. There appears to be no discernible
trend for $\Phi$ as a function of the twist angle, but we show that the phase
angle depends strongly on the electric field. While not affecting the energy
dispersion, nonzero $\Phi$ adds a radial component to the typical tangential
Rashba texture and should be considered when interpreting spin transport
experiments in twisted heterostructures. In particular, the spin accumulation
due to electric current need not be perpendicular to the current.
Unfortunately, the phase angles appear too sensitive to the band offsets and
transverse electric fields, preventing us to make reliable quantitative
predictions on their magnitudes. It is natural to assume that each twisted
heterostructure is in this sense unique, exhibiting $\Phi$ based on the twist
angle, unintended doping, interlayer distance, and applied gate. An
experimental observation of $\Phi$ and its tunability is currently
outstanding.

As required by symmetry, the valley-Zeeman SOC vanishes at 30$\degree$. This
clear case should serve as a useful tool for experiments relying on valley-
Zeeman coupling. For other twist angles, we can conclude that the coupling
does not change sign between 0$\degree$ and 30$\degree$, and that
$\lambda_{\text{VZ}}$ for Mo-based TMDCs has a (global) maximum at roughly
20°. For W-based heterostructures the maximum coupling appears to be at
0$\degree$. The Rashba coupling is, in general, weaker than the valley-Zeeman
coupling, except close to 30$\degree$ where the latter vanishes.

We have also studied the effects of structural relaxation on the proximity
effects. We found that rippling of graphene is the main factor significantly
enhancing the staggered potential $\Delta$ and the Kane-Mele coupling
$\lambda_{\text{KM}}$ (which is otherwise negligible for unrelaxed
structures), but does not affect strongly the valley-Zeeman and Rashba
couplings. Relative lateral shifts of the graphene and TMDC layers too do not
have a significant effect on the proximity spin-orbit parameters.

###### Acknowledgements.

This work was funded by the International Doctorate Program Topological
Insulators of the Elite Network of Bavaria, the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) SFB 1277 (Project-ID
314695032), SPP 2244 (project no. 443416183), and by the European Union
Horizon 2020 Research and Innovation Program under contract number 881603
(Graphene Flagship). M.G. acknowledges VEGA 1/0105/20.

Table 4: Parameters extracted from the band structure calculations. For all four TMDCs and all used angles, we list the band offset $E_{V}$ ($E_{C}$) of the Dirac cone with respect to the valence band (conduction band) and the extracted model Hamiltonian (see Eq. (9)) parameters (staggered potential $\Delta$, intrinsic SOC for site A and B $\lambda_{\text{I}}^{A}$ and $\lambda_{\text{I}}^{B}$ and Rashba SOC $\lambda_{\text{R}}$). We denote the offsets and parameters after correction with the electric field with a bar, for example $\bar{\lambda}_{I}^{A}$. The electric field is defined as positive, if it points from the TMDC layer to the graphene layer. $|\theta|[\degree]$ | $\epsilon$ | $\Delta$ | $\lambda_{\text{I}}^{A}$ | $\lambda_{\text{I}}^{B}$ | $|\lambda_{\text{R}}|$ | $E_{V}$ | $E_{C}$ | E-field | $\bar{\Delta}$ | $\bar{\lambda}_{I}^{A}$ | $\bar{\lambda}_{I}^{B}$ | $|\bar{\lambda}_{R}|$ | $\bar{E}_{V}$ | $\bar{E}_{C}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
| [%] | [meV] | [meV] | [meV] | [meV] | [eV] | [eV] | [$\text{V}/\text{nm}$] | [meV] | [meV] | [meV] | [meV] | [eV] | [eV]
MoS2 |  |  |  |  |  |  |  |  |  |  |  |  |  |
0 | -2.9 | 0.014 | 0.346 | -0.346 | 0.085 | 1.598 | 0.004 | 1.499 | 0.009 | 0.331 | -0.331 | 0.081 | 1.341 | 0.262
0 | 3.58 | -0.005 | 0.402 | -0.408 | 0.043 | 1.066 | 0.533 | -1.748 | 0.130 | 0.445 | -0.442 | 0.016 | 1.358 | 0.236
5.2 | -4.99 | 0.022 | 0.455 | -0.457 | 0.165 | 1.773 | -0.169 | - | - | - | - | - | - | -
6.6 | 2.89 | 0.054 | 0.341 | -0.341 | 0.104 | 1.118 | 0.482 | -1.301 | 0.032 | 0.384 | -0.384 | 0.111 | 1.334 | 0.264
9.5 | 7.1 | 0.023 | 0.225 | -0.225 | 0.139 | 0.838 | 0.761 | -2.962 | 0.026 | 0.356 | -0.358 | 0.146 | 1.322 | 0.262
13.9 | 7.73 | 0.017 | 0.119 | -0.120 | 0.168 | 0.809 | 0.791 | -3.183 | 0.010 | 0.291 | -0.293 | 0.165 | 1.331 | 0.258
13.9 | 16.7 | 0.037 | -0.270 | 0.270 | 0.265 | 0.377 | 1.230 | -6.154 | -0.327 | 0.176 | -0.174 | 0.157 | 1.381 | 0.256
13.9 | -6.64 | 0.124 | 8.707 | -8.754 | 1.251 | 1.828 | -0.222 | - | - | - | - | - | - | -
19.1 | -2.13 | 0.023 | 1.223 | -1.230 | 0.241 | 1.532 | 0.075 | 1.107 | 0.031 | 0.593 | -0.602 | 0.159 | 1.336 | 0.258
22.7 | 7.1 | -0.005 | -0.040 | 0.040 | 0.047 | 0.859 | 0.745 | -2.564 | -0.004 | 0.052 | -0.052 | 0.092 | 1.242 | 0.315
27 | -4.99 | -0.004 | 0.459 | -0.458 | 0.531 | 1.773 | -0.168 | - | - | - | - | - | - | -
30 | 12.13 | 0.000 | 0.000 | 0.000 | 0.075 | 0.596 | 1.009 | -6.582 | 0.059 | 0.002 | 0.002 | 0.057 | 1.359 | 0.235
30 | -10.3 | 0.000 | -0.062 | -0.062 | 1.774 | 1.868 | -0.418 | - | - | - | - | - | - | -
WS2 |  |  |  |  |  |  |  |  |  |  |  |  |  |
0 | -3.05 | 0.021 | 1.094 | -1.095 | 0.154 | 1.248 | 0.325 | 0.951 | 0.014 | 1.068 | -1.072 | 0.168 | 1.082 | 0.492
0 | 3.41 | -0.097 | 1.201 | -1.193 | 0.077 | 0.717 | 0.853 | -1.969 | -0.047 | 1.32 | -1.312 | 0.024 | 1.053 | 0.514
5.2 | -5.14 | 0.000 | 1.148 | -1.152 | 0.287 | 1.459 | 0.114 | 2.234 | 0.008 | 1.011 | -1.012 | 0.233 | 1.056 | 0.518
13.9 | 7.56 | 0.866 | -0.076 | 0.068 | 0.386 | 0.463 | 1.110 | -3.416 | -0.263 | 0.394 | -0.397 | 0.309 | 1.035 | 0.531
19.1 | -2.28 | 0.004 | 1.305 | -1.316 | 0.526 | 1.185 | 0.389 | 0.662 | 0.005 | 0.938 | -0.947 | 0.454 | 1.067 | 0.507
27 | -5.14 | -0.052 | 0.608 | -0.606 | 0.959 | 1.459 | 0.116 | 2.195 | 0.010 | 0.171 | -0.170 | 0.559 | 1.073 | 0.503
30 | 11.95 | 0.000 | 0.004 | 0.004 | 0.214 | 0.237 | 1.327 | -4.973 | -0.001 | 0.004 | 0.004 | 0.069 | 1.067 | 0.504
MoSe2 |  |  |  |  |  |  |  |  |  |  |  |  |  |
0 | 1.19 | 0.034 | 0.418 | -0.416 | 0.425 | 0.535 | 0.817 | -0.863 | 0.017 | 0.404 | -0.402 | 0.370 | 0.673 | 0.679
0 | 7.93 | -0.625 | 0.282 | -0.271 | 0.717 | 0.098 | 1.252 | -3.754 | 0.101 | 0.463 | -0.457 | 0.358 | 0.667 | 0.678
5.2 | -1 | 0.024 | 0.412 | -0.412 | 0.368 | 0.708 | 0.645 | 0.209 | 0.016 | 0.413 | -0.413 | 0.378 | 0.674 | 0.679
6.6 | 7.22 | 0.069 | -0.023 | 0.030 | 0.834 | 0.142 | 1.209 | -3.453 | -0.006 | 0.311 | -0.310 | 0.466 | 0.672 | 0.676
9.5 | 11.6 | 0.038 | -1.099 | 1.061 | 1.380 | -0.051 | 1.386 | -5.101 | 0.038 | 0.239 | -0.239 | 0.555 | 0.687 | 0.657
10.9 | -11.67 | 0.244 | 4.906 | -4.632 | 1.011 | 1.719 | -0.363 | - | - | - | - | - | - | -
13.9 | 12.26 | 0.029 | -0.994 | 0.923 | 1.621 | -0.069 | 1.433 | -5.150 | -0.051 | 0.285 | -0.281 | 0.676 | 0.675 | 0.672
13.9 | -2.71 | 0.024 | 0.694 | -0.696 | 0.227 | 0.850 | 0.503 | 1.080 | -0.016 | 0.624 | -0.626 | 0.269 | 0.676 | 0.677
19.1 | 1.99 | -0.015 | 0.873 | -0.873 | 0.573 | 0.496 | 0.856 | -1.048 | 0.000 | 0.871 | -0.872 | 0.490 | 0.655 | 0.678
22.7 | 11.6 | -0.013 | -0.028 | 0.029 | 0.881 | -0.026 | 1.370 | -4.958 | 0.033 | 0.360 | -0.360 | 0.411 | 0.696 | 0.655
23.4 | -7.14 | -0.002 | -0.209 | 0.213 | 0.253 | 1.284 | 0.071 | - | - | - | - | - | - | -
27 | -1 | -0.012 | 0.306 | -0.305 | 0.446 | 0.721 | 0.633 | 0.288 | -0.012 | 0.313 | -0.313 | 0.465 | 0.677 | 0.677
30 | 16.84 | 0.000 | 0.005 | 0.005 | 0.509 | -0.178 | 1.544 | -6.515 | -0.236 | 0.003 | 0.003 | 0.121 | 0.675 | 0.655
30 | -6.53 | 0.000 | -0.004 | -0.004 | 0.213 | 1.224 | 0.131 | - | - | - | - | - | - | -
WSe2 |  |  |  |  |  |  |  |  |  |  |  |  |  |
0 | 1.19 | 0.042 | 1.651 | -1.647 | 0.671 | 0.186 | 1.081 | -0.636 | 0.056 | 1.635 | -1.633 | 0.621 | 0.289 | 0.978
0 | 7.93 | 0.219 | 1.537 | -1.519 | 0.671 | -0.018 | 1.282 | -3.445 | 0.084 | 1.728 | -1.722 | 0.445 | 0.286 | 0.977
5.2 | -1 | 0.034 | 1.495 | -1.494 | 0.675 | 0.360 | 0.907 | 0.449 | 0.033 | 1.491 | -1.490 | 0.703 | 0.291 | 0.977
6.6 | 7.22 | 0.090 | -0.161 | 0.169 | 1.219 | -0.164 | 1.429 | -3.137 | 0.073 | 1.052 | -1.050 | 0.740 | 0.285 | 0.979
10.9 | -11.67 | 0.014 | 3.185 | -3.211 | 1.051 | 1.417 | -0.147 | - | - | - | - | - | - | -
13.9 | -2.71 | 0.023 | 1.612 | -1.616 | 0.642 | 0.504 | 0.764 | 1.327 | -0.022 | 1.387 | -1.388 | 0.709 | 0.290 | 0.978
19.1 | 1.99 | 0.009 | 0.640 | -0.631 | 1.051 | 0.138 | 1.123 | -0.920 | 0.014 | 0.836 | -0.835 | 0.950 | 0.294 | 0.974
23.4 | -7.14 | -0.004 | 0.641 | -0.644 | 0.911 | 0.943 | 0.326 | 3.085 | 0.000 | 0.588 | -0.588 | 0.906 | 0.427 | 0.783
27 | -1 | 0.000 | 0.347 | -0.343 | 0.868 | 0.375 | 0.894 | 0.555 | 0.000 | 0.328 | -0.323 | 0.892 | 0.294 | 0.975
30 | 16.84 | 0.001 | 0.018 | 0.019 | 0.198 | -0.349 | 1.633 | -6.065 | 0.003 | 0.009 | 0.009 | 0.225 | 0.313 | 0.901

##  Appendix A Computational details

All electronic structure calculations are performed by density functional
theory (DFT) (Hohenberg and Kohn, 1964) with Quantum ESPRESSO (Giannozzi and
et al., 2009). Self-consistent calculations are carried out with a $k$-point
sampling of $n_{k}\times n_{k}\times 1$. The number $n_{k}$ is listed in Table
5 for all twist angles. We use energy cutoffs listed in Tab. 5 for the scalar
relativistic pseudopotential with the projector augmented wave method (Kresse
and Joubert, 1999) with the Perdew-Burke-Ernzerhof exchange correlation
functional (Perdew  _et al._ , 1996). Graphene’s $d$-orbitals are not included
in the calculations. We used Van der Waals corrections (Grimme, 2006; Grimme
_et al._ , 2010; Barone  _et al._ , 2009).

The electric fields are implemented in the DFT calculations using a sawtooth
potential in $z$-direction within the quasi 2D unit cell. The electric
potential increases linearly in the area of the heterostructure and then falls
rapidly in the vacuum.

Table 5: Computational details: Used charge density cutoff energy $E_{\rho}$, wave function kinetic energy cutoff $E_{\text{wfc}}$ and $k$-grid density (we used a $n_{k}\times n_{k}$ grid) for the calculations listed in Tab. 4. E-field calculations are listed with a tilde over the material. The supercells are in the same order as in Tab. 1. | MoS2 | $\widetilde{\text{MoS}}_{2}$ | WS2 | $\widetilde{\text{WS}}_{2}$ | MoSe2 | $\widetilde{\text{MoSe}}_{2}$ | WSe2 | $\widetilde{\text{WSe}}_{2}$
---|---|---|---|---|---|---|---|---
$E_{\rho}$[Ry] | 55 | 55 | 70 | 70 | 60 | 60 | 65 | 65
$E_{\text{wfc}}$[Ry] | 350 | 350 | 500 | 500 | 350 | 350 | 550 | 550
$n_{k}$ |  |  |  |  |  |  |  |
0° | 21 | 21 | 15 | 15 | 21 | 15 | 21 | 15
0° | 3 | 9 | 3 | 9 | 3 | 9 | 3 | 9
5.2° | 21 | - | 21 | 15 | 21 | 15 | 21 | 21
6.6° | 21 | 21 | - | - | 18 | 12 | 18 | 12
9.5° | 21 | 18 | - | - | 18 | 18 | - | -
10.9° | - | - | - | - | 21 | - | 21 | -
13.9° | 21 | 21 | 3 | 9 | 21 | 15 | - | -
13.9° | 9 | 15 | - | - | - | - | - | -
13.9° | 21 | - | - | - | 6 | 12 | 3 | 3
19.1° | 30 | 30 | 30 | 30 | 21 | 21 | 21 | 21
22.7° | 21 | 18 | - | - | 18 | 12 | - | -
23.4° | - | - | - | - | 18 | - | 21 | 12
27° | 21 | - | 15 | 15 | 21 | 15 | 21 | 21
30° | 39 | 39 | 21 | 36 | 21 | 21 | 21 | 21
30° | 21 | - | - | - | 6 | - | - | -

##  Appendix B Effect of lateral shift between the layers

For incommensurate heterostructures consisting of two materials whose lattice
constants are no integer multiples of each other, the lateral shifting degree
of freedom does not exist. If the sample can be assumed to be infinite in
$x$\- and $y$-direction, every shifting configuration exists somewhere on the
sample. The physical properties of the configurations will then average out,
when considering the properties of the whole material.

However, the structures we consider in our DFT calculations are commensurate
to be computationally viable and the lattice constants are forced by strain to
be compatible. Therefore, different lateral relative shifts might ensue
different physical properties including different proximity induced SOC.
Naturally, this effect is less relevant for larger supercells, for which the
averaging over the different shifts can occur within the supercell.

In order to estimate the effects of lateral shifts we look at different
configurations for one of our supercells: the 19.1° supercell of MoS2. As one
can see in Fig. 6, all three supercells have very similar spin-orbit
parameters and only the staggered potential $\Delta$ changes. While we cannot
make a sweeping conclusion based on one twist angle and a few lateral shifts
(and computational power limits one’s capabilities here), these results
indicated that the sizes of our supercells are already sufficient to yield
quantitatively reasonable results for the magnitudes of the proximity valley-
Zeeman and Rashba spin-orbit couplings.

Figure 6: Three different options for lateral shifting for the MoS2 19.1°
supercell and their (electric field corrected) model fit parameters.

##  Appendix C Effects of structural relaxation

All the results in the main text were obtained without prior relaxation of the
atomic structure. To study how atomic relaxation affects proximity SOC, we
performed calculations on selected supercells which were relaxed in different
ways. The results are gathered in Tab. 6. Our main finding is that, while the
valley-Zeeman SOC $\lambda_{\text{VZ}}$ is largely unaffected, the staggered
potential $\Delta$ and the Kane-Mele SOC $\lambda_{\text{KM}}$, which are both
very small for the unrelaxed structures, are greatly enhanced by the
relaxation. This can be traced to the rippling of graphene, which results from
the relaxation. We checked this by flattening the relaxed graphene, which
exhibited again very small $\Delta$ and $\lambda_{\text{KM}}$. Our
calculations indicate that even very small ripplings can induce a staggered
potential $\Delta$ on the meV scale. Since the unrelaxed structures show
almost zero $\Delta$ or $\lambda_{\text{KM}}$ (on the order of few tens of
$\mu$eV), the rippling can be considered the sole cause for the sublattice-
asymmetric parameters $\Delta$ and $\lambda_{\text{KM}}$ to be non-zero in
graphene/TMDC heterostructures..

Table 6: Relaxation and its effect on the Dirac cone. We examine the 0° and the 19.1° supercells of MoS2 and the 19.1° supercell of WSe2. We additionally consider the 19.1° MoS2 supercell, with strained TMDC instead of strained graphene (marked with *). For each of those cases we distinguish between calculations without prior relaxation (’fixed’), with prior relaxation (’relaxed’) and with prior relaxation, but with flattened Graphene (’flat Gr’). The latter was achieved by setting the C atom’s $z$-coordinate to an average value after the relaxation. For each case we list the important structural parameters (rippling, interlayer distance $d$ and chalcogen-chalcogen distance $d_{\text{XX}}$) and the extracted parameters from the model fit ($\Delta$, $\lambda_{A}$, $\lambda_{B}$ and $\lambda_{\text{R}}$). All calculations were performed without the electric field corrections.  | rippling[mÅ] | d[Å] |  $d_{\text{XX}}[$Å] |  $\Delta$[meV] |  $\lambda_{A}$[meV] |  $\lambda_{B}$[meV] |  $\lambda_{\text{R}}$[meV]
---|---|---|---|---|---|---|---
MoS-0° |  |  |  |  |  |  |
fixed | 0.00 | 3.30 | 3.138 | 0.014 | 0.35 | -0.35 | 0.085
relaxed | 368.20 | 3.11 | 3.128 | 2.985 | 0.33 | -0.14 | 0.206
flat Gr | 0.00 | 3.33 | 3.128 | 0.003 | 0.30 | -0.30 | 0.074
MoS-19.1° |  |  |  |  |  |  |
fixed | 0.00 | 3.30 | 3.138 | 0.023 | 1.22 | -1.23 | 0.241
relaxed | 23.63 | 3.31 | 3.126 | -1.016 | 1.38 | -0.96 | 0.252
flat Gr | 0.00 | 3.31 | 3.126 | -0.175 | 1.15 | -1.10 | 0.229
MoS-19.1°* |  |  |  |  |  |  |
fix | 0.00 | 3.30 | 3.138 | 0.018 | 0.86 | -0.87 | 0.223
relaxed | 8.64 | 3.37 | 3.079 | -0.294 | 0.59 | -0.53 | 0.163
flat Gr | 0.00 | 3.31 | 3.312 | -0.082 | 0.65 | -0.64 | 0.198
WSe-19.1 |  |  |  |  |  |  |
fix | 0.00 | 3.30 | 3.364 | 0.009 | 0.64 | -0.63 | 1.051
relaxed | 6.56 | 3.37 | 3.347 | 0.295 | 0.60 | -0.51 | 0.763

##  Appendix D Rotation of the Rashba term: Fit results

As mentioned in the main text, to accurately describe the spin-x and spin-y
expectation values of graphene in twisted heterostructures without a mirror
reflection symmetry David  _et al._ (2019); Li and Koshino (2019), one needs
to implement a rotation of the Rashba term $H_{\text{SOC,R}}$ around the
$z$-axis. We show the parameter $\Phi$, as extracted from our fits, for
selected heterostructures in Tab. 7 and Fig. 7 (c). For $\Theta=0\degree$ and
$\Theta=30\degree$ the reflection symmetry of the system demands
$\Phi=0\degree$ or $\Phi=180\degree$. Our DFT results confirm this.

However, there appears to be no distinct dependence of $\Phi$ on the twist
angle. What we can establish is that $\Phi$ depends strongly on the applied
transverse electric field, see Fig. 7, but is rather insensitive to the
lateral shift of the graphene and TMDC layers at a given twist angle, see Fig.
6.

Figure 7: The parameter $\Phi$, as obtained by fitting the spin-x and spin-y expectation values (for better readability, $\Phi$ is always plotted modulo 180°): (a) Spin-x and spin-y expectation values of the spin-up (red arrows) and spin-down (blue arrows) valence band of the graphene Dirac cone proximitized by WS2 with a twist angle of 0°. The k-path (dotted line) goes along a circular path around the K-point. Here the reflection symmetry along the $x$-axis demands $\Phi=0\degree$ or $\Phi=180\degree$. (b) shows the same as (a), but for a 5.2° supercell. (c) The angle $\Phi$ extracted from the fits for several heterostructures (see Tab. 7) plotted against the twist angle $|\theta|$. All results shown here are done with the electric field corrections. For 0° and 30°, we extract $\Phi=0$ and $\Phi=180\degree$ as expected. (d) Electric field dependence of $\Phi$ for four WS2/graphene heterostructures. Note that in c) and d) for the 5.2∘ and the 27∘ cases, we used supercells with $\theta<0$. Since the Rashba phase angle also changes sign, when changing the sign of the twist angle, the reported phase angles are misleading. We expect that for the given twist angles with a uniform sign, the Rashba phase angles will also have uniform sign. Material | twist angle $\theta$ [°] | E-field [V/nm] |  $\Phi$[°]
---|---|---|---
WS2 | 0 | 0 | 0
WS2 | 0 | 0.951 | 0
WS2 | -5.2 | 0 | -50
WS2 | -5.2 | 0.520 | -45
WS2 | -5.2 | 2.234 | -29
WS2 | 19.1 | 0 | 16
WS2 | 19.1 | -1.039 | 27
WS2 | 19.1 | 0.662 | 4
WS2 | -27 | 0 | -11
WS2 | -27 | 0.553 | -11
WS2 | -27 | 2.195 | -10
WS2 | 30 | 0 | 180
WS2 | 30 | -5.142 | 180
WS2 | 30 | -4.973 | 180
MoSe2 | 0 | -0.863 | 0
MoSe2 | -5.2 | 0.209 | -3
MoSe2 | 19.1 | -1.048 | 13
MoSe2 | -27 | 0.288 | -27
MoSe2 | 30 | -6.515 | 0
MoS2 | 19.1 | 0 | 41
MoS2 | 19.1 | 1.107 | 25
Table 7: Extracted parameter $\Phi$ for selected twisted graphene/TMDC
heterostructures.

##  Appendix E Effect of electric field

To correct the band offsets of our DFT calculations, we use electric fields in
the $z$-direction. This is important to obtain the right spin-orbit
parameters, since the energetic vicinity of the Dirac cone to the valence or
the conductance band of the TMDC can massively influence the effectivity of
the proximity effects. In Fig. 8 we depict this effect by exposing the 19.1°
MoS2 supercell to different electric fields. Without any electric field, the
Dirac cone of this supercell lies very close to the TMDC conductance band and
therefore has very strong proximity SOC. By applying a electric field, we can
bring the Dirac cone down in energy and exponentially reduce both the Rashba
spin orbit-coupling $\lambda_{\text{R}}$ and the valley-Zeeman SOC
$\lambda_{\text{VZ}}$. In Fig. 9 we show how Fig. 5 would look without the
electric field corrections.

Figure 8: The effect of the correcting electric field on the spin-orbit
parameters for the MoS2 19.1° supercell. The electric field increases the
conduction band offset $E_{C}$ and therefore reduces the proximity induced
SOC. Figure 9: The uncorrected spin-orbit parameters $\lambda_{\text{R}}$ (red
squares) and $\lambda_{\text{VZ}}$ (blue dots) of the graphene Dirac cone as
function of the twist angle $\Theta$ for the four TMDCs MoS2, MoSe2, WSe2 and
WS2 similar to Fig. 5. The transparent data points come from supercells with
built-in strain $\epsilon>10\%$ and are therefore less reliable. We removed
all data points with negative band offsets. The dotted line is a guide for the
eyes connecting reliable data points.

## References

  * Gmitra  _et al._ (2009) M. Gmitra, S. Konschuh, C. Ertler, C. Ambrosch-Draxl, and J. Fabian, “Band-structure topologies of graphene: Spin-orbit coupling effects from first principles,” Phys. Rev. B 80, 235431 (2009).
  * Sichau  _et al._ (2019) J. Sichau, M. Prada, T. Anlauf, T. J. Lyon, B. Bosnjak, L. Tiemann, and R. H. Blick, “Resonance microwave measurements of an intrinsic spin-orbit coupling gap in graphene: A possible indication of a topological state,” Phys. Rev. Lett. 122, 046403 (2019).
  * Han  _et al._ (2014) Wei Han, Roland K Kawakami, Martin Gmitra, and Jaroslav Fabian, “Graphene spintronics,” Nat. Nanotechnol. 9, 794 (2014).
  * Sierra  _et al._ (2021) Juan F. Sierra, Jaroslav Fabian, Roland K. Kawakami, Stephan Roche, and Sergio O. Valenzuela, “Van der waals heterostructures for spintronics and opto-spintronics,” Nature Nanotechnology 16, 856–868 (2021).
  * Žutić  _et al._ (2004) Igor Žutić, Jaroslav Fabian, and S. Das Sarma, “Spintronics: Fundamentals and applications,” Review of Modern Physics 76, 323 (2004).
  * Avsar  _et al._ (2020) A. Avsar, H. Ochoa, F. Guinea, B. Özyilmaz, B. J. van Wees, and I. J. Vera-Marun, “Colloquium: Spintronics in graphene and other two-dimensional materials,” Rev. Mod. Phys. 92, 021003 (2020).
  * Kane and Mele (2005a) C. L. Kane and E. J. Mele, “Quantum Spin hall effect in graphene,” Phys. Rev. Lett. 95, 226801 (2005a).
  * Kane and Mele (2005b) C. L. Kane and E. J. Mele, “Z2 topological order and the quantum spin hall effect,” Phys. Rev. Lett. 95, 146802 (2005b).
  * Qiao  _et al._ (2010) Zhenhua Qiao, Shengyuan A Yang, Wanxiang Feng, Wang-Kong Tse, Jun Ding, Yugui Yao, Jian Wang, and Qian Niu, “Quantum Anomalous Hall Effect in Graphene from Rashba and Exchange Effects,” Phys. Rev. B 82, 161414(R) (2010).
  * Ren  _et al._ (2016) Yafei Ren, Zhenhua Qiao, and Qian Niu, “Topological phases in two-dimensional materials: a review,” Rep. Prog. Phys. 79, 066501 (2016).
  * Frank  _et al._ (2018) Tobias Frank, Petra Högl, Martin Gmitra, Denis Kochan, and Jaroslav Fabian, “Protected pseudohelical edge states in ${\mathbb{z}}_{2}$-trivial proximitized graphene,” Phys. Rev. Lett. 120, 156402 (2018).
  * Högl  _et al._ (2020) Petra Högl, Tobias Frank, Klaus Zollner, Denis Kochan, Martin Gmitra, and Jaroslav Fabian, “Quantum anomalous hall effects in graphene from proximity-induced uniform and staggered spin-orbit and exchange coupling,” Phys. Rev. Lett. 124, 136403 (2020).
  * Gmitra and Fabian (2015) Martin Gmitra and Jaroslav Fabian, “Graphene on transition-metal dichalcogenides: A platform for proximity spin-orbit physics and optospintronics,” Phys. Rev. B 92, 155403 (2015).
  * Wang  _et al._ (2015a) Zhe Wang, Dong-Keun Ki, Hua Chen, Helmuth Berger, Allan H MacDonald, and Alberto F Morpurgo, “Strong interface-induced spin-orbit interaction in graphene on WS2,” Nat. Commun. 6, 8339 (2015a).
  * Gmitra  _et al._ (2016) 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, 155104 (2016).
  * Song  _et al._ (2018) Kenan Song, David Soriano, Aron W. Cummings, Roberto Robles, Pablo Ordejón, and Stephan Roche, “Spin Proximity Effects in Graphene/Topological Insulator Heterostructures,” Nano Lett. 18, 2033 (2018).
  * Zollner and Fabian (2019) Klaus Zollner and Jaroslav Fabian, “Single and bilayer graphene on the topological insulator ${\mathrm{bi}}_{2}{\mathrm{se}}_{3}$: Electronic and spin-orbit properties from first principles,” Phys. Rev. B 100, 165141 (2019).
  * Zollner and Fabian (2021a) Klaus Zollner and Jaroslav Fabian, “Heterostructures of graphene and topological insulators bi2se3, bi2te3, and sb2te3,” physica status solidi (b) 258, 2000081 (2021a).
  * Khokhriakov  _et al._ (2018) Dmitrii Khokhriakov, Aron W. Cummings, Kenan Song, Marc Vila, Bogdan Karpiak, André Dankert, Stephan Roche, and Saroj P. Dash, “Tailoring emergent spin phenomena in Dirac material heterostructures,” Sci. Adv. 4, eaat9349 (2018).
  * Avsar  _et al._ (2014) A Avsar, J Y Tan, T Taychatanapat, J Balakrishnan, G. K. W. Koon, Y Yeo, J Lahiri, A Carvalho, A S Rodin, E. C. T. O’Farrell, G Eda, A H Castro Neto, and B. Özyilmaz, “Spin–orbit proximity effect in graphene,” Nat. Commun. 5, 4875 (2014).
  * Yang  _et al._ (2016) Bowen Yang, Min-Feng Tu, Jeongwoo Kim, Yong Wu, Hui Wang, Jason Alicea, Ruqian Wu, Marc Bockrath, and Jing Shi, “Tunable spin–orbit coupling and symmetry-protected edge states in graphene/WS 2,” 2D Mater. 3, 031012 (2016).
  * Wang  _et al._ (2016) Zhe Wang, Dong-Keun Ki, Jun Yong Khoo, Diego Mauro, Helmuth Berger, Leonid S. Levitov, and Alberto F. Morpurgo, “Origin and Magnitude of ’Designer’ Spin-Orbit Interaction in Graphene on Semiconducting Transition Metal Dichalcogenides,” Phys. Rev. X 6, 041020 (2016).
  * Völkl  _et al._ (2017) Tobias Völkl, Tobias Rockinger, Martin Drienovsky, Kenji Watanabe, Takashi Taniguchi, Dieter Weiss, and Jonathan Eroms, “Magnetotransport in heterostructures of transition metal dichalcogenides and graphene,” Phys. Rev. B 96, 125405 (2017).
  * Zihlmann  _et al._ (2018) Simon Zihlmann, Aron W. Cummings, Jose H. Garcia, Máté Kedves, Kenji Watanabe, Takashi Taniguchi, Christian Schönenberger, and Péter Makk, “Large spin relaxation anisotropy and valley-zeeman spin-orbit coupling in ${\mathrm{wse}}_{2}$/graphene/$h$-bn heterostructures,” Phys. Rev. B 97, 075434 (2018).
  * Garcia  _et al._ (2018) Jose H. Garcia, Marc Vila, Aron W. Cummings, and Stephan Roche, “Spin transport in graphene/transition metal dichalcogenide heterostructures,” Chem. Soc. Rev. 47, 3359 (2018).
  * Safeer  _et al._ (2019) CK Safeer, Josep Ingla-Aynés, Franz Herling, José H Garcia, Marc Vila, Nerea Ontoso, M Reyes Calvo, Stephan Roche, Luis E Hueso, and Fèlix Casanova, “Room-temperature spin Hall effect in graphene/MoS2 van der Waals heterostructures,” Nano Lett. 19, 1074 (2019).
  * Herling  _et al._ (2020) Franz Herling, C. K. Safeer, Josep Ingla-Aynés, Nerea Ontoso, Luis E. Hueso, and Fèlix Casanova, “Gate tunability of highly efficient spin-to-charge conversion by spin hall effect in graphene proximitized with wse2,” APL Materials 8, 071103 (2020).
  * Khoo  _et al._ (2017) Jun Yong Khoo, Alberto F. Morpurgo, and Leonid Levitov, “On-Demand Spin-Orbit Interaction from Which-Layer Tunability in Bilayer Graphene,” Nano Lett. 17, 7003 (2017).
  * Omar and van Wees (2018) S. Omar and B. J. van Wees, “Spin transport in high-mobility graphene on ${\mathrm{ws}}_{2}$ substrate with electric-field tunable proximity spin-orbit interaction,” Phys. Rev. B 97, 045414 (2018).
  * Omar and van Wees (2017) S. Omar and B. J. van Wees, “Graphene-${\mathrm{ws}}_{2}$ heterostructures for tunable spin injection and spin transport,” Phys. Rev. B 95, 081404 (2017).
  * Fülöp  _et al._ (2021) Bálint Fülöp, Albin Márffy, Simon Zihlmann, Martin Gmitra, Endre Tóvári, Bálint Szentpéteri, Máté Kedves, Kenji Watanabe, Takashi Taniguchi, Jaroslav Fabian, Christian Schönenberger, Péter Makk, and Szabolcs Csonka, “Boosting proximity spin–orbit coupling in graphene/wse2 heterostructures via hydrostatic pressure,” npj 2D Materials and Applications 5, 82 (2021).
  * Gmitra and Fabian (2017) Martin Gmitra and Jaroslav Fabian, “Proximity effects in bilayer graphene on monolayer ${\mathrm{wse}}_{2}$: Field-effect spin valley locking, spin-orbit valve, and spin transistor,” Phys. Rev. Lett. 119, 146401 (2017).
  * Island  _et al._ (2019) J. O. Island, X. Cui, C. Lewandowski, J. Y. Khoo, E. M. Spanton, H. Zhou, D. Rhodes, J. C. Hone, T. Taniguchi, K. Watanabe, L. S. Levitov, M. P. Zaletel, and A. F. Young, “Spin–orbit-driven band inversion in bilayer graphene by the van der Waals proximity effect,” Nature 571, 85 (2019).
  * Tiwari  _et al._ (2021) Priya Tiwari, Saurabh Kumar Srivastav, and Aveek Bid, “Electric-field-tunable valley zeeman effect in bilayer graphene heterostructures: Realization of the spin-orbit valve effect,” Phys. Rev. Lett. 126, 096801 (2021).
  * Amann  _et al._ (2021) Julia Amann, Tobias Völkl, Denis Kochan, Kenji Watanabe, Takashi Taniguchi, Jaroslav Fabian, Dieter Weiss, and Jonathan Eroms, “Gate-tunable spin-orbit-coupling in bilayer graphene-wse2-heterostructures,” arXiv:2012.05718 (2021).
  * Cummings  _et al._ (2017) Aron W. Cummings, Jose H. Garcia, Jaroslav Fabian, and Stephan Roche, “Giant spin lifetime anisotropy in graphene induced by proximity effects,” Phys. Rev. Lett. 119, 206601 (2017).
  * Ghiasi  _et al._ (2017) Talieh S. Ghiasi, Josep Ingla-Aynés, Alexey A. Kaverzin, and Bart J. van Wees, “Large proximity-induced spin lifetime anisotropy in transition-metal dichalcogenide/graphene heterostructures,” Nano Letters 17, 7528–7532 (2017).
  * Benítez  _et al._ (2018) L. Antonio Benítez, Juan F. Sierra, Williams Savero Torres, Aloïs Arrighi, Frédéric Bonell, Marius V. Costache, and Sergio O. Valenzuela, “Strongly anisotropic spin relaxation in graphene-transition metal dichalcogenide heterostructures at room temperature,” Nat. Phys. 14, 303 (2018).
  * Cao  _et al._ (2018) Yuan Cao, Valla Fatemi, Shiang Fang, Kenji Watanabe, Takashi Taniguchi, Efthimios Kaxiras, and Pablo Jarillo-Herrero, “Unconventional superconductivity in magic-angle graphene superlattices,” Nature 556, 43–50 (2018).
  * Arora  _et al._ (2020) Harpreet Singh Arora, Robert Polski, Yiran Zhang, Alex Thomson, Youngjoon Choi, Hyunjin Kim, Zhong Lin, Ilham Zaky Wilson, Xiaodong Xu, Jiun-Haw Chu, _et al._ , “Superconductivity in metallic twisted bilayer graphene stabilized by wse 2,” Nature 583, 379–384 (2020).
  * Stepanov  _et al._ (2020) Petr Stepanov, Ipsita Das, Xiaobo Lu, Ali Fahimniya, Kenji Watanabe, Takashi Taniguchi, Frank HL Koppens, Johannes Lischner, Leonid Levitov, and Dmitri K Efetov, “Untying the insulating and superconducting orders in magic-angle graphene,” Nature 583, 375–378 (2020).
  * Balents  _et al._ (2020) Leon Balents, Cory R Dean, Dmitri K Efetov, and Andrea F Young, “Superconductivity and strong correlations in moiré flat bands,” Nature Physics 16, 725–733 (2020).
  * Carr  _et al._ (2017) Stephen Carr, Daniel Massatt, Shiang Fang, Paul Cazeaux, Mitchell Luskin, and Efthimios Kaxiras, “Twistronics: Manipulating the electronic properties of two-dimensional layered structures through their twist angle,” Phys. Rev. B 95, 075420 (2017).
  * Hennighausen and Kar (2021) Zachariah Hennighausen and Swastik Kar, “Twistronics: a turning point in 2d quantum materials,” Electronic Structure 3, 014004 (2021).
  * Ribeiro-Palau  _et al._ (2018) Rebeca Ribeiro-Palau, Changjian Zhang, Kenji Watanabe, Takashi Taniguchi, James Hone, and Cory R. Dean, “Twistable electronics with dynamically rotatable heterostructures,” Science 361, 690–693 (2018).
  * Carr  _et al._ (2020) Stephen Carr, Shiang Fang, and Efthimios Kaxiras, “Electronic-structure methods for twisted moiré layers,” Nature Reviews Materials 5, 748–763 (2020).
  * Zollner and Fabian (2021b) Klaus Zollner and Jaroslav Fabian, “Engineering proximity exchange by twisting: Reversal of ferromagnetic and emergence of antiferromagnetic dirac bands in graphene/cr2ge2te6,” arXiv:2108.03984 (2021b).
  * Yan  _et al._ (2021) Shiming Yan, Shengmei Qi, Dunhui Wang, and Wenbo Mi, “Novel electronic structures and magnetic properties in twisted two-dimensional graphene/janus 2h–vsete heterostructures,” Physica E: Low-dimensional Systems and Nanostructures 134, 114854 (2021).
  * David  _et al._ (2019) 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, 085412 (2019).
  * Li and Koshino (2019) Yang Li and Mikito Koshino, “Twist-angle dependence of the proximity spin-orbit coupling in graphene on transition-metal dichalcogenides,” Phys. Rev. B 99, 075438 (2019).
  * Pezo  _et al._ (2021) Armando Pezo, Zeila Zanolli, Nils Wittemeier, Pablo Ordejón, Adalberto Fazzio, Stephan Roche, and Jose H Garcia, “Manipulation of spin transport in graphene/transition metal dichalcogenide heterobilayers upon twisting,” 2D Materials 9, 015008 (2021).
  * Koshino (2015) Mikito Koshino, “Interlayer interaction in general incommensurate atomic layers,” New Journal of Physics 17, 015014 (2015).
  * Hou  _et al._ (2017) Siyao Hou, Lihong Han, Liyuan Wu, Ruge Quhe, and Pengfei Lu, “Robust quasi-ohmic contact against angle rotation in noble transition-metal-dichalcogenide/graphene heterobilayers,” RSC Adv. 7, 45896–45901 (2017).
  * Felice  _et al._ (2017) D Di Felice, E Abad, C González, A Smogunov, and Y J Dappe, “Angle dependence of the local electronic properties of the graphene/MoS2interface determined byab initiocalculations,” Journal of Physics D: Applied Physics 50, 17LT02 (2017).
  * Koda  _et al._ (2016) Daniel S Koda, Friedhelm Bechstedt, Marcelo Marques, and Lara K Teles, “Coincidence lattices of 2d crystals: heterostructure predictions and applications,” The Journal of Physical Chemistry C 120, 10895–10908 (2016).
  * Wang  _et al._ (2015b) Zilu Wang, Qian Chen, and Jinlan Wang, “Electronic structure of twisted bilayers of graphene/mos2 and mos2/mos2,” The Journal of Physical Chemistry C 119, 4752–4758 (2015b), https://doi.org/10.1021/jp507751p .
  * Zollner  _et al._ (2019) Klaus Zollner, Paulo E. Faria Junior, and Jaroslav Fabian, “Strain-tunable orbital, spin-orbit, and optical properties of monolayer transition-metal dichalcogenides,” Phys. Rev. B 100, 195126 (2019).
  * Si  _et al._ (2016) Chen Si, Zhimei Sun, and Feng Liu, “Strain engineering of graphene: a review,” Nanoscale 8, 3207–3217 (2016).
  * Choi  _et al._ (2010) Seon-Myeong Choi, Seung-Hoon Jhi, and Young-Woo Son, “Effects of strain on electronic properties of graphene,” Phys. Rev. B 81, 081407 (2010).
  * Bahn and Jacobsen (2002) S. R. Bahn and K. W. Jacobsen, “An object-oriented scripting interface to a legacy electronic structure code,” Comput. Sci. Eng. 4, 56 (2002).
  * Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, “Inhomogeneous electron gas,” Phys. Rev. 136, B864 (1964).
  * Giannozzi and et al. (2009) Paolo Giannozzi and et al., “Quantum espresso: a modular and open-source software project for quantum simulations of materials,” J. Phys.: Cond. Mat. 21, 395502 (2009).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, “From ultrasoft pseudopotentials to the projector augmented-wave method,” Phys. Rev. B 59, 1758 (1999).
  * Perdew  _et al._ (1996) John P. Perdew, Kieron Burke, and Matthias Ernzerhof, “Generalized gradient approximation made simple,” Phys. Rev. Lett. 77, 3865 (1996).
  * Grimme (2006) Stefan Grimme, “Semiempirical gga-type density functional constructed with a long-range dispersion correction,” J. Comput. Chem. 27, 1787 (2006).
  * Grimme  _et al._ (2010) Stefan Grimme, Jens Antony, Stephan Ehrlich, and Helge Krieg, “A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu,” J. Chem. Phys. 132, 154104 (2010).
  * Barone  _et al._ (2009) Vincenzo Barone, Maurizio Casarin, Daniel Forrer, Michele Pavone, Mauro Sambi, and Andrea Vittadini, “Role and effective treatment of dispersive forces in materials: Polyethylene and graphite crystals as test cases,” J. Comput. Chem. 30, 934 (2009).