# Twisted bilayer graphene reveals its flat bands under spin pumping

Sonia Haddad1,2<EMAIL_ADDRESS>Takeo Kato2 Lassaad Mandhour1
1Laboratoire de Physique de la Matière Condensée, Faculté des Sciences de
Tunis, Université Tunis El Manar, Campus Universitaire 1060 Tunis, Tunisia
2 Institute for Solid State Physics, University of Tokyo, Kashiwa, Chiba
277-8581, Japan

(August 27, 2024)

###### Abstract

The salient property of the electronic band structure of twisted bilayer
graphene (TBG), at the so-called magic angle (MA), is the emergence of flat
bands around the charge neutrality point. These bands are associated with the
observed superconducting phases and the correlated insulating states occurring
at the MA. Scanning tunneling microscopy combined with angle resolved
photoemission spectroscopy are usually used to visualize the flatness of the
band structure of the TGB at the MA. Here, we theoretically argue that spin
pumping (SP) provides a direct probe of the flat bands of TBG and an accurate
determination of the MA. We consider a junction separating a ferromagnetic
insulator and an heterostructure of TBG adjacent to a monolayer of a
transition metal dichalcogenide. We show that the Gilbert damping of the
ferromagnetic resonance experiment through this junction depends on the twist
angle of TBG, and shows a sharp drop at the MA. Our results open the way to a
twist switchable spintronics in twisted van der Waals heterostructures.

Introduction. –  Stacking two graphene layers with a relative twist angle
$\theta$ results into a moiré superstructure which is found to host, in the
vicinity of the so-called magic angle (MA) $\theta_{M}\sim 1.1^{\circ}$,
unconventional superconductivity and strongly correlated insulating states
Herrero1 ; Herrero2 ; Yank . There is a general consensus that such strong
electronic correlations originate from the moiré flat bands emerging at the MA
around the charge neutrality point Volovik ; Senthil ; Wu ; Roy ; Bernevig ;
Efetov ; Young ; Herrero3 . The tantalizing signature of the flat bands have
been experimentally demonstrated by probing the corresponding peaks of the
density of states using transport Herrero2 ; Pablo18 ; Dean ; Efetov19 ;
Dean19 , electronic compressibility measurements Pablo19 ; Pablo20 , scanning
tunneling microscopy (STM) and spectroscopy (STS) Eva ; Kerelsky ; Yazdani19 ;
Nadj19 ; Yazdani20 ; Nadj21 ; Eva2 ; Yazdani22 . The direct evidence of these
flat bands has been reported by angle resolved photoemission spectroscopy
(ARPES) measurements combined to different imaging techniques Utama ; Efetov21
; Sato . However, spectroscopic measurements on magic-angle TBG raise many
technical challenges related to the need of an accurate control of the twist
angle, and the necessity to have non-encapsulated samples which can degrade in
air. Efetov21 .

Here we propose a non-invasive method to probe the flat bands of the TBG and
accurately determine the MA. This method is based on spin pumping (SP) induced
by ferromagnetic resonance (FMR) Bauer1 ; Bauer2 ; SP-book ; Hellman , where
the increase in the FMR linewidth, given by the Gilbert damping (GD)
coefficient, provides insight into the spin excitations of the nonmagnetic
(NM) material adjacent to the ferromagnet Qiu2016 ; Yang2018 ; Han2020 . SP is
expected to be efficient if the NM has high spin–orbit coupling (SOC) strength
Hait .
In our work, we consider spin injection from a ferromagnetic insulator (FI)
into a TBG aligned on a monolayer of a transition metal dichalcogenides (TMD)
which are considered as good substrate candidates to induce relatively strong
SOC in graphene and TBG Castro ; Morpu15 ; Morpu16 ; Bock16 ; Casa ; Shi ;
Wees ; Eroms ; Eroms2 ; Makk ; BouchitaPRL ; BLG ; Omar ; Valen ; Roche ;
Zaletel ; Wang ; Bouchiat19 ; David ; Zaletel20 ; Bouchiat21 ; Lin ; Alex ;
Bhowmick .

We theoretically study a planar junction of a FI and a TBG adjacent to
$\mathrm{WSe}_{2}$ (TBG/$\mathrm{WSe}_{2}$) as depicted in figure 1. We
consider the case where a microwave of a frequency $\Omega$ is applied to this
junction, and focus on the role of the twist angle on the FMR linewidth of the
FI. We start by deriving the the continuum model of the heterostructure
TBG/$\mathrm{WSe}_{2}$. Using a second-order perturbative theory, with respect
to the interfacial exchange interaction Ohnuma ; Matsuo18 ; Kato19 ; Kato20 ;
Matsuo-JP ; Matsuo20 ; Yama , we determine the twist angle dependence of the
correction to the GD coefficient.

Figure 1: Schematic representation of the junction between a ferromagnetic
insulator (FI) and an heterostructure of TBG adjacent to a monolayer of
$\mathrm{WSe}_{2}$. The labels (1) and (2) denote the graphene layers of TBG
represented by the red and the blue lines. The red arrow indicates the spin
orientation of the FI characterized by an average spin, written in the
coordinate frame of the FI magnetization as
$\langle\mathbf{S}_{FI}\rangle=\left(S_{0},0,0\right)$.

Continuum model. –  In TBG with a twist angle $\theta$, the Hamiltonian of a
graphene layer $l$ ($l=1,2$), rotated at an angle $\theta_{l}$, is

$\displaystyle
h_{l,\text{rot}}(\mathbf{k})=e^{i\frac{\theta_{l}}{2}\sigma_{z}}h_{l}(\mathbf{k})e^{-i\frac{\theta_{l}}{2}\sigma_{z}}$
(1)

where $\theta_{2}=-\theta_{1}=\frac{\theta}{2}$, $h_{l}(\mathbf{k})$ is the
monolayer Hamiltonian. $h_{1}(\mathbf{k})$ reduces to the graphene Hamiltonian
given, in the continuum limit, by $h_{1}(\mathbf{k})=-\hbar
v_{F}\mathbf{k}\cdot\mathbf{\sigma}^{\ast}$, where $v_{F}$ is the Fermi
velocity, $\xi$ is the valley index,
$\mathbf{\sigma}^{\ast}=\left(\xi\sigma_{x},\sigma_{y}\right)$ and
$\sigma_{i}$ ($i=x,y,z$) are the sublattice-Pauli matrices. The layer (2), in
contact with $\mathrm{WSe}_{2}$ monolayer, is descried by the Hamiltonian Alex

$\displaystyle
h_{2}(\mathbf{k})=h_{1}(\mathbf{k})+h_{\text{SOC}}+\frac{m}{2}\sigma_{z}$ (2)

where $h_{\text{SOC}}$ is given by

$\displaystyle h_{\text{SOC}}=\frac{\lambda_{I}}{2}\xi
s_{z}+\frac{\lambda_{R}}{2}\left(\xi\sigma_{x}s_{y}-\sigma_{y}s_{x}\right)+\frac{\lambda_{KM}}{2}\xi\sigma_{z}s_{z},$
(3)

$s_{i}$ ($i=x,y,z$) are the spin-Pauli matrices, $\lambda_{I}$, $\lambda_{R}$
and $\lambda_{KM}$ correspond, respectively, to the Ising, Rashba and Kane-
Mele SOC parameters Alex . The variation ranges of these parameters are
$\lambda_{I}\sim 1-5\;\mathrm{meV}$, $\lambda_{R}\sim 1-15\;\mathrm{meV}$,
while $\lambda_{KM}$ is expected to be small SOC1 ; SOC2 ; SOC3 ; SOC4 ; SOC5
; SOC6 ; SOC7 . The last term in Eq. 2 is due to the inversion symmetry
breaking induced by the TMD layer. Hereafter, we neglect this term regarding
the small value of $m$ compared to the SOC parameters Alex .
As in the case of TBG Mc11 , the low-energy Hamiltonian of
TBG/$\mathrm{WSe}_{2}$ reduces, at the valley $\xi$, to

$\displaystyle
H_{\xi,SOC}(\mathbf{k})=\begin{pmatrix}h_{1}(\mathbf{k})&T_{1}&T_{2}&T_{3}\\\
T^{\dagger}_{1}&h_{2,1}(\mathbf{k})&0&0\\\
T^{\dagger}_{2}&0&h_{2,2}(\mathbf{k})&0\\\
T^{\dagger}_{3}&0&0&h_{2,3}(\mathbf{k})\\\ \end{pmatrix}.$ (4)

$H_{\xi,\text{SOC}}(\mathbf{k})$ is written in the basis
$\Psi=\left(\psi_{0}(\mathbf{k}),\psi_{1}(\mathbf{k}),\psi_{2}(\mathbf{k}),\psi_{3}(\mathbf{k})\right)$
constructed on the four-component spin-sublattice spinor
$\psi_{0}(\mathbf{k})$ and $\psi_{j}(\mathbf{k})$, ($j=1,2,3$) corresponding,
respectively, to layer $(1)$ and layer $(2)$ (see sections I and II of the
supplemental material supp ). The momentum $\mathbf{k}$ is measured relatively
to the Dirac point $\mathbf{K}_{1\xi}$ of layer (1). In Eq. 4,
$h_{2,j}\left(\mathbf{k}\right)=h_{2}\left(\mathbf{k}+\mathbf{q}_{j\xi}\right)$,
($j=1,2,3$) where $\mathbf{q}_{j\xi}$ are the vectors connecting
$\mathbf{K}_{1\xi}$ to its three neighboring Dirac points $\mathbf{K}_{2\xi}$
of layer (2) in the moiré Brillouin zone (mBZ) Mc11 , and are given by
$\mathbf{q}_{1\xi}=\mathbf{K}_{1\xi}-\mathbf{K}_{2\xi}$,
$\mathbf{q}_{2\xi}=\mathbf{q}_{1\xi}+\xi\mathbf{G}^{M}_{1}$,
$\mathbf{q}_{3\xi}=\mathbf{q}_{1\xi}+\xi\left(\mathbf{G}^{M}_{1}+\mathbf{G}^{M}_{2}\right)$,
where $\left(\mathbf{G}^{M}_{1},\mathbf{G}^{M}_{2}\right)$ is the mBZ basis
(see section I of the supplemental material supp ).
We consider AA stacked moiré domains of TBG, where the flat bands are expected
to emerge. The $T_{j}$ are the spin-independent interlayer coupling matrices
given, in the unrelaxed TBG, by
$T_{1}=w\left(\mathbb{I}_{\sigma}+\sigma_{x}\right)$,
$T_{2}=w\left(\mathbb{I}_{\sigma}-\frac{1}{2}\sigma_{x}+\xi\frac{\sqrt{3}}{2}\sigma_{y}\right)$
and
$T_{3}=w\left(\mathbb{I}_{\sigma}-\frac{1}{2}\sigma_{x}-\xi\frac{\sqrt{3}}{2}\sigma_{y}\right)$
supp , where $w\sim 118\,\mathrm{meV}$ w  is the interlayer tunneling
amplitude and $\mathbb{I}_{\sigma}$ is the identity matrix acting on the
sublattice indices.

Using the perturbative approach of Ref. Mc11 , we derive, from Eq. 4, the
effective low-energy Hamiltonian $H_{\xi,\text{SOC}}^{(1)}(\mathbf{k})$ of
TBG/$\mathrm{WSe_{2}}$ (see section II of the supplemental material supp ). To
the leading order in $\mathbf{k}$, $H_{\xi,\text{SOC}}^{(1)}(\mathbf{k})$
reads as supp

$\displaystyle
H_{\xi,\text{SOC}}^{(1)}(\mathbf{k})=\frac{\langle\Psi|H_{\xi,\text{SOC}}|\Psi\rangle}{\langle\Psi|\Psi\rangle}=\psi^{\dagger}_{0}\left[h_{\text{eff}}\left(\mathbf{k}\right)+h_{\text{eff}}^{\text{SOC}}\right]\psi_{0},$
(5) $\displaystyle h_{\text{eff}}\left(\mathbf{k}\right)=-\frac{\hbar
v}{\langle\Psi|\Psi\rangle}\left\\{k_{x}\left[\left(1-3\alpha^{2}\right)\xi\sigma_{x}\mathbb{I}_{s}-\frac{3\alpha^{2}}{\hbar
vq_{0}}\left(\xi\lambda_{I}\sigma_{y}s_{z}+\lambda_{R}\left(\xi\sigma_{y}s_{y}-\sigma_{x}s_{x}\right)\right)\right]\right.$
$\displaystyle\hskip
79.6678pt\left.+k_{y}\left[\left(1-3\alpha^{2}\right)\sigma_{y}\mathbb{I}_{s}-\frac{3\alpha^{2}}{\hbar
vq_{0}}\left(-\lambda_{I}\sigma_{x}s_{z}+\lambda_{R}\left(\sigma_{x}s_{y}+\xi\sigma_{y}s_{x}\right)\right)\right]\right\\},$
(6) $\displaystyle
h_{\text{eff}}^{\text{SOC}}=\frac{3\alpha^{2}}{\langle\Psi|\Psi\rangle}\left[\xi\lambda_{I}s_{z}\mathbb{I}_{\sigma}+\frac{\lambda_{R}}{2}\left(s_{x}\sigma_{y}-\xi
s_{y}\sigma_{x}\right)\right],$ (7)

$\langle\Psi|\Psi\rangle\sim 1+6\alpha^{2}$, $\alpha=\frac{w}{\hbar
v_{F}q_{0}}$, $q_{0}=|\mathbf{q}_{j\xi}|=\frac{4\pi}{3a}\theta$, $a$ is the
graphene lattice constant and $\sigma_{i}$, ($i=x,y,z$) act now on the band
indices $\sigma=\pm$ of the eigenergies of $H_{\xi,\text{SOC}}^{(1)}$ denoted
$E_{\sigma,\pm}$ and given, to the leading orders in $\mathbf{k}$ and
$\frac{\lambda_{I,R}}{\hbar v_{F}q_{0}}$, by

$\displaystyle
E(\mathbf{k}))_{\sigma,\pm}=\frac{\sigma}{\langle\Psi|\Psi\rangle}\sqrt{f_{1}(\mathbf{k})\pm
6\alpha^{2}\sqrt{f_{2}(\mathbf{k})}}$ (8) $\displaystyle
f_{1}(\mathbf{k})=(\hbar
v)^{2}\left(1-3\alpha^{2}\right)^{2}||{\mathbf{k}}||^{2}+\frac{9}{2}\alpha^{4}\left(2\lambda_{I}^{2}+\lambda_{R}^{2}\right)$
$\displaystyle f_{2}(\mathbf{k})=(\hbar
v)^{2}\left(1-3\alpha^{2}\right)^{2}||{\mathbf{k}}||^{2}\left(\lambda_{I}^{2}+\frac{1}{4}\lambda_{R}^{2}\right)+\frac{9}{16}\alpha^{4}\lambda_{R}^{4}.$

Equation 7 shows that the SOC parameters $\lambda_{I}$ and $\lambda_{R}$ are
renormalized by the moiré structure of TBG to

$\displaystyle\tilde{\lambda}_{I}=\frac{6\alpha^{2}}{\langle\Psi|\Psi\rangle}\lambda_{I},\;\tilde{\lambda}_{R}=\frac{3\alpha^{2}}{\langle\Psi|\Psi\rangle}\lambda_{R}.$
(9)

The expression of $H_{\xi,\text{SOC}}^{(1)}$ (Eq. 5) can be taken as a
starting point to unveil the role of SOC in the emergence of the stable
superconducting phase observed, at $\theta\sim 0.8^{\circ}$, in TBG adjacent
to $\mathrm{WSe}_{2}$ Alex .

To derive the GD, one needs to determine the magnon self-energy of the FI,
resulting from the interfacial exchange interaction between the FI and the
heterostructure TBG/$\mathrm{WSe}_{2}$.

Gilbert damping. –  In the absence of a junction, the magnon Green function of
the FI is defined as Yama

$\displaystyle
G_{0}\left(\mathbf{q},i\omega_{n}\right)=\frac{2S_{0}/\hbar}{i\omega_{n}-\omega_{\mathbf{q}}+i\alpha_{G}\omega_{n}}$
(10)

where $\omega_{n}=2\pi n/\hbar\beta$ are the Matsubara frequencies for bosons,
$S_{0}$ is the amplitude of the average spin per site, $\alpha_{G}\sim
10^{-4}-10^{-3}$ is the GD strength. In the presence of the interfacial
coupling, the magnon Green function is given by the Dyson equation (see
section IV of the supplemental material)

$\displaystyle
G_{0}\left(\mathbf{q},i\omega_{n}\right)=\frac{1}{G_{0}(\mathbf{q},i\omega_{n})^{-1}-\Sigma(\mathbf{q},i\omega_{n})},$
(11)

where $\Sigma\left(\mathbf{q},i\omega_{n}\right)$ is the magnon self-energy
resulting from the interfacial exchange interaction. For the sake of
simplicity, we neglect the real part of
$\Sigma\left(\mathbf{q},i\omega_{n}\right)$ related to the shift of the FMR.
We also assume that the interface is clean which turns out to consider that,
in the FMR experiment, the microwave irradiation induces a uniform spin
precession, which limits the magnon self-energy to the processes with
$\mathbf{q}=0$ Funato .

In the second order perturbation with respect to the interfacial exchange
interaction $T_{\mathbf{q}}$, the magnon self-energy (Eq. 11) is written as
Yama

$\displaystyle\Sigma(\mathbf{q},i\omega_{n})$ $\displaystyle=$
$\displaystyle\frac{|T_{\mathbf{q}}|^{2}}{4\beta}\sum_{\mathbf{k},i\omega_{m}}\mathrm{Tr}\left[\sigma_{s}^{x^{\prime},-}\;\hat{g}(\mathbf{k},i\omega_{m})\right.$
(12) $\displaystyle\times$
$\displaystyle\left.\sigma_{s}^{x^{\prime},+}\;\hat{g}(\mathbf{k}+\mathbf{q},i\omega_{m}+i\omega_{n})\right].$

$\sigma_{s}^{x^{\prime},\pm}$ are the electronic spin ladder operators written
in the coordinate system $\left(x^{\prime},y^{\prime},z^{\prime}\right)$ of
the FI magnetization characterized by an average spin
$\langle\mathbf{S}_{FI}\rangle=\left(S_{0},0,0\right)$ (see section IV of the
supplemental material supp ). $\hat{g}(\mathbf{k},i\omega_{m})$ is the
electronic Matsubara Green function given by

$\displaystyle\hat{g}(\mathbf{k},i\omega_{n})=\left[i\omega_{n}\mathbb{I}-H_{\text{SOC}}^{(1)}(\mathbf{k})\right]^{-1}$
(13)

where $\omega_{n}=(2n+1)\pi/\hbar\beta$ are the fermionic Matsubara
frequencies.

In the basis of the spin-band four-component spinor
$\Psi=\left(\psi_{+,\uparrow},\psi_{+,\downarrow},\psi_{-,\uparrow},\psi_{-,\downarrow}\right)$
Alex , $\hat{g}(\mathbf{k},i\omega_{n})$ reads as

$\displaystyle\hat{g}(\mathbf{k},i\omega_{n})=\hat{g}_{0}(\mathbf{k},i\omega_{n})\mathbb{I}_{s}+\mathbf{\hat{g}}(\mathbf{k},i\omega_{n})\cdot\mathbf{s}$
(14)

where $\mathbf{s}=(s_{x},s_{y},s_{z})$ are the spin-Pauli matrices,
$\mathbf{\hat{g}}=(\hat{g}_{x},\hat{g}_{y},\hat{g}_{z})$, $\hat{g}_{0}$ and
$\hat{g}_{i}$ ($i=x,y,z$) are expressed, to the leading order in the SOC, in
terms of the band-Pauli matrices $\sigma_{i}$ (see section III of the
supplemental material supp ).

Using the analytical continuation $i\omega_{n}\rightarrow\omega+i\delta$ in
Eq. 11, we obtain the magnon retarded Green function

$\displaystyle
G_{0}^{R}\left(\mathbf{q},i\omega_{n}\right)=\frac{2S_{0}/\hbar}{\omega-\omega_{\mathbf{q=0}}+i\left(\alpha_{G}+\delta\alpha_{G}\right)\omega},$
(15)
$\displaystyle\delta\alpha_{G}\left(\omega\right)\equiv-\frac{2S_{0}}{\hbar\omega}\mathrm{Im}\Sigma^{R}\left(\mathbf{q=0},\omega\right)$
(16)

$\delta\alpha_{G}$ is the correction to the GD due to the interfacial exchange
interaction between the FI and the TBG/$\mathrm{WSe_{2}}$. In general,
$\delta\alpha_{G}$ depends on $\omega$. However, since the ferromagnetic peak,
given by $\mathrm{Im}G_{0}^{R}$, is sharp enough, namely
$\alpha_{G}+\delta\alpha_{G}\ll 1$, one can replace $\omega$ by the FMR
frequency $\Omega$

$\displaystyle\delta\alpha_{G}\sim-\frac{2S_{0}}{\hbar\Omega}\mathrm{Im}\Sigma^{R}\left(\mathbf{q=0},\Omega\right)$
(17)

Carrying out the summation over $\omega_{m}$ in Eq. 12, we obtain analytical
expressions of the magnon self-energy (see section IV of the supplemental
material supp ). The sum over the electronic states
$\mathbf{k}=\left(k,\varphi_{\mathbf{k}}\right)$ runs over the states included
within a cutoff, $k_{c}=q_{0}/2$, on the momentum amplitude $k$, where the
Dirac cone description of the monolayer Hamiltonian is expected to hold.

In the Following, we discuss the behavior of the normalized GD coefficient

$\displaystyle\delta\alpha_{G}/\alpha_{G}^{0}=\left(\frac{\lambda}{\hbar\Omega}\right)^{2}\tilde{\Sigma}\left(\mathbf{q}=\mathbf{0},\Omega\right),$
(18)

where $\delta\alpha_{G}$ is given by Eq. 17, $\tilde{\Sigma}$ is a
dimensionless function depending on the twist angle $\theta$, temperature $T$,
the chemical potential $\mu$ and the orientation of the FI magnetization,
$\alpha_{G}^{0}=2S_{0}\left(\frac{|T_{\mathbf{0}}|}{\lambda}\right)^{2}$ and
$\lambda=\frac{\lambda_{I}+\lambda_{R}}{2}$ is the average SOC (for details,
see section IV of the supplemental material supp ).

Discussion. –  In figure 2, we plot $\delta\alpha_{G}/\alpha_{G}^{0}$ (Eq.18),
as a function of the twist angle $\theta$, for the undoped TBG, at different
temperatures and for a fixed FMR energy $\hbar\Omega=0.06\,\mathrm{meV}$ which
corresponds to the yttrium iron garnet. The SOC parameters are
$\lambda_{I}=3\,\mathrm{meV}$ and $\lambda_{I}=4\,\mathrm{meV}$ as in
Ref.Alex, .

Figure 2: Normalized GD, $\delta\alpha_{G}/\alpha_{G}^{0}$ (Eq.18), as a
function of the twist angle at different temperature ranges. Calculations are
done for $\lambda_{I}=3\;\mathrm{meV}$, $\lambda_{R}=4\;\mathrm{meV}$, $\mu=0$
and for a FMR energy $\hbar\Omega=0.06\;\mathrm{meV}$.

Figure 2 shows that, regardless of the temperature range, $\delta\alpha_{G}$
increases by decreasing $\theta$ but drops sharply at the MA where it exhibits
a relatively small peak, which is smeared out at low temperature.

Putting aside its drop at the MA, the enhancement of $\delta\alpha_{G}$, by
decreasing $\theta$, results from the dependence of the magnon self-energy
(Eq. 12) on the effective SOC (Eq. 9) which increase as $\frac{1}{\theta^{2}}$
by decreasing $\theta$ (see section II of the supplemental material supp ).

To understand the behavior of $\delta\alpha_{G}$ at the MA one needs to go
back to the band structure, $E_{\sigma,\pm}(\mathbf{k})$ (Eq. 8), of the
continuum Hamiltonian of TBG/$\mathrm{WSe_{2}}$, which is depicted in Fig. 3
at different twist angles. The arrows indicate the out-of-plane electronic
spin projection $\langle s_{z}\rangle$ which has been numerically calculated,
in Ref.Alex, , within the perturbative approach Mc11 , and which can be
deduced from Eq. 7.

$\begin{array}[]{cccc}\includegraphics[width=108.405pt]{energy-plot-
theta=0,5-MA}\includegraphics[width=108.405pt]{energy-plot-
MA-.eps}\includegraphics[width=108.405pt]{energy-plot-
MA+.eps}\includegraphics[width=108.405pt]{energy-plot-
theta=1,2-MA}\end{array}$

Figure 3: Band structure of TBG/$\mathrm{WSe_{2}}$ in the continuum limit (Eq.
8) at (a) $\theta=0.5^{\circ}$, (b) $\theta=\theta^{-}_{M}=1.043^{\circ}$, (c)
$\theta=\theta^{+}_{M}=1.058^{\circ}$ and $\theta=1.2^{\circ}$. The dashed
lines represent the bands at the MA ($\theta_{M}=1.05^{\circ}$). The red
(blue) arrows correspond to the out-of-plane electronic spin projection
$\langle s_{z}\rangle=+1$ ($\langle s_{z}\rangle=-1$) Alex .

Away from the MA, the band dispersion gets larger as $\theta$ decreases and,
in particular, the separation between bands with opposite $\langle
s_{z}\rangle$, involved in the SP process, increases. This behavior is due to
the angle dependence of the effective Fermi velocity of TBG $v^{\ast}\sim
v_{F}\frac{1-3\alpha^{2}}{1+6\alpha^{2}}$ (see sections I and II of the
supplemental Material supp ).
However, around the MA at $\theta^{-}_{M}$ and $\theta^{+}_{M}$, the bands
$E_{+,-}$ and $E_{-,-}$ get closer, compared to those at $\theta_{M}$.
The expression of the GD (Eq. 18) includes transitions between bands with
opposite $\langle s_{z}\rangle$ (see section IV of the supplemental material
supp ). These transitions depend on the statistical weight $\Delta
f(E)=f(E_{\langle s_{z}\rangle})-f(E_{-\langle s_{z}\rangle})$ where $f(x)$ is
the Fermi-Dirac function and $E_{\langle s_{z}\rangle}$ is the energy band
with a spin orientation $\langle s_{z}\rangle$.

Figure 4: Schematic representation of the band structure $E_{\sigma,\pm}$ (Eq.
8) and the Fermi-Dirac distribution $f(E)$. The bands in dashed and green
lines correspond, respectively, to the MA and to an angle $\theta$ far from
the MA. The red (blue) arrows represent the projection of the out-of-plane
spin projection $\langle S_{z}\rangle=+1$ ($\langle S_{z}\rangle=-1$). Around
the MA, the bands are almost flat and the statistical weights $\Delta f(E)$,
corresponding to the transitions between $E_{-,+}\rightarrow E_{+,+}$ and
$E_{-,-}\rightarrow E_{+,-}$, are small compared to the case of the twist
angle away from the MA, where the band dispersion is larger. As a consequence,
the correction to the Gilbert damping $\delta\alpha_{G}$, which depends on
$\Delta f(E)$, is reduced at the MA (for details, see section IV of the
supplemental material).

In Fig. 4, we plot a schematic representation of the band structure of the
continuum model (Eq.8) and the Fermi-Dirac distribution $f(E)$ at a given
temperature $T$. The dashed (green) line corresponds to a twist angle $\theta$
close to (far from) the MA (for details, see section IV of the supplemental
material). The band dispersion gets larger as $\theta$ moves away from the MA
(Fig. 3) and the separation between the bands with opposite $\langle
S_{z}\rangle$ increases. As a consequence, the corresponding statistical
weight $\Delta f(E)$ is enhanced compared to the case around the MA. This
behavior explains the drop of the GD at the MA.

For twist angles around the MA ($\theta^{+}_{M}$ and $\theta^{-}_{M}$), the
separating between the bands $E_{-,-}$ and $E_{+,-}$ is smaller compared to
the case at the MA (Fig. 3), which reduces the statistical weight $\Delta
f(E)$. This behavior gives rise to the small peak at the MA (Fig. 2), which
vanishes at low temperature ($k_{B}T<\lambda$) where bands around the MA have
the same statistical weight $\Delta f(E)=1$.

It comes out that the twist angle dependence of $\delta\alpha_{G}$ is a direct
probe of the emergence of the flat bands in the TBG. On the other hand, the
temperature dependence of the fine structure around the MA provides, not only,
an accurate measurement of the MA but also an estimation of the SOC induced in
TBG adjacent to a monolayer of TMD.

Conclusion. –  To conclude, we have proposed an experiment to probe the flat
bands of TBG and to measure accurately its MA. The experiment is based on a
spin pumping measurement through a junction separating a FI and a TBG adjacent
to a monolayer of $\mathrm{WSe_{2}}$. We first derived the continuum model of
TBG with SOC, which constitutes a first step to develop an analytical
understanding of the emergence of a stable superconducting state at small
twist angle observed in TBG in proximity to $\mathrm{WSe_{2}}$ Alex . Our
results show that, the twist angle dependence of the GD correction,
$\delta\alpha_{G}$ exhibits a drop at the MA with a temperature-dependent fine
structure. The latter provides an accurate determination of the MA and an
estimation of the SOC induced in TBG by its proximity to the TMD layer. Our
work opens the gate to a twist tunable spintronics in twisted layered
heterostructures.

Acknowledgment. –  We thank Prof. Mamoru Matsuo for stimulating discussions.
S. H. acknowledges the kind hospitality of ISSP where this work was carried
out. S. H. was supported by the ISSP visiting professor programme.

## References

  * (1) Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature, 556, 43 (2018).
  * (2) Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, Nature, 556, 80 (2018).
  * (3) M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, C. R. Dean, Science 363, eaav1910 (2019).
  * (4) N. B. Kopnin, T. T. Heikkila, and G. E. Volovik, Phys. Rev. B 83, 220503(R) (2011).
  * (5) H. C. Po, L. Zou, A. Vishwanath, and T. Senthil, Phys. Rev. X 8, 031089 (2018).
  * (6) F. Wu, A. H. MacDonald, and I. Martin, Phys. Rev. Lett. 121, 257001 (2018).
  * (7) B. Roy and V. Juricic, Phys. Rev. B 99, 121407(R) (2019).
  * (8) B. Lian, Z. Wang, and B. A. Bernevig, Phys. Rev. Lett. 122, 257002 (2019).
  * (9) P. Stepanov, I. Das, X. Lu, A. Fahimniya, K. Watanabe, T. Taniguchi, F. H. L. Koppens, J. Lischner, L. Levitov and D. K. Efetov, Nature 583, 375 (2020).
  * (10) Y. Saito, J. Ge, K. Watanabe, T. Taniguchi, A. F. Young, Nature Physics 16, 926 (2020).
  * (11) Y. Cao, D. Rodan-Legrain, J. M. Park, F. N. Yuan, K. Watanabe, T. Taniguchi, R. M. Fernandes, L. Fu, P. Jarillo-Herrero, Science, 372, 264 (2021).
  * (12) Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras and P. Jarillo-Herrero, Nature 556, 43 (2018).
  * (13) M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Science 363, 1059 (2019).
  * (14) X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir, I. Das, C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang, A. Bachtold, A. H. MacDonald, and D. K. Efetov, Nature 574, 653 (2019).
  * (15) H. Polshyn, M. Yankowitz, S. Chen, Y. Zhang, K. Watanabe, T. Taniguchi, C. R. Dean, and A. F. Young, Nat. Phys. bf 15, 1011 (2019).
  * (16) Y. Cao, D. Chowdhury, D. Rodan-Legrain, O. Rubies-Bigorda, K. Watanabe, T. Taniguchi, T. Senthil, and P. Jarillo-Herrero, Phys. Rev. Lett. 124, 076801 (2020).
  * (17) S. L. Tomarken, Y. Cao, A. Demir, K. Watanabe, T. Taniguchi, P. Jarillo-Herrero, and R. C. Ashoori Phys. Rev. Lett. 123, 046601 (2019).
  * (18) G. Li, A. Luican, J. M. B. Lopes dos Santos, A. H. Castro Neto, A. Reina, J. Kong and E. Y. Andrei, Nature Physics, 6, 109 (2010).
  * (19) A. Kerelsky, L. J. McGilly, D. M. Kennes, L. Xian, M. Yankowitz, S. Chen, K. Watanabe, T. Taniguchi, J. Hone, C. Dean, A. Rubio and A. N. Pasupathy, Nature, 572, 95 (2019).
  * (20) Y. Xie, B. Lian, B. Jäck, X. Liu, C.-L. Chiu, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Nature 572, 101 (2019).
  * (21) Y. Choi, J. Kemmer, Y. Peng, A. Thomson, H. Arora, R. Polski, Y. Zhang, H. Ren, J. Alicea, G. Refael, F. von Oppen, K. Watanabe, T. Taniguchi, and S. Nadj-Perge, Nat. Phys. 15, 1174 (2019).
  * (22) D. Wong, K. P. Nuckolls, M. Oh, B. Lian, Y. Xie, S. Jeon, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Nature 582, 198 (2020).
  * (23) Y. Choi, H. Kim, Y. Peng, A. Thomson, C. Lewandowski, R. Polski, Y. Zhang, H. S. Arora, K. Watanabe, T. Taniguchi, J. Alicea, and S. Nadj-Perge, Nature 589, 536 (2021).
  * (24) N. Tilak, X. Lai, S. Wu, Z. Zhang, M. Xu, R. de Almeida Ribeiro, P. C. Canfield and E. Y. Andrei, Nat Commun bf 12, 4180 (2021).
  * (25) D. Călugăru, N. Regnault, M. Oh, K. P. Nuckolls, D. Wong, R. L. Lee, A. Yazdani, O. Vafek, and B. A. Bernevig, Phys. Rev. Lett. 129, 117602 (2022).
  * (26) M. I. B. Utama, R. J. Koch, K. Lee, et al. Nat. Phys. 17, 184 (2021).
  * (27) S. Lisi, X. Lu, T. Benschop, et al., Nat. Phys. bf 17, 189 (2021).
  * (28) K. Sato, N. Hayashi, T. Ito,et al., Commun Mater 2, 117 (2021).
  * (29) Y. Tserkovnyak, A. Brataas, and G. E. W. Bauer, Phys. Rev. Lett. 88, 117601 (2002).
  * (30) Y. Tserkovnyak, A. Brataas, G. E. W. Bauer, and B. I. Halperin, Rev. Mod. Phys. 77, 1375 (2005).
  * (31) M. Sadamichi and others (eds.), Spin Current, 1st edn, Series on Semiconductor Science and Technology (Oxford, 2012; online edn, Oxford Academic, 17 Dec. 2013).
  * (32) F. Hellman, A. Hoffmann, Y. Tserkovnyak, G. S. D. Beach, E. E. Fullerton, C. Leighton, A. H. MacDonald, D. C. Ralph, D. A. Arena, H. A. Dürr, P. Fischer, J. Grollier, J. P. Heremans, T. Jungwirth, A. V. Kimel, B. Koopmans, I. N. Krivorotov, S. J. May, A. K. Petford-Long, J. M. Rondinelli, N. Samarth, I. K. Schuller, A. N. Slavin, M. D. Stiles, O. Tchernyshyov, A. Thiaville, and B. L. Zink, Rev. Mod. Phys. 89, 025006 (2017).
  * (33) Z. Qiu, J. Li, D. Hou, E. Arenholz, A. T. N’Diaye, A. Tan, K.-i. Uchida, K. Sato, S. Okamoto, Y. Tserkovnyak, Z. Q. Qiu, and E. Saitoh, Nat. Commun. 7, 12670 (2016).
  * (34) F. Yang and P. C. Hammel, J. Phys. D Appl. Phys. 51, 253001 (2018).
  * (35) W. Han, S. Maekawa, and X. Xie, Nat. Mater. 19, 139 (2020).
  * (36) S. Hait et al., ACS Appl. Mater. Interfaces 14, 37182 (2022).
  * (37) 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, Nat. Commun. 5, 4875 (2014).
  * (38) Z. Wang, D.-K. Ki, H. Chen, H. Berger, A. H. MacDonald, and A. F. Morpurgo, Nat. Commun. 6, 8339 (2015).
  * (39) Z. Wang, D.-K. Ki, J. Y. Khoo, D. Mauro, H. Berger, L. S. Levitov, and A. F. Morpurgo, Phys. Rev. X 6, 041020 (2016).
  * (40) B. Yang, M.-F. Tu, J. Kim, Y. Wu, H. Wang, J. Alicea, R. Wu, M. Bockrath, and J. Shi, 2D Mater. 3, 031012 (2016).
  * (41) W. Yan, O. Txoperena, R. Llopis, H. Dery, L. E. Hueso, and F. Casanova, Nat. Commun. 7, 13372 (2016).
  * (42) B. Yang, M. Lohmann, D. Barroso, I. Liao, Z. Lin, Y. Liu, L. Bartels, K. Watanabe, T. Taniguchi, and J. Shi, Phys. Rev. B 96, 041409(R) (2017).
  * (43) T. S. Ghiasi, J. Ingla-Aynés, A. A. Kaverzin, and B. J. van Wees, Nano Lett. 17, 7528 (2017).
  * (44) A. Dankert and S. P. Dash, Nat. Commun. 8, 16093 (2017).
  * (45) T. Völkl, T. Rockinger, M. Drienovsky, K. Watanabe, T. Taniguchi, D. Weiss, and J. Eroms, Phys. Rev. B 96, 125405 (2017).
  * (46) S. Zihlmann, A. W. Cummings, J. H. Garcia, M. Kedves, K. Watanabe, T. Taniguchi, C. Schönenberger, and P. Makk, Phys. Rev. B 97, 075434 (2018).
  * (47) T. Wakamura, F. Reale, P. Palczynski, S. Guéron, C. Mattevi, and H. Bouchiat, Phys. Rev. Lett. 120, 106802 (2018).
  * (48) J. C. Leutenantsmeyer, J. Ingla-Aynés, J. Fabian, and B. J. van Wees, Phys. Rev. Lett. 121, 127702 (2018).
  * (49) S. Omar and B. J. van Wees, Phys. Rev. B 97, 045414 (2018).
  * (50) L. A. Benítez, J. F. Sierra, W. S. Torres, A. Arrighi, F. Bonell, M. V. Costache, and S. O. Valenzuela, Nat. Phys. 14, 303 (2018).
  * (51) C. K. Safeer, J. Ingla-Aynés, F. Herling, J. H. Garcia, M. Vila, N. Ontoso, M. R. Calvo, S. Roche, L. E. Hueso, and F. Casanova, Nano Lett. 19, 1074 (2019).
  * (52) 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, Nature 571, 85 (2019).
  * (53) D. Wang, S. Che, G. Cao, R. Lyu, K. Watanabe, T. Taniguchi, C. Ning Lau, and M. Bockrath, Nano Lett. 19, 7028 (2019).
  * (54) T. Wakamura, F. Reale, P. Palczynski, M. Q. Zhao, A. T. C. Johnson, S. Guéron, C. Mattevi, A. Ouerghi, and H. Bouchiat, Phys. Rev. B 99, 245402 (2019).
  * (55) A. David, P. Rakyta, A. Kormányos, and G. Burkard Phys. Rev. B bf 100, 085412(2019).
  * (56) T. Wang, N. Bultinck, and M. P. Zaletel Phys. Rev. B 102, 235146 (2020).
  * (57) For a rieview see, T. Wakamura, S. Guéron and H. Bouchiat, Comptes Rendus. Physique, 22, 145 (2021).
  * (58) H. S. Arora, R. Polski, Y. Zhang, A. Thomson, Y. Choi, H. Kim, Z. Lin, I. Z. Wilson, X. Xu, J. -H. Chu, K. Watanabe, T. Taniguchi, J. Alicea and S. Nadj-Perge, Nature 583 379 (2020).
  * (59) J.-X. Lin, Y.-H. Zhang, E. Morissette, Z. Wang, S. Liu, D. Rhodes, K. Watanabe, T. Taniguchi, J. Hone, J. I. A. Li, Science 375, 437 (2022).
  * (60) S. Bhowmik, B. Ghawri, Y. Park, D. Lee, S. Datta, R. Soni, K. Watanabe, T. Taniguchi, A. Ghosh, J. Jung, U. Chandni, arXiv:2211.01251.
  * (61) Y. Ohnuma, H. Adachi, E. Saitoh, and S. Maekawa Phys. Rev. B 89, 174417 (2014).
  * (62) M. Matsuo, Y. Ohnuma, T. Kato, and S. Maekawa, Phys. Rev. Lett. 120, 037201 (2018).
  * (63) T. Kato, Y. Ohnuma, M. Matsuo, J. Rech, T. Jonckheere, and T. Martin, Phys. Rev. B 99, 144411 (2019).
  * (64) T. Kato, Y. Ohnuma, and M. Matsuo, Phys. Rev. B 102, 094437 (2020).
  * (65) Y. Ominato and M. Matsuo, J. Phys. Soc. Jpn. 89, 053704 (2020).
  * (66) Y. Ominato, J. Fujimoto, and M. Matsuo, Phys. Rev. Lett. 124, 166803 (2020).
  * (67) M. Yama, M. Tatsuno, T. Kato, and M. Matsuo Phys. Rev. B 104, 054410 (2021).
  * (68) M. Gmitra, and J. Fabian, Phys. Rev. B 92, 155403 (2015).
  * (69) M. Gmitra, D. Kochan, P. Hogl, and J. Fabian, Phys. Rev. B 93, 155104 (2016).
  * (70) S. Zihlmann, et al., Phys. Rev. B 97, 075434 (2018).
  * (71) J. O. Island, et al., Nature 571, 85 (2019).
  * (72) D. Wang, et al., Nano Lett. 19, 7028 (2019).
  * (73) Z. Wang et al., Phys. Rev. X 6, 041020 (2016).
  * (74) A. W. Cummings, J. H. Garcia, J. Fabian, J. and S. Roche, Phys. Rev. Lett. 119, 206601 (2017).
  * (75) R. Bistritzer and A. H. MacDonald, Proc. Natl. Acad. Sci.U.S.A., 108, 12233 (2011).
  * (76) For details, see the Supplemental Material.
  * (77) G. Cantele, D.Alfè, F. Conte, V. Cataudella, D. Ninno, and P. Lucignano, Phys. Rev. Research 2, 043127 (2020).
  * (78) T. Funato, T. Kato, and M. Matsuo, Phys. Rev. B 106, 144418 (2022).
  * (79) M. Koshino, N. F. Q. Yuan, T. Koretsune, M. Ochi, K. Kuroki, and L. Fu, Phys. Rev. X 8, 031087 (2018).
  * (80) D. A. Ruiz-Tijerina and V. I. Fal’ko, Phys. Rev. B 99, 125424 (2019).
  * (81) S. Zhang, A. Song, L. Chen, C. Jiang, C. Chen, L. Gao, Y. Hou, L. Liu, T. Ma, H. Wang, X.-Q. Feng and Q. Li, Sc. Adv. eabc5555 6 (2020).
  * (82) A. C. Gadelha, D. A. A. Ohlberg, C. Rabelo, et al., Nature 590, 405 (2021).
  * (83) Z. Bi, N. F. Q. Yuan and L. Fu, Phys. Rev. B 100, 035448 (2019).

Supplemental material for Twisted bilayer graphene reveals its flat bands
under spin pumping

##  S.1 I. Derivation of the low-energy Hamiltonian of TBG without SOC

We start by a brief overview of the perturbative approach proposed by
Bistritzer and MacDonald Mc11  to derive the continuum model of TBG. We
consider a TBG where the two layers $l=1,2$ are rotated oppositely
$\theta_{2}=-\theta_{1}=\frac{\theta}{2}$. The Hamiltonian of a graphene layer
$l$ rotated at an angle $\theta_{l}$ is

$\displaystyle
h_{l,rot}(\mathbf{k})=e^{i\frac{\theta_{l}}{2}\sigma_{z}}h_{l}(\mathbf{k})e^{-i\frac{\theta_{l}}{2}\sigma_{z}}$
(S.1)

where $h_{l}(\mathbf{k})$ is the Hamiltonian of the unrotated layer $(l)$
given in the continuum limit, by

$\displaystyle h_{l}(\mathbf{k})=-\hbar
v_{F}\mathbf{k}\cdot\mathbf{\sigma}^{\ast},$ (S.2)

where the momentum $\mathbf{k}$ is written relatively to $\mathbf{K}_{1,\xi}$,
$v_{F}$ is the Fermi velocity, $\xi$ is the valley index,
$\mathbf{\sigma}^{\ast}=\left(\xi\sigma_{x},\sigma_{y}\right)$ and
$\sigma_{i}$ ($i=x,y,z$) are the sublattice-Pauli matrices.

The leading contribution of the interlayer tunneling can be limited to three
nearest hopping processes in the momentum space connecting states
$|\mathbf{k}\rangle_{1}$, around the Dirac point $\mathbf{K}_{1,\xi}$ of layer
$(1)$, to the states $|\mathbf{k+q_{j\xi}}\rangle_{2}$ around
$\mathbf{K}_{2,\xi}$, the Dirac point of layer $(2)$. The $\mathbf{q_{j\xi}}$
vectors are given by Mc11

$\displaystyle\mathbf{q}_{1\xi}=\xi
k_{\theta}\left(0,1\right),\;\mathbf{q}_{2\xi}=\mathbf{q}_{1\xi}+\xi\mathbf{G}^{M}_{1}=\xi
k_{\theta}\left(-\frac{\sqrt{3}}{2},-\frac{1}{2}\right),$
$\displaystyle\mathbf{q_{3}}=\mathbf{q}_{1\xi}+\xi\left(\mathbf{G}^{M}_{1}+\mathbf{G}^{M}_{2}\right)=\xi
k_{\theta}\left(\frac{\sqrt{3}}{2},-\frac{1}{2}\right),$ (S.3)

where $k_{\theta}=2k_{D}\sin\frac{\theta}{2}\sim\theta k_{D}$ and
$k_{D}=|\mathbf{K}_{1,\xi}|=|\mathbf{K}_{2,\xi}|=\frac{4\pi}{3a}$, $a$ being
the graphene lattice parameter. The
$\left(\mathbf{G}^{M}_{1},\mathbf{G}^{M}_{2}\right)$ is the moiré BZ basis
given by $\mathbf{G}^{M}_{i}=\mathcal{R}_{t}^{T}\mathbf{G}_{i}$,
$\mathbf{G}_{i}$ are the lattice basis vectors of the monolayer reciprocal
lattice $\mathbf{G}_{1}=\frac{2\pi}{a}\left(1,-1/\sqrt{3}\right)$ and
$\mathbf{G}_{2}=\frac{2\pi}{a}\left(0,2/\sqrt{3}\right)$. $\mathcal{R}_{t}$ is
the rotation tensor written, in the sublattice basis, at a small twist angle
as

$\displaystyle R(\theta)=\begin{pmatrix}0&-\theta\\\ \theta&0\end{pmatrix}.$
(S.4)

In the basis
$\left\\{|\mathbf{k}\rangle_{1},|\mathbf{k+q_{j,\xi}}\rangle_{2}\right\\}$,
the Hamiltonian, at the valley $\xi$, reads as Mc11

$\displaystyle
H(\mathbf{k})=\begin{pmatrix}h_{1}(\mathbf{k})&T_{1}&T_{2}&T_{3}\\\
T^{\dagger}_{1}&h_{2,1}(\mathbf{k})&0&0\\\
T^{\dagger}_{2}&0&h_{2,2}(\mathbf{k})&0\\\
T^{\dagger}_{3}&0&0&h_{2,3}(\mathbf{k})\\\ \end{pmatrix},$ (S.5)

For the relaxed TBG the $T_{j}$ matrices are given by Koshino18

$\displaystyle T_{1}=\begin{pmatrix}w&w^{\prime}\\\
w^{\prime}&w^{\prime\prime}\
\end{pmatrix},T_{2}=e^{i\xi\mathbf{G}^{M}_{1}\cdot\mathbf{r}}\begin{pmatrix}w&w^{\prime}e^{-i\xi\Phi}\\\
w^{\prime}e^{i\xi\Phi}&w^{\prime\prime}\\\
\end{pmatrix},T_{3}=e^{i\xi\left(\mathbf{G}^{M}_{1}+\mathbf{G}^{M}_{2}\right)\cdot\mathbf{r}}\begin{pmatrix}w&w^{\prime}e^{i\xi\Phi}\\\
w^{\prime}e^{-i\xi\Phi}&w^{\prime\prime}\ \end{pmatrix},$ (S.6)

Here $h_{1}(\mathbf{k})$ is given by Eq. S.2 and $h_{2,j}(\mathbf{k})\equiv
h_{2}(\mathbf{k}+\mathbf{q}_{j\xi})=h_{1}(\mathbf{k}+\mathbf{q}_{j\xi})$,
where the momentum $\mathbf{k}$ is written relatively to $\mathbf{K}_{1,\xi}$.
$\Phi=\frac{2\pi}{3}$ and $\mathbf{r}$ is the shortest inplane shifts between
carbon atoms of the two layers Falko19 . We consider the AA stacking domains
hosting the flat bands of TGB Zhang ; Gadelha . In this configuration, each A
atom of the top layer is directly located above an A atom of the bottom layer
and $\mathbf{r}=\mathbf{0}$ Falko19 . The parameters $w$, $w^{\prime}$ and
$w^{\prime\prime}$ are the tunneling amplitudes which have the same value
$w=w^{\prime}=w^{\prime\prime}\sim 118\;\mathrm{meV}$ in the rigid TBG Bi . In
the TBG relaxed lattice, these amplitudes are no more equal $w\sim
w^{\prime\prime}\sim 90$ meV and $w^{\prime}=117$ meV. In the present work, we
do not consider the lattice relaxation effect, since the SOC parameters
$\lambda_{I},\;\lambda_{R}\sim 2\,\mathrm{meV}$ are small compared to the
difference between the interlayer amplitudes $\Delta w=w-w^{\prime}\sim
20\,\mathrm{meV}$. In the unrelaxed AA domains, the $T_{j}$ matrices can be
written as

$\displaystyle
T_{1}=w\left(\mathbb{I}_{\sigma}+\sigma_{x}\right),T_{2}=w\left(\mathbb{I}_{\sigma}-\frac{1}{2}\sigma_{x}+\xi\frac{\sqrt{3}}{2}\sigma_{y}\right),T_{3}=w\left(\mathbb{I}_{\sigma}-\frac{1}{2}\sigma_{x}-\xi\frac{\sqrt{3}}{2}\sigma_{y}\right).$
(S.7)

Here, $\mathbb{I}_{\sigma}$ is the identity matrix acting on the sublattice
indices.

Considering, in Eq. S.5, the $\mathbf{k}$ dependent term as a perturbation,
the effective Hamiltonian can be written, to the leading order in
$\mathbf{k}$, as

$\displaystyle
H^{(1)}\left(\mathbf{k}\right)=\frac{\langle\Psi|H(\mathbf{k})|\Psi\rangle}{\langle\Psi|\Psi\rangle},$
(S.8)

where
$\Psi=\left(\psi_{0}(\mathbf{k}),\psi_{1}(\mathbf{k}),\psi_{2}(\mathbf{k}),\psi_{3}(\mathbf{k})\right)$
is the zero energy eigenstate of $H\left(\mathbf{k}=\mathbf{0}\right)$. $\Psi$
is constructed on the two-component sublattice spinor $\psi_{0}(\mathbf{k})$
($\psi_{j}(\mathbf{k})$) of layer $1$ (layer $2$) taken at the momentum
$\mathbf{k}$ ($\mathbf{k}+\mathbf{q}_{j\xi}$) around the Dirac point
$\mathbf{K}_{1,\xi}$ ($\mathbf{K}_{2,\xi}$) at the valley $\xi$. $\psi_{0}$ is
the zero energy eigenstate of $h_{1}$. Then, the $\Psi$ components satisfy

$\displaystyle
h_{1}\psi_{0}+\sum_{j}T_{j}\psi_{j}=0,\,\mathrm{and}\;T^{\dagger}_{j}\psi_{0}+h_{j}\psi_{j}=0,\;\mathrm{with}\;h_{1}\psi_{0}=0,$
(S.9)

where $h_{j}\equiv h_{2}(\mathbf{q}_{j,\xi})$. Then

$\displaystyle\psi_{j}=-h^{-1}_{j}T^{\dagger}_{j}\psi_{0},\mathrm{and}\;\sum_{j}T_{j}h^{-1}_{j}T^{\dagger}_{j}=0.$
(S.10)

To the leading order in $\mathbf{k}$, $H^{(1)}\left(\mathbf{k}\right)$, takes
the following form

$\displaystyle
H^{(1)}\left(\mathbf{k}\right)=\frac{\langle\Psi|H(\mathbf{k})|\Psi\rangle}{\langle\Psi|\Psi\rangle}=\frac{1}{\langle\Psi|\Psi\rangle}\left[\psi^{\dagger}_{0}h_{1}\left(\mathbf{k}\right)\psi_{0}+\psi^{\dagger}_{0}\sum_{j}T_{j}h^{-1}_{j}h_{j}\left(\mathbf{k}\right)h^{-1}_{j}T^{\dagger}_{j}\psi_{0}\right],$
(S.11) $\displaystyle H^{(1)}\left(\mathbf{k}\right)=-\hbar
v^{\ast}\psi^{\dagger}_{0}\,\mathbf{k}\cdot\mathbf{\sigma^{\ast}}\,\psi_{0},$
(S.12)

Eq. S.12 is obtained by neglecting $\theta$ in $h_{j}\left(\mathbf{k}\right)$,
which turns out to take $\mathbf{q}_{j,\xi}=\mathbf{0}$ in
$h_{j}\left(\mathbf{k}\right)$ Mc11 .
$v^{\ast}$ is the effective velocity of the energy band of TBG around the zero
energy which vanishes at the MA $\theta_{m}$, and is given by Mc11

$\displaystyle v^{\ast}=v_{F}\frac{1-3\alpha^{2}}{1+6\alpha^{2}}$ (S.13)

where $\alpha=\frac{w}{\hbar v_{F}q_{0}}$,
$q_{0}=|\mathbf{q}_{j\xi}|\sim\frac{4\pi}{3a}\theta$.
In our numerical calculations, we take $w=118\;\mathrm{meV}$ and $\hbar
v_{F}/a\sim 2.68\;\mathrm{eV}$ which corresponds to $\theta_{m}=1.05^{\circ}$
for the first MA Mc11 .

##  S.2 II. Derivation of the low-energy Hamiltonian of TBG with SOC

We now consider the heterostructure consisting of TBG adjacent to a monolayer
of $\mathrm{WSe_{2}}$ as shown in Fig.1 of the main text, where we denote the
graphene layer in contact with the TMD by layer (2). This layer is subject to
a SOC induced by proximity effect by the TMD, and the corresponding
Hamiltonian can be written as

$\displaystyle
h_{2}(\mathbf{k})=h_{1}(\mathbf{k})+h_{\text{SOC}}+\frac{m}{2}\sigma_{z}$
(S.14)

where $h_{\text{SOC}}$ is given by Eq. 3 of the main text.

To derive the continuum model of TBG/$\mathrm{WSe}_{2}$, we follow the
perturbative approach of Ref. Mc11  presented in the previous section. Now,
the basis
$\Psi=\left(\psi_{0}(\mathbf{k}),\psi_{1}(\mathbf{k}),\psi_{2}(\mathbf{k}),\psi_{3}(\mathbf{k})\right)$
is constructed on the four-component spin-sublattice spinor
$\psi_{0}(\mathbf{k})$ and $\psi_{j}(\mathbf{k})$, ($j=1,2,3$) corresponding,
respectively, to layer $(1)$ and layer $(2)$. $\psi_{0}(\mathbf{k})$ is
written as
$\psi_{0}(\mathbf{k})^{T}=\left(\psi_{0,A\uparrow},\psi_{0,A\downarrow},\psi_{0,B\uparrow},\psi_{0,B\downarrow}\right)$.
In this basis, the Hamiltonian of TBG/$\mathrm{WSe}_{2}$ takes the form

$\displaystyle
H_{\xi,\text{SOC}}(\mathbf{k})=\begin{pmatrix}h_{1}(\mathbf{k})&T_{1}&T_{2}&T_{3}\\\
T^{\dagger}_{1}&h_{2,1}(\mathbf{k})&0&0\\\
T^{\dagger}_{2}&0&h_{2,2}(\mathbf{k})&0\\\
T^{\dagger}_{3}&0&0&h_{2,3}(\mathbf{k})\\\ \end{pmatrix}.$ (S.15)

The momentum $\mathbf{k}$ is measured relatively to the Dirac point
$\mathbf{K}_{1\xi}$ of layer (1),
$h_{2,j}\left(\mathbf{k}\right)=h_{2}\left(\mathbf{k}+\mathbf{q}_{j\xi}\right)$,
($j=1,2,3$) and $h_{2}\left(\mathbf{k}\right)$ includes now the SOC terms (Eq.
S.14).
We consider AA stacked moiré domains with spin independent interlayer hopping.
The $T_{j}$ matrices are written as the tensor product of those given by Eq.
S.7, with the $2\times 2$ identity spin-matrix $\mathbb{I}_{s}$.

Regarding the small values of the SOC, we assume that
$H_{\xi,\text{SOC}}(\mathbf{k})$ has a zero eigenenergy and the corresponding
eigenstate $\Psi$ satisfies the condition given by Eq. S.10.

Following the same procedure as in the previous section, we derive from Eq.
S.11 the effective low energy Hamiltonian
$H_{\xi,\text{SOC}}^{(1)}(\mathbf{k})$ of TBG/$\mathrm{WSe_{2}}$ by
substituting $h_{j}\left(\mathbf{k}\right)$ by the Hamiltonian of layer (2),
rotated at $\theta/2$ and including SOC as

$\displaystyle
h_{\text{2,rot}}(\mathbf{k})=-\hbar\mathbf{k}\cdot\mathbf{\sigma}^{\ast}\mathbb{I}_{s}+\frac{\lambda_{I}}{2}\xi
s_{z}+\frac{\lambda_{R}}{2}\left(\xi\sigma_{x}s_{y}-\sigma_{y}s_{x}\right)+\frac{\lambda_{KM}}{2}\xi\sigma_{z}s_{z}-\frac{\lambda_{R}}{2}\theta\left(\xi\sigma_{y}s_{y}+\sigma_{x}s_{x}\right),$
(S.16)

Hereafter, we neglect the Kane and Mele term whose contribution, to the
leading order in $\mathbf{k}$, is found to vanish. We also disregard the last
term in Eq. S.16, which results into a higher order correction in $\theta$.

To the first order in the SOC coupling, we obtain the continuum model of
TBG/$\mathrm{WSe_{2}}$ described by the Hamiltonian
$H_{\xi,\text{SOC}}^{(1)}(\mathbf{k})$ given by Eq. 5 in the main text. This
Hamiltonian contain a SOC term ($h_{\text{eff}}^{\text{SOC}}$) with
renormalized Ising and Rashba interactions

$\displaystyle\tilde{\lambda}_{I}=\frac{6\alpha^{2}}{\langle\Psi|\Psi\rangle}\lambda_{I},\;\tilde{\lambda}_{R}=\frac{3\alpha^{2}}{\langle\Psi|\Psi\rangle}\lambda_{R},\;\mathrm{and}\;\langle\Psi|\Psi\rangle\sim
1+6\alpha^{2}$ (S.17)

$\tilde{\lambda}_{I}$ and $\tilde{\lambda}_{R}$ are enhanced by decreasing the
twist angle from the MA.

Figure S.1: Energy bands of TBG/$\mathrm{WSe}_{2}$ around zero energy as
function of the dimensionless momentum amplitude $k/q_{0}$ at different twist
angles. The bands are represented up to the cutoff $k_{c}=q_{0}/2$.
Calculations are done for $\lambda_{I}=3\;\mathrm{meV}$ and
$\lambda_{R}=4\;\mathrm{meV}$ Alex .

To the leading order in $\mathbf{k}$, the four eigenergies of the Hamiltonian
$H_{\xi,\text{SOC}}^{(1)}$ (Eq. 5 of the main text), denoted
$E(\mathbf{k})_{\sigma,\pm}$, are given by

$\displaystyle
E(\mathbf{k}))_{\sigma,\pm}=\frac{\sigma}{\langle\Psi|\Psi\rangle}\sqrt{f_{1}(\mathbf{k}))\pm
6\alpha^{2}\sqrt{f_{2}(\mathbf{k}))}},$ (S.18) $\displaystyle
f_{1}(\mathbf{k}))=(\hbar
v)^{2}\left(1-3\alpha^{2}\right)^{2}||{\mathbf{k}}||^{2}+\frac{9}{2}\alpha^{4}\left(2\lambda_{I}^{2}+\lambda_{R}^{2}\right),$
$\displaystyle f_{2}(\mathbf{k}))=(\hbar
v)^{2}\left(1-3\alpha^{2}\right)^{2}||{\mathbf{k}}||^{2}\left(\lambda_{I}^{2}+\frac{1}{4}\lambda_{R}^{2}\right)+\frac{9}{16}\alpha^{4}\lambda_{R}^{4},$
(S.19)

where $\sigma=\pm$ is the band index. $E(\mathbf{k})_{\sigma,\pm}$ are
depicted in Fig. S.1 at different twist angles.

At the Dirac point, the eigenenergies reduce to

$\displaystyle E_{\pm,+}=\frac{\pm
3\alpha^{2}}{\langle\Psi|\Psi\rangle}\sqrt{\lambda_{I}^{2}+\lambda_{R}^{2}}\;,E_{\pm,-}=\frac{\pm
3\alpha^{2}}{\langle\Psi|\Psi\rangle}\lambda_{I}.$ (S.20)

##  S.3 Electronic Green function

The Matsubara Green function associated to the effective Hamiltonian
$H_{\xi,\text{SOC}}^{(1)}$ (Eq. 5) is

$\displaystyle\hat{g}(\mathbf{k},i\omega_{n})=\left[i\omega_{n}\mathbb{I}_{S}\mathbb{I}_{\sigma}-H_{\xi,\text{SOC}}(\mathbf{k})\right]^{-1},$
(S.21)

where $\mathbb{I}_{S}$ and $\mathbb{I}_{\sigma}$ are the $2\times 2$ spin and
band identity matrices, respectively. $\hat{g}(\mathbf{k},i\omega_{n})$ can be
expressed as

$\displaystyle\hat{g}(\mathbf{k},i\omega_{n})=\hat{g}_{0}(\mathbf{k},i\omega_{n})\mathbb{I}_{s}+\mathbf{\hat{g}}\cdot\mathbf{s}.$
(S.22)

$\mathbf{s}=(s_{x},s_{y},s_{z})$ are spin-Pauli matrices,
$\hat{g}_{0}(\mathbf{k},i\omega_{n})$ and the components $\hat{g}_{i}$
($i=x,y,z$) of $\mathbf{\hat{g}}$ are written as

$\displaystyle\hat{g}_{0}(\mathbf{k},i\omega_{n})=A_{0}(\mathbf{k},i\omega_{n})+C_{0}(\mathbf{k},i\omega_{n}),\,\hat{g}_{x}(\mathbf{k},i\omega_{n})=B_{x}(\mathbf{k},i\omega_{n})+D_{x}(\mathbf{k},i\omega_{n}),\,$
$\displaystyle
g_{y}(\mathbf{k},i\omega_{n})=B_{y}(\mathbf{k},i\omega_{n})+D_{y}(\mathbf{k},i\omega_{n}),\,\hat{g}_{z}(\mathbf{k},i\omega_{n})=A_{z}(\mathbf{k},i\omega_{n})+C_{z}(\mathbf{k},i\omega_{n}).$
(S.23)

The $A,B,C$ and $D$ operators functions are written in terms of the band-Pauli
matrices $\sigma_{x,y,z}$ and the corresponding identity matrix
$\mathbb{I}_{\sigma}$

$\displaystyle
A_{i}(\mathbf{k},i\omega_{n})=A_{i1}(\mathbf{k},i\omega_{n})\mathbb{I}_{\sigma}+A_{iz}(\mathbf{k},i\omega_{n})\sigma_{z}$
$\displaystyle
C_{i}(\mathbf{k},i\omega_{n})=C_{ix}(\mathbf{k},i\omega_{n})\sigma_{x}+C_{iy}(\mathbf{k},i\omega_{n})\sigma_{y}$
$\displaystyle
B_{j}(\mathbf{k},i\omega_{n})=B_{j1}(\mathbf{k},i\omega_{n})\mathbb{I}_{\sigma}+B_{jz}(\mathbf{k},i\omega_{n})\sigma_{z}$
$\displaystyle
D_{j}(\mathbf{k},i\omega_{n})=D_{jx}(\mathbf{k},i\omega_{n})\sigma_{x}+D_{jy}(\mathbf{k},i\omega_{n})\sigma_{y}$
(S.24)

here $j=x,y$ and $i=0,z$.

In the limit of small SOC couplings $\lambda_{R},\lambda_{I}\ll\hbar vq_{0}$,
A, B, C, and D become

$\displaystyle A_{01}(\mathbf{k},i\omega_{m})$
$\displaystyle=\frac{Y_{00}-E_{1}^{2}}{2(E_{3}^{2}-E_{1}^{2})}\left[\frac{1}{i\hbar\omega_{m}-E_{1}}+\frac{1}{i\hbar\omega_{m}-E_{2}}\right]+\frac{E_{3}^{2}-Y_{00}}{2(E_{3}^{2}-E_{1}^{2})}\left[\frac{1}{i\hbar\omega_{m}-E_{3}}+\frac{1}{i\hbar\omega_{m}-E_{4}}\right]$
$\displaystyle Y_{00}=\frac{(\hbar
v)^{2}(1-3\alpha^{2})^{2}}{\langle\Psi|\Psi\rangle^{2}}||\mathbf{k}||^{2}+\frac{1}{2}\left(\frac{3\alpha^{2}\lambda_{R}}{\langle\Psi|\Psi\rangle}\right)^{2}+\left(\frac{3\alpha^{2}\lambda_{I}}{\langle\Psi|\Psi\rangle}\right)^{2}$
$\displaystyle A_{0z}(\mathbf{k},i\omega_{m})$
$\displaystyle=\frac{Y_{0}}{2\left(E_{1}^{2}-E_{3}^{2}\right)}\left\\{\frac{1}{E_{1}}\left(\frac{1}{i\hbar\omega_{m}-E_{1}}-\frac{1}{i\hbar\omega_{m}-E_{2}}\right)-\frac{1}{E_{3}}\left(\frac{1}{i\hbar\omega_{m}-E_{3}}-\frac{1}{i\hbar\omega_{m}-E_{4}}\right)\right\\}$
$\displaystyle
Y_{0}=-\frac{3}{2}\left(\frac{3\alpha^{2}\lambda_{R}}{\langle\Psi|\Psi\rangle}\right)^{2}\frac{\alpha^{2}\lambda_{I}}{\langle\Psi|\Psi\rangle}$
$\displaystyle C_{0x}(\mathbf{k},i\omega_{m})$
$\displaystyle=\frac{1}{2\left(E_{1}^{2}-E_{3}^{2}\right)}\left\\{Y_{2}\left[\frac{1}{E_{1}}\left(\frac{1}{i\hbar\omega_{m}-E_{1}}-\frac{1}{i\hbar\omega_{m}-E_{2}}\right)-\frac{1}{E_{3}}\left(\frac{1}{i\hbar\omega_{m}-E_{3}}-\frac{1}{i\hbar\omega_{m}-E_{4}}\right)\right]\right.$
$\displaystyle\left.-Y_{1}\left[E_{1}\left(\frac{1}{i\hbar\omega_{m}-E_{1}}-\frac{1}{i\hbar\omega_{m}-E_{2}}\right)-E_{3}\left(\frac{1}{i\hbar\omega_{m}-E_{3}}-\frac{1}{i\hbar\omega_{m}-E_{4}}\right)\right]\right\\}\xi
k_{x}$ $\displaystyle Y_{1}=\frac{\hbar
v(1-3\alpha^{2})}{\langle\Psi|\Psi\rangle},\;Y_{2}=-Y_{1}\left(\frac{3\alpha^{2}\lambda_{I}}{\langle\Psi|\Psi\rangle}\right)^{2}$
$\displaystyle C_{0y}(\mathbf{k},i\omega_{m})$
$\displaystyle=\frac{1}{2\left(E_{1}^{2}-E_{3}^{2}\right)}\left\\{Y_{2}\left[\frac{1}{E_{1}}\left(\frac{1}{i\hbar\omega_{m}-E_{1}}-\frac{1}{i\hbar\omega_{m}-E_{2}}\right)-\frac{1}{E_{3}}\left(\frac{1}{i\hbar\omega_{m}-E_{3}}-\frac{1}{i\hbar\omega_{m}-E_{4}}\right)\right]\right.$
(S.28)
$\displaystyle\left.-Y_{1}\left[\left(E_{1}\frac{1}{i\hbar\omega_{m}-E_{1}}-\frac{1}{i\hbar\omega_{m}-E_{2}}\right)-E_{3}\left(\frac{1}{i\hbar\omega_{m}-E_{3}}-\frac{1}{i\hbar\omega_{m}-E_{4}}\right)\right]\right\\}k_{y}$
$\displaystyle A_{z1}(\mathbf{k},i\omega_{m})$
$\displaystyle=\xi\frac{1}{2\left(E_{1}^{2}-E_{3}^{2}\right)}\left\\{Y_{4}\left[\frac{1}{E_{1}}\left(\frac{1}{i\hbar\omega_{m}-E_{1}}-\frac{1}{i\hbar\omega_{m}-E_{2}}\right)-\frac{1}{E_{3}}\left(\frac{1}{i\hbar\omega_{m}-E_{3}}-\frac{1}{i\hbar\omega_{m}-E_{4}}\right)\right]\right.$
$\displaystyle\left.-Y_{3}\left[E_{1}\left(\frac{1}{i\hbar\omega_{m}-E_{1}}-\frac{1}{i\hbar\omega_{m}-E_{2}}\right)-E_{3}\left(\frac{1}{i\hbar\omega_{m}-E_{3}}-\frac{1}{i\hbar\omega_{m}-E_{4}}\right)\right]\right\\}$
$\displaystyle
Y_{3}=-\frac{3\alpha^{2}\lambda_{I}}{\langle\Psi|\Psi\rangle},\;Y_{4}=-Y_{3}\left[\frac{(\hbar
v)^{2}}{\langle\Psi|\Psi\rangle^{2}}(1-3\alpha^{2})^{2}||\mathbf{k}||^{2}-\frac{1}{2}\left(\frac{3\alpha^{2}\lambda_{R}}{\langle\Psi|\Psi\rangle}\right)^{2}-\left(\frac{3\alpha^{2}\lambda_{I}}{\langle\Psi|\Psi\rangle}\right)^{2}\right]$
$\displaystyle A_{zz}(\mathbf{k},i\omega_{m})$
$\displaystyle=\xi\frac{Y_{5}}{2\left(E_{3}^{2}-E_{1}^{2}\right)}\left\\{\frac{1}{i\hbar\omega_{m}-E_{1}}+\frac{1}{i\hbar\omega_{m}-E_{2}}-\frac{1}{i\hbar\omega_{m}-E_{3}}-\frac{1}{i\hbar\omega_{m}-E_{4}}\right\\}$
$\displaystyle
Y_{5}=\frac{1}{2}\left(\frac{\lambda_{R}}{\langle\Psi|\Psi\rangle}\right)^{2}$
$\displaystyle C_{zx}(\mathbf{k},i\omega_{m})$
$\displaystyle=\frac{Y_{6}}{2\left(E_{3}^{2}-E_{1}^{2}\right)}\left\\{\frac{1}{i\hbar\omega_{m}-E_{1}}+\frac{1}{i\hbar\omega_{m}-E_{2}}-\frac{1}{i\hbar\omega_{m}-E_{3}}-\frac{1}{i\hbar\omega_{m}-E_{4}}\right\\}k_{x}$
$\displaystyle Y_{6}=\hbar
v(1-3\alpha^{2})\frac{6\alpha^{2}\lambda_{I}}{\langle\Psi|\Psi\rangle^{2}}$
$\displaystyle C_{zy}(\mathbf{k},i\omega_{m})$
$\displaystyle=\frac{Y_{6}}{2\left(E_{3}^{2}-E_{1}^{2}\right)}\left\\{\frac{1}{i\hbar\omega_{m}-E_{1}}+\frac{1}{i\hbar\omega_{m}-E_{2}}-\frac{1}{i\hbar\omega_{m}-E_{3}}-\frac{1}{i\hbar\omega_{m}-E_{4}}\right\\}\xi
k_{y}$ $\displaystyle B_{x1}(\mathbf{k},i\omega_{m})$
$\displaystyle=\frac{Y_{7}}{2\left(E_{3}^{2}-E_{1}^{2}\right)}\left\\{\frac{1}{i\hbar\omega_{m}-E_{1}}+\frac{1}{i\hbar\omega_{m}-E_{2}}-\frac{1}{i\hbar\omega_{m}-E_{3}}-\frac{1}{i\hbar\omega_{m}-E_{4}}\right\\}k_{y}$
$\displaystyle Y_{7}=\hbar
v(1-3\alpha^{2})\frac{3\alpha^{2}\lambda_{R}}{\langle\Psi|\Psi\rangle^{2}}$
$\displaystyle B_{xz}(\mathbf{k},i\omega_{m})$
$\displaystyle=Y_{8}\frac{1}{2\left(E_{1}^{2}-E_{3}^{2}\right)}\,k_{y}\times$
$\displaystyle\left\\{\frac{1}{E_{1}}\left[\frac{1}{i\hbar\omega_{m}-E_{1}}-\frac{1}{i\hbar\omega_{m}-E_{2}}\right]+\frac{1}{E_{3}}\left[\frac{1}{i\hbar\omega_{m}-E_{3}}-\frac{1}{i\hbar\omega_{m}-E_{4}}\right]\right\\}$
$\displaystyle Y_{8}=-\hbar
v(1-3\alpha^{2})\frac{3\alpha^{2}\lambda_{R}}{\langle\Psi|\Psi\rangle^{2}}\frac{3\alpha^{2}\lambda_{I}}{\langle\Psi|\Psi\rangle}$
$\displaystyle D_{xx}(\mathbf{k},i\omega_{m})$
$\displaystyle=Y_{9}\frac{1}{2\left(E_{1}^{2}-E_{3}^{2}\right)}\,\xi
k_{x}k_{y}\times$
$\displaystyle\left\\{\frac{1}{E_{1}}\left[\frac{1}{i\hbar\omega_{m}-E_{1}}-\frac{1}{i\hbar\omega_{m}-E_{2}}\right]-\frac{1}{E_{3}}\left[\frac{1}{i\hbar\omega_{m}-E_{3}}-\frac{1}{i\hbar\omega_{m}-E_{4}}\right]\right\\}$
$\displaystyle Y_{9}=(\hbar
v)^{2}(1-3\alpha^{2})^{2}\frac{3\alpha^{2}\lambda_{R}}{\langle\Psi|\Psi\rangle}$
$\displaystyle D_{xy}(\mathbf{k},i\omega_{m})=$
$\displaystyle\frac{1}{2\left(E_{1}^{2}-E_{3}^{2}\right)}\left\\{Y_{11}\left[\frac{1}{E_{1}}\left(\frac{1}{i\hbar\omega_{m}-E_{1}}-\frac{1}{i\hbar\omega_{m}-E_{2}}\right)-\frac{1}{E_{3}}\left(\frac{1}{i\hbar\omega_{m}-E_{3}}-\frac{1}{i\hbar\omega_{m}-E_{4}}\right)\right]\right.$
$\displaystyle\left.-Y_{10}\left[E_{1}\left(\frac{1}{i\hbar\omega_{m}-E_{1}}-\frac{1}{i\hbar\omega_{m}-E_{2}}\right)-E_{3}\left(\frac{1}{i\hbar\omega_{m}-E_{3}}-\frac{1}{i\hbar\omega_{m}-E_{4}}\right)\right]\right\\}$
$\displaystyle
Y_{10}=-\frac{3\alpha^{2}\lambda_{R}}{2\langle\Psi|\Psi\rangle},\;Y_{11}=Y_{10}\left[(\hbar
v)^{2}(1-3\alpha^{2})^{2}(k_{x}^{2}-k_{y}^{2})+\left(\frac{3\alpha^{2}\lambda_{I}}{\langle\Psi|\Psi\rangle}\right)^{2}\right]$
$\displaystyle
B_{y1}(\mathbf{k},i\omega_{m})=-Y_{7}\frac{1}{2\left(E_{3}^{3}-E_{1}^{2}\right)}\left\\{\frac{1}{i\hbar\omega_{m}-E_{1}}+\frac{1}{i\hbar\omega_{m}-E_{2}}-\frac{1}{i\hbar\omega_{m}-E_{3}}-\frac{1}{i\hbar\omega_{m}-E_{4}}\right\\}k_{x}$
(S.37)

$\displaystyle
B_{yz}(\mathbf{k},i\omega_{m})=-Y_{8}\frac{1}{2\left(E_{1}^{2}-E_{3}^{2}\right)}\,k_{x}\left\\{\frac{1}{E_{1}}\left[\frac{1}{i\hbar\omega_{m}-E_{1}}-\frac{1}{i\hbar\omega_{m}-E_{2}}\right]-\frac{1}{E_{3}}\left[\frac{1}{i\hbar\omega_{m}-E_{3}}-\frac{1}{i\hbar\omega_{m}-E_{4}}\right]\right\\}$
(S.38) $\displaystyle D_{yx}(\mathbf{k},i\omega_{m})$ $\displaystyle=$
$\displaystyle\xi\frac{1}{2\left(E_{1}^{2}-E_{3}^{2}\right)}\left\\{Y_{12}\left[\frac{1}{E_{1}}\left(\frac{1}{i\hbar\omega_{m}-E_{1}}-\frac{1}{i\hbar\omega_{m}-E_{2}}\right)-\frac{1}{E_{3}}\left(\frac{1}{i\hbar\omega_{m}-E_{3}}-\frac{1}{i\hbar\omega_{m}-E_{4}}\right)\right]\right.$
$\displaystyle\left.+Y_{10}\left[E_{1}\left(\frac{1}{i\hbar\omega_{m}-E_{1}}-\frac{1}{i\hbar\omega_{m}-E_{2}}\right)-E_{3}\left(\frac{1}{i\hbar\omega_{m}-E_{3}}-\frac{1}{i\hbar\omega_{m}-E_{4}}\right)\right]\right\\}$
$\displaystyle
Y_{12}=\frac{3\alpha^{2}\lambda_{R}}{3\langle\Psi|\Psi\rangle}\left[-(\hbar
v)^{2}(1-3\alpha^{2})^{2}(k_{x}^{2}-k_{y}^{2})+\left(\frac{3\alpha^{2}\lambda_{I}}{\langle\Psi|\Psi\rangle}\right)^{2}\right]$
$\displaystyle D_{yy}(\mathbf{k},i\omega_{m})$
$\displaystyle=-Y_{9}\frac{1}{\left(E_{1}^{2}-E_{3}^{2}\right)}k_{x}k_{y}\left\\{\frac{1}{E_{1}}\left[\frac{1}{i\hbar\omega_{m}-E_{1}}-\frac{1}{i\hbar\omega_{m}+E_{1}}\right]-\frac{1}{E_{3}}\left[\frac{1}{i\hbar\omega_{m}-E_{3}}-\frac{1}{i\hbar\omega_{m}-E_{4}}\right]\right\\}$
(S.40)

##  S.4 III. Magnon Green function and Gilbert damping

###  S.4.1 Interfacial exchange coupling between a ferro- magnetic insulator
(FI) and a TBG

We consider the Hamiltonian of the ferromagnetic insulator (FI) in the
independent magnon approximation as

$\displaystyle H_{\rm
FI}=\sum_{\bm{k}}\hbar\omega_{\bm{k}}b_{\bm{k}}^{\dagger}b_{\bm{k}},$ (S.41)

where $\hbar\omega_{\bm{k}}\simeq{\cal D}|{\bm{k}}|^{2}+\hbar\gamma h_{\rm
dc}$ is a dispersion of magnons, ${\cal D}$ is a spin stiffness, $\gamma$ is
the gyromagnetic ratio, $h_{\rm dc}$ is a static magnetic field. In the spin
pumping setup, only the static part associated to $\bm{k=0}$ is relevant.
Considering only the uniform spin precession, the Hamiltonian of the FI can be
simply written as

$\displaystyle H_{\rm FI}=\hbar\omega_{\bm{0}}b_{\bm{0}}^{\dagger}b_{\bm{0}},$
(S.42)

where $b_{\bm{q}}$ is the Fourier transformation of the site representation
$b_{i}$ defined as

$\displaystyle b_{i}$ $\displaystyle=\frac{1}{\sqrt{N_{\rm
FI}}}\sum_{\bm{q}}e^{i{\bm{q}}\cdot{\bm{r}}_{i}}b_{\bm{q}}\simeq\frac{1}{\sqrt{N_{\rm
FI}}}b_{\bm{0}},$ (S.43) $\displaystyle b_{i}^{\dagger}$
$\displaystyle=\frac{1}{\sqrt{N_{\rm
FI}}}\sum_{\bm{q}}e^{-i{\bm{q}}\cdot{\bm{r}}_{i}}b^{\dagger}_{\bm{q}}\simeq\frac{1}{\sqrt{N_{\rm
FI}}}b^{\dagger}_{\bm{0}},$ (S.44)

where $N_{\rm FI}$ is the number of unit cells in the FI.

We consider the (retarded) magnon Green function as

$\displaystyle G^{R}({\bm{q}},\omega)$ $\displaystyle=\int
dt\,G^{R}({\bm{q}},t)e^{i\omega t},$ (S.45) $\displaystyle G^{R}({\bm{q}},t)$
$\displaystyle=-\frac{i}{\hbar}\theta(t)\langle[S^{+}_{\bm{q}}(t),S^{-}_{\bm{q}}(0)]\rangle$
$\displaystyle=-\frac{2iS_{0}}{\hbar}\theta(t)\langle[b_{\bm{q}}(t),b_{\bm{q}}^{\dagger}(0)]\rangle.$
(S.46)

In the absence of the junction, the magnon Green function is

$\displaystyle
G_{0}^{R}({\bm{q}},\omega)=\frac{2S_{0}/\hbar}{\omega-\omega_{\bm{q}}+i\delta}.$
(S.47)

Usually, we introduce spin relaxation of bulk FI as

$\displaystyle
G_{0}^{R}({\bm{q}},\omega)=\frac{2S_{0}/\hbar}{\omega-\omega_{\bm{q}}+i\alpha_{\rm
G}\omega},$ (S.48)

where $\alpha_{\rm G}$ is a dimensionless strength of the Gilbert damping,
which is of order of $10^{-4}$–$10^{-3}$. We note that the FMR line shape is
proportional to ${\rm Im}\,G_{0}^{R}({\bm{q}}={\bm{0}},\omega)$.
In the presence of the interfacial coupling and for a uniform spin precession,
the magnon Green function is given by the Dyson equation as

$\displaystyle
G_{0}^{R}({\bm{q=0}},\omega)=\frac{2S_{0}/\hbar}{\omega-\omega_{\bm{q}=0}+i\alpha_{\rm
G}\omega-(2S_{0}/\hbar)\Sigma^{R}(\bm{q=0},\omega)},$ (S.49)

where $\Sigma^{R}(\omega)$ is the self-energy. Although the real part of
$\Sigma^{R}(\omega)$ is related to the shift of the ferromagnetic resonance,
we neglect it for simplicity. Then, the magnon Green function is rewritten as

$\displaystyle G_{0}^{R}({\bm{0}},\omega)$
$\displaystyle=\frac{2S_{0}/\hbar}{\omega-\omega_{\bm{q}=0}+i(\alpha_{\rm
G}+\delta\alpha_{\rm G})}\omega,$ (S.50) $\displaystyle\delta\alpha_{\rm
G}(\omega)$
$\displaystyle=-\frac{2S_{0}}{\hbar\omega}\rm{Im}\Sigma^{R}(\bm{q=0},\omega).$
(S.51)

We note that $\delta\alpha_{\rm G}(\omega)$ depends on $\omega$ in general.
However, since the ferromagnetic resonance peak is sharp enough ($\alpha_{\rm
G}+\delta\alpha_{\rm G}\ll 1$), we can replace $\omega$ with $\omega_{\rm
0}=\Omega$ (the peak position of the ferromagnetic resonance):

$\displaystyle\delta\alpha_{\rm G}$
$\displaystyle\simeq-\frac{2S_{0}}{\hbar\Omega}\Sigma^{R}(\bm{q=0},\Omega).$
(S.52)

The Hamiltonian of the interfacial coupling is given as

$\displaystyle H_{\rm int}=\sum_{\langle
i,j\rangle}T_{ij}(S^{+}_{i}s^{-}_{j}+{\rm h.c.}).$ (S.53)

Here, $S^{\pm}_{i}$ is a spin ladder operator of the FI and is described by
magnon annihilation/creation operators ($b_{i}$ and $b_{i}^{\dagger}$) as

$\displaystyle S^{+}_{i}=\sqrt{2S_{0}}b_{i},\quad
S^{-}_{i}=\sqrt{2S_{0}}b_{i}^{\dagger},$ (S.54)

where $S_{0}$ is an amplitude of the localized spin in the FI. $s^{\pm}_{j}$
is a spin ladder operator of electrons in two-dimensional electron systems
(twisted bilayer graphene) and is described by the electron
annihilation/creation operators ($c_{j\sigma}$ and $c_{j\sigma}^{\dagger}$) as

$\displaystyle s^{+}_{j}=c_{j\uparrow}^{\dagger}c_{j\downarrow},\quad
s^{-}_{j}=c_{j\downarrow}^{\dagger}c_{j\uparrow}.$ (S.55)

We define the Fourier transformation as

$\displaystyle c_{j\sigma}$
$\displaystyle=\frac{1}{\sqrt{N}}\sum_{\bm{k}}e^{i{\bm{k}}\cdot{\bm{r}}_{j}}c_{{\bm{k}}\sigma},$
(S.56) $\displaystyle c^{\dagger}_{j\sigma}$
$\displaystyle=\frac{1}{\sqrt{N}}\sum_{\bm{k}}e^{-i{\bm{k}}\cdot{\bm{r}}_{j}}c^{\dagger}_{{\bm{k}}\sigma},$
(S.57)

where $N$ is the number of unit cells and ${\bm{r}}_{j}$ is the position of
the site $j$ in TBG. Then, we obtain

$\displaystyle s^{+}_{j}$
$\displaystyle=\frac{1}{N}\sum_{{\bm{k}},{\bm{k}}^{\prime}}e^{-i{\bm{k}}\cdot{\bm{r}}_{j}+i{\bm{k}}^{\prime}\cdot{\bm{r}}_{j}}c_{{\bm{k}}\uparrow}^{\dagger}c_{{\bm{k}}^{\prime}\downarrow},$
(S.58) $\displaystyle s^{-}_{j}$
$\displaystyle=\frac{1}{N}\sum_{{\bm{k}},{\bm{k}}^{\prime}}e^{-i{\bm{k}}\cdot{\bm{r}}_{j}+i{\bm{k}}^{\prime}\cdot{\bm{r}}_{j}}c_{{\bm{k}}\downarrow}^{\dagger}c_{{\bm{k}}^{\prime}\uparrow}.$
(S.59)

We define the Fourier transformation of $s^{\pm}_{j}$ as

$\displaystyle s^{+}_{j}$
$\displaystyle=\frac{1}{N}\sum_{\bm{q}}e^{i{\bm{q}}\cdot{\bm{r}}_{j}}s^{+}_{\bm{q}},$
(S.60) $\displaystyle s^{-}_{j}$
$\displaystyle=\frac{1}{N}\sum_{\bm{q}}e^{i{\bm{q}}\cdot{\bm{r}}_{j}}s^{-}_{\bm{q}}.$
(S.61)

From $s^{-}_{j}=(s^{+}_{j})^{\dagger}$, we obtain the relation

$\displaystyle s^{-}_{\bm{q}}=(s^{+}_{-{\bm{q}}})^{\dagger}.$ (S.62)

Note that $(s^{-}_{\bm{q}})^{\dagger}\neq s^{+}_{{\bm{q}}}$. The inverse
Fourier transformation is given as

$\displaystyle s^{+}_{\bm{q}}$
$\displaystyle=\sum_{j}e^{-i{\bm{q}}\cdot{\bm{r}}_{j}}s^{+}_{j},$ (S.63)
$\displaystyle s^{-}_{\bm{q}}$
$\displaystyle=(s^{+}_{-{\bm{q}}})=\sum_{j}e^{-i{\bm{q}}\cdot{\bm{r}}_{j}}s^{-}_{j}.$
(S.64)

Especially for ${\bm{q}}={\bm{0}}$, we obtain

$\displaystyle s^{+}_{\bm{0}}=\sum_{j}s^{+}_{j},\quad
s^{-}_{\bm{0}}=\sum_{j}s^{-}_{j},$ (S.65)

Using Eqs. (S.56) and (S.57), we obtain

$\displaystyle s^{+}_{\bm{q}}$
$\displaystyle=\frac{1}{N}\sum_{j}e^{-i{\bm{q}}\cdot{\bm{r}}_{j}}\sum_{{\bm{k}},{\bm{k}}^{\prime}}e^{-i{\bm{k}}\cdot{\bm{r}}_{j}+i{\bm{k}}^{\prime}\cdot{\bm{r}}_{j}}c_{{\bm{k}}\uparrow}^{\dagger}c_{{\bm{k}}^{\prime}\downarrow}$
$\displaystyle=\sum_{{\bm{k}}}c_{{\bm{k}}\uparrow}^{\dagger}c_{{\bm{k}}+{\bm{q}}\downarrow},$
(S.66) $\displaystyle s^{-}_{\bm{q}}$
$\displaystyle=\frac{1}{N}\sum_{j}e^{-i{\bm{q}}\cdot{\bm{r}}_{j}}\sum_{{\bm{k}},{\bm{k}}^{\prime}}e^{-i{\bm{k}}\cdot{\bm{r}}_{j}+i{\bm{k}}^{\prime}\cdot{\bm{r}}_{j}}c_{{\bm{k}}\downarrow}^{\dagger}c_{{\bm{k}}^{\prime}\uparrow}$
$\displaystyle=\sum_{{\bm{k}}}c_{{\bm{k}}\downarrow}^{\dagger}c_{{\bm{k}}+{\bm{q}}\uparrow}.$
(S.67)

For a clean interface, we can set $T_{ij}=T$. Then, using Eqs. (S.43) and
(S.44), the Hamiltonian of the interface is written as

$\displaystyle H_{\rm int}$ $\displaystyle=\frac{T\sqrt{2S_{0}}}{\sqrt{N_{\rm
FI}}}\sum_{\langle
i,j\rangle}(b_{\bm{0}}s_{j}^{-}+b_{\bm{0}}^{\dagger}s_{j}^{+})$
$\displaystyle\simeq\frac{T\sqrt{2S_{0}}}{\sqrt{N_{\rm
FI}}}\left[b_{\bm{0}}\Bigl{(}\sum_{j}s_{j}^{-}\Bigr{)}+b_{\bm{0}}^{\dagger}\Bigl{(}\sum_{j}s_{j}^{+}\Bigr{)}\right]$
$\displaystyle=\sqrt{2S_{0}}b_{\bm{0}}\tilde{s}^{-}+\sqrt{2S_{0}}b_{\bm{0}}\tilde{s}^{+},$
(S.68)

where $\tilde{s}^{\pm}$ is defined as

$\displaystyle\tilde{s}^{-}$ $\displaystyle=\frac{T}{\sqrt{N_{\rm
FI}}}s_{\bm{0}}^{-},$ (S.69) $\displaystyle\tilde{s}^{+}$
$\displaystyle=\frac{T}{\sqrt{N_{\rm FI}}}s_{\bm{0}}^{+}.$ (S.70)

By the second-order perturbation, the self-energy of the magnon at
${\bm{q}}={\bm{0}}$ is calculated as

$\displaystyle\Sigma^{R}(\omega)$ $\displaystyle=\int
dt\,\Sigma^{R}(t)e^{i\omega t},$ (S.71) $\displaystyle\Sigma^{R}(t)$
$\displaystyle=-\frac{i}{\hbar}\theta(t)\langle[\tilde{s}^{+}(t),\tilde{s}^{-}(0)]\rangle.$
(S.72)

The self-energy can be related to a retarded component of the dynamic spin
susceptibility per unit area as

$\displaystyle\Sigma^{R}(\omega)$ $\displaystyle=-\frac{T^{2}N}{N_{\rm
FI}}\chi({\bm{0}},\omega),$ (S.73) $\displaystyle\chi({\bm{q}},\omega)$
$\displaystyle=\int dt\,\chi({\bm{q}},t)e^{i\omega t},$ (S.74)
$\displaystyle\chi({\bm{q}},t)$
$\displaystyle=\frac{i}{N\hbar}\theta(t)\langle[s_{-{\bm{q}}}^{+}(t),s_{\bm{q}}^{-}(0)]\rangle.$
(S.75)

$\chi({\bm{q}},t)$ is calculated for one-band of TBG without spin-orbit
interaction as

$\displaystyle\chi({\bm{q}},t)=\frac{1}{N}\sum_{\bm{q}}\frac{f(\xi_{\bm{k}})-f(\xi_{{\bm{k}}+{\bm{q}}})}{\hbar\omega+i\delta+\xi_{{\bm{k}}+{\bm{q}}}-\xi_{\bm{k}}},$
(S.76)

where $\xi_{\bm{k}}=\epsilon_{\bm{k}}-\mu$, $\epsilon_{\bm{k}}$ is a
dispersion of electrons, $\mu$ is a chemical potential. This is just a
Lindhard function.
We note that $\chi({\bm{q}},t)$ is independent of the system size (area). For
systems with spin-orbit interaction, we have to extend the Lindhard function
into the spin-dependent one.

Then, the enhancement of the Gilbert damping is written as

$\displaystyle\delta\alpha_{\rm G}$
$\displaystyle=-\frac{2S_{0}}{\hbar\Omega}\,{\rm
Im}\Sigma^{R}(\rm{q=0},\Omega)$ $\displaystyle=\frac{2S_{0}T^{2}N}{N_{\rm
FI}\hbar\omega_{\bm{0}}}\chi(\rm{q=0},\Omega).$ (S.77)

We note that the number of the unit cell of twisted bilayer graphene is
written as $N=S/A$ where $S$ is a area of the junction and $A$ is an area of a
unit cell of twisted bilayer graphene. We also note that the number of the
unit cell of the FI is written as $N_{\rm FI}=Sd/a^{3}$ where $d$ is a
thickness of the FI, $a$ is a lattice constant of the FI. Using these
parameters, we obtain

$\displaystyle\delta\alpha_{\rm G}$
$\displaystyle=\frac{2S_{0}T^{2}a^{3}}{Ad\hbar\omega_{\bm{0}}}\chi(\rm{q=0},\Omega).$
(S.78)

We note that $\delta\alpha_{\rm G}$ is proportional to $1/d$ in consistent
with experimental results.

Table 1: Experimental parameters. Microwave frequency | $\omega_{\bm{0}}$ | $1\,{\rm GHz}$
---|---|---
Amplitude of spins of FI | $S_{0}$ | 10
Lattice constant of FI | $a$ |  $12.376\,$Å
Thickness of FI | $d$ | $\geq 10\,{\rm nm}$
Interfacial exchange coupling | $J$ |  $\sim 1\,{\rm K}$ (not known)

If YIG (Yttrium Iron Garnet) is chosen as the ferromagnet insulator, the
parameter is given in the Table.

###  S.4.2 Electronic spins in the FI magnetization frame

Regarding the dependence of the electronic Hamiltonian (Eq. 5 of the main
text) on $s_{z}$, one should consider a $3D$ FI magnetization as in Ref.
Funato . The average spin vector is along the orthoradial spherical vector
$\langle\mathbf{S}_{FI}\rangle=\langle\mathbf{S}_{FI}\rangle\,\mathbf{u}_{x^{\prime}}$.
The radial vector $\mathbf{u}_{z^{\prime}}$ forms an angle $\theta_{m}$ with
the $z$ axis perpendicular to the interface. The third axis $y^{\prime}$ is in
the $(xoy)$ plane and its unit vector is the orthoradial inplane vector
$\mathbf{u}_{y^{\prime}}=-\sin\Phi_{m}\;\mathbf{u}_{x}+\cos\Phi_{m}\;\mathbf{u}_{y}$
as shown in Fig.S.2.

Figure S.2: magnetization-fixed coordinate frame
$\left(x^{\prime},y^{\prime},z^{\prime}\right)$ with respect to the Laboratory
frame $\left(x,y,z\right)$.

In the FI spin frame ($x^{\prime},y^{\prime},z^{\prime}$), the components of
the electronic spin operators are given by:

$\displaystyle
s^{x^{\prime}}_{\mathbf{k}}=\mathbf{s}_{\mathbf{k}}\cdot\mathbf{u}_{x^{\prime}}=\cos\theta_{m}\cos\Phi_{m}\,s^{x}_{\mathbf{k}}+\cos\theta_{m}\sin\Phi_{m}\,s^{y}_{\mathbf{k}}-\sin\theta_{m}\,s^{z}_{\mathbf{k}}$
$\displaystyle
s^{y^{\prime}}_{\mathbf{k}}=\mathbf{s}_{\mathbf{k}}\cdot\mathbf{u}_{y^{\prime}}=-\sin\Phi_{m}\,s^{x}_{\mathbf{k}}+\cos\Phi_{m}\,s^{y}_{\mathbf{k}}$
$\displaystyle
s^{z^{\prime}}_{\mathbf{k}}=\mathbf{s}_{\mathbf{k}}\cdot\mathbf{u}_{z^{\prime}}=\sin\theta_{m}\cos\Phi_{m}\,s^{x}_{\mathbf{k}}+\sin\theta_{m}\sin\Phi_{m}\,s^{y}_{\mathbf{k}}+\cos\theta_{m}\,s^{z}_{\mathbf{k}}$

We define the ladder electronic spin operators as

$\displaystyle s^{x^{\prime},\pm}_{\mathbf{k}}$ $\displaystyle=$
$\displaystyle s^{y^{\prime}}_{\pm\mathbf{k}}\pm
is^{z^{\prime}}_{\pm\mathbf{k}}=\frac{1}{2}\sum_{\sigma,\sigma^{\prime},\mathbf{k}^{\prime}}c^{\dagger}_{\mathbf{k}^{\prime},\sigma}\left(\sigma_{s}^{x^{\prime},\pm}\right)_{\sigma,\sigma^{\prime}}c^{\dagger}_{\mathbf{k}^{\prime}\pm\mathbf{k},\sigma^{\prime}}$
(S.80)

where $s^{x^{\prime},\pm}=s_{x}\left(-\sin\Phi_{m}\pm
i\sin\theta_{m}\cos\Phi_{m}\right)+s_{y}\left(\cos\Phi_{m}\pm
i\sin\theta_{m}\sin\Phi_{m}\right)\pm i\cos\theta_{m}s_{z}$.

###  S.4.3 Magnon self-energy

In the second order perturbation with respect to the interfacial exchange
interaction $T_{\mathbf{q}}$ the magnon self-energy is given by Yama

$\displaystyle\Sigma(\mathbf{q},i\omega_{n})=\frac{|T_{\mathbf{q}}|^{2}}{4\beta}\sum_{\mathbf{k},i\omega_{m}}\text{Tr}\left[\sigma_{s}^{x^{\prime},-}\;\hat{g}(\mathbf{k},i\omega_{m})\sigma_{s}^{x^{\prime},+}\;\hat{g}(\mathbf{k}+\mathbf{q},i\omega_{m}+i\omega_{n})\right]$
(S.81)

where $\hat{g}(\mathbf{k},i\omega_{m})$ is the electronic Green function given
by Eq. S.22.

The trace term is of the form:

$\displaystyle\text{Tr}\left[\mathbf{a}^{\ast}\cdot\mathbf{\sigma}_{s}\left(\hat{g}_{0}+\hat{\mathbf{g}}\cdot\mathbf{\sigma}_{s}\right)\mathbf{a}\cdot\mathbf{\sigma}_{s}\left(\hat{g}^{\prime}_{0}+\hat{\mathbf{g}}^{\prime}\cdot\mathbf{\sigma}_{s}\right)\right]$
(S.82)

where the vector
$\mathbf{a}=\left(-\sin\Phi_{m}+i\sin\theta_{m}\cos\Phi_{m},\cos\Phi_{m}+i\sin\theta_{m}\sin\Phi_{m},i\cos\theta_{m}\right)$
is written in the laboratory frame $(x,y,z)$.

We set $\hat{g}=\hat{g}(\mathbf{k},i\omega_{m})$ and
$\hat{g}^{\prime}=\hat{g}(\mathbf{k},i\omega_{m}+i\omega_{n})$. Taking into
account the operator character of $\hat{g}$ one could use the identity

$\displaystyle\left(\mathbf{a}\cdot\mathbf{\sigma}_{s}\right)\left(\mathbf{b}\cdot\mathbf{\sigma}_{s}\right)=\left(\mathbf{a}\cdot\mathbf{b}\right)\mathbb{I}+i\left(\mathbf{a}\times\mathbf{b}\right)\cdot\mathbf{\sigma}_{s}$
(S.83)

Given the expressions of $\hat{g}$ and $\hat{g}^{\prime}$ in Eq. S.23, the
trace term (Eq. S.82) reduces to:

$\displaystyle\text{Tr}\left[\mathbf{a}^{\ast}\cdot\mathbf{\sigma}_{s}\left(\hat{g}_{0}+\hat{\mathbf{g}}\cdot\mathbf{\sigma}_{s}\right)\mathbf{a}\cdot\mathbf{\sigma}_{s}\left(\hat{g}^{\prime}_{0}+\hat{\mathbf{g}}^{\prime}\cdot\mathbf{\sigma}_{s}\right)\right]=\sum_{i=0,1,2}F_{i}(\mathbf{k},i\omega_{m},i\omega_{n})$
(S.84)

where

$\displaystyle
F_{0}(\mathbf{k},i\omega_{m},i\omega_{n})=4\left(A_{01}A^{\prime}_{01}+A_{0z}A^{\prime}_{0z}+C_{0x}C^{\prime}_{0x}+C_{0y}C^{\prime}_{0y}\right)$
$\displaystyle
F_{1}(\mathbf{k},i\omega_{m},i\omega_{n})=-2\left\\{\cos\theta_{m}\cos\Phi_{m}\right.$
$\displaystyle\left(A_{01}B^{\prime}_{x1}-B_{x1}A^{\prime}_{01}+A_{0z}B^{\prime}_{xz}-A^{\prime}_{0z}B_{xz}+C_{0x}D^{\prime}_{xx}-C^{\prime}_{0x}D_{xx}+C_{0y}D^{\prime}_{xy}-C^{\prime}_{0y}D_{xy}\right)$
$\displaystyle+\cos\theta_{m}\sin\Phi_{m}\left(A_{01}B^{\prime}_{y1}-B_{y1}A^{\prime}_{01}+A_{0z}B^{\prime}_{yz}-A^{\prime}_{0z}B_{yz}+C_{0x}D^{\prime}_{yx}-C^{\prime}_{0x}D_{yx}+C_{0y}D^{\prime}_{yy}-C^{\prime}_{0y}D_{yy}\right)$
$\displaystyle\left.-\sin\theta_{m}\left(A_{01}A^{\prime}_{z1}-A^{\prime}_{01}A_{z1}+A_{0z}A^{\prime}_{zz}-A^{\prime}_{0z}A_{zz}+C_{0x}C^{\prime}_{zx}-C^{\prime}_{0x}C_{zx}+C_{0y}C^{\prime}_{zy}-C^{\prime}_{zy}C_{0y}\right)\right\\}$
$\displaystyle
F_{2}(\mathbf{k},i\omega_{m},i\omega_{n})=-2\cos^{2}\theta_{m}\cos^{2}\Phi_{m}\left(B_{x1}B^{\prime}_{x1}+B_{xz}B^{\prime}_{xz}+D_{xx}D^{\prime}_{xx}+D_{xy}D^{\prime}_{xy}\right)$
$\displaystyle-2\cos^{2}\theta_{m}\sin^{2}\Phi_{m}\left(B_{y1}B^{\prime}_{y1}+B_{yz}B^{\prime}_{yz}+D_{yx}D^{\prime}_{yx}+D_{yy}D^{\prime}_{yy}\right)$
$\displaystyle-2\sin^{2}\theta_{m}\left(A_{z1}A^{\prime}_{z1}+A_{zz}A^{\prime}_{zz}+C_{zx}C^{\prime}_{zx}+C_{zy}C^{\prime}_{zy}\right)$
$\displaystyle-\cos^{2}\theta_{m}\sin
2\Phi_{m}\left(B_{x1}B^{\prime}_{y1}+B_{xz}B^{\prime}_{yz}+D_{xx}D^{\prime}_{yx}+D_{xy}D^{\prime}_{yy}+B^{\prime}_{x1}B_{y1}+B^{\prime}_{xz}B_{yz}+D^{\prime}_{xx}D_{yx}+D^{\prime}_{xy}D_{yy}\right)$
$\displaystyle+\cos\Phi_{m}\sin
2\theta_{m}\left(B_{x1}A^{\prime}_{z1}+B_{xz}A^{\prime}_{zz}+D_{xx}C^{\prime}_{zx}+D_{xy}C^{\prime}_{zy}B^{\prime}_{x1}A_{z1}+B^{\prime}_{xz}A_{zz}+D^{\prime}_{xx}C_{zx}+D^{\prime}_{xy}C_{zy}\right)$
$\displaystyle+\sin\Phi_{m}\sin
2\theta_{m}\left(B_{y1}A^{\prime}_{z1}+B_{yz}A^{\prime}_{zz}+D_{yx}C^{\prime}_{zx}+D_{yy}C^{\prime}_{zy}B^{\prime}_{y1}A_{z1}+B^{\prime}_{yz}A_{zz}+D^{\prime}_{yx}C_{zx}+D^{\prime}_{yy}C_{zy}\right)$
(S.85)

The terms with a prime are expressed in terms of
$i\omega_{n}^{\prime}=i\omega_{n}+i\omega_{m}$.

Regarding the $\mathbf{k}$ dependence of the $A,B,C$ and $D$ operators (Eqs.
S.3\- S.40), only $F_{0}$, the last term in $F_{1}$ and the three first terms
in $F_{2}$ give non-vanishing contributions after summing over $\mathbf{k}$ in
Eq. 12.

On the other hand, the terms between parentheses in the first and second line
in $F_{2}$ expression give the same contribution. As a result, the GD is found
to be independent of the azimutal angle $\Phi_{m}$, which expresses isotropy
of the electronic band structure $E_{\sigma,\mp}$ (Eq.S.18). However, the GD
depends on the out-of-plane orientation of the FI magnetization via the angle
$\theta_{m}$.

According to Eq. S.81, the terms to calculate are of the form

$\displaystyle\sum_{\mathbf{k}}F(\mathbf{k})\sum_{\omega_{m}}\frac{1}{i\hbar\omega_{m}-E_{i}}\frac{1}{i\hbar\omega_{m}^{\prime}-E_{j}},$
(S.86)

where $\omega_{m}^{\prime}=\omega_{m}+\omega_{n}$ and $F(\mathbf{k})$ is a
function of $\mathbf{k}=\left(k,\varphi_{\mathbf{k}}\right)$.

The summation over $\omega_{m}$ in Eq. S.86 can be written as

$\displaystyle\sum_{\omega_{m}}\frac{1}{i\hbar\omega_{m}-E_{i}}\frac{1}{i\hbar\omega_{m}^{\prime}-E_{j}}$
$\displaystyle=\frac{1}{i\hbar\omega_{n}-(E_{j}-E_{i})}\sum_{\omega_{m}}\left[\frac{1}{i\hbar\omega_{m}-E_{i}}-\frac{1}{i\hbar\omega_{m}^{\prime}-E_{j}}\right]$
(S.87)
$\displaystyle=-\frac{1}{i\hbar\omega_{n}-(E_{j}-E_{i})}\frac{1}{k_{B}T}\int_{c}\frac{dz}{2\pi
i}h(z)f(z)$

where $h(z)=\frac{1}{z-E_{i}}-\frac{1}{z+i\hbar\omega_{n}-E_{j}}$, $f(z)$ is
the Fermi-Dirac function, $C$ is clockwise contour around the poles
$z=i\hbar\omega_{m}$.

Equation S.86 reduces, then, to

$\displaystyle\sum_{\omega_{m}}\frac{1}{i\hbar\omega_{m}-E_{i}}\frac{1}{i\hbar\omega_{m}^{\prime}-E_{j}}=\frac{f(E_{j})-f(E_{i})}{i\hbar\omega_{n}-(E_{j}-E_{i})}$
(S.88)

Taking the analytic continuation $i\hbar\omega_{n}=\hbar\omega+i\eta$, Eq.
S.86 becomes

$\displaystyle\lim_{\eta\to
0}\sum_{\mathbf{k}}F(\mathbf{k})\eta\;\frac{f(E_{j})-f(E_{i})}{(\hbar\omega-
E_{j}+E_{i})^{2}+\eta^{2}}=\lim_{\eta\to
0}\sum_{\mathbf{k}}F(\mathbf{k})(f(E_{j})-f(E_{i}))L(\hbar\omega-(E_{j}-E_{i})),$
(S.89)

$L(x)=\frac{\eta}{x^{2}+\eta^{2}}$ being the Lorentzian function. The sum over
$\mathbf{k}=(k,\varphi_{\mathbf{k}})$ in Eq. S.86 reduces to
$\frac{A}{(2\pi)^{2}}\int_{0}^{k_{c}}kdk\int_{0}^{2\pi}d\varphi_{\mathbf{k}}$,
where $A$ is the moiré superlattice area, $k_{c}$ is a cutoff on the momentum
amplitude $k$, below which the continuum model for the monolayer is justified.
We take $k_{c}=q_{0}/2$, where $q_{0}=\frac{4\pi}{3a}\theta$ is the separation
between the Dirac points $\mathbf{K}_{1,\xi}$ and $\mathbf{K}_{2,\xi}$ of,
respectively, layer (1) and layer (2) at a given monolayer valley $\xi$.

###  S.4.4 Gilbert damping correction

For a uniform spin precession, the Gilbert damping correction
$\delta\alpha_{G}$, at the FMR frequency $\Omega$, can be expressed as Yama

$\displaystyle\delta\alpha_{G}=-\frac{2S_{0}}{\hbar\Omega}\mathrm{Im}\Sigma(\mathbf{q=0},\Omega)$
(S.90)

The imaginary part of the magnon self-energy is of the form
$\mathrm{Im}\Sigma(\mathbf{q=0},\Omega)=\frac{T_{0}^{2}}{\hbar\Omega}\tilde{\Sigma}(\mathbf{q=0},\Omega)$,
where $\tilde{\Sigma}(\mathbf{q=0},\Omega)$ is a dimensionless integral.
Introducing the average SOC $\lambda=\frac{1}{2}(\lambda_{I}+\lambda_{R})$,
$\delta\alpha_{G}$ can be written as

$\displaystyle\delta\alpha_{G}=-\alpha_{G}^{0}\left(\frac{\lambda}{\hbar\Omega}\right)^{2}\tilde{\Sigma}(\mathbf{q=0},\Omega)$
(S.91)

where $\alpha_{G}^{0}=2S_{0}\frac{T_{0}^{2}}{\lambda^{2}}$.

In Fig. S.3, we plot the normalized Gilbert damping correction
$\delta\alpha_{G}/\alpha_{G}^{0}$ as a function of the twist angle $\theta$
and the FMR energy $\hbar\Omega$ at different temperatures. The bottom panels
are a zoomed representation around the MA. Fig. S.3 shows that, the GD
increases by decreasing the twist angle and sharply drops at the MA,
regardless of the temperature range and the FMR frequency .

At high temperature ($k_{B}T>\lambda$), the GD exhibits, around the MA, a fine
structure characterized by a peak which disappears at low temperature. The
origin of this peak is, as discussed in the main text, due to the dispersion
of the energy bands of TBG with the SOC and their corresponding thermal
weights.

Figure S.3: Color plot of the normalized Gilbert damping correction
$\delta\alpha_{G}/\alpha_{G}^{0}$ as a function of the twist angle $\theta$
and the FMR energy $\hbar\Omega$ at $k_{B}T=0.1\;\mathrm{meV}$ ((a) and (d)),
$k_{B}T=1\;\mathrm{meV}$ ((b) and (e)) and $k_{B}T=25\;\mathrm{meV}$ ((c) and
(f)). The bottom panels show the behavior of the GD around the MA.
Calculations are done for $\mu=0$, $\lambda_{I}=3\;\mathrm{meV}$ and
$\lambda_{R}=4\,\mathrm{meV}$.

According to Eq.S.89, the GD depends on the thermal weight $\Delta
f(E)=f(E_{\langle S_{z}\rangle})-f(E_{-\langle S_{z}\rangle})$ (see main
text). In Fig. S.4 we plot $\Delta f(E)$ corresponding to the transitions
between $E_{-,+}\rightarrow E_{+,+}$ and $E_{-,-}\rightarrow E_{+,-}$ in the
case of the undoped system.

Figures S.4 (a) and (b) show that, at high temperature
($k_{B}T>\lambda_{I},\lambda_{R}$), $\Delta f(E)$ increases as the band
dispersion gets larger and reaches its minimal value at the MA. This behavior
explains the drop of the GD at the MA and its enhancement at small twist
angles.

In figure S.4 (c), we plot $\Delta f(E)$ around the MA, at relatively high
thermal energy compared to the SOC, where the GD exhibits a peak at the MA
(Fig. 2 of the main text). In this case, $\Delta f(E)$ is maximal at the MA
and decreases at the angles $\theta^{+}_{M}$ and $\theta^{-}_{M}$ close to the
MA. This feature results from the decrease of the energy separation between
$E_{-,-}$ and $E_{+,-}$ at $\theta^{+}_{M}$ and $\theta^{-}_{M}$, compared to
that at $\theta_{M}$ (Fig. 3 of the main text). At low temperature, and around
the MA, one gets $\Delta f(E)=1$ for the transitions $E_{-,-}\rightarrow
E_{+,-}$ and $E_{-,+}\rightarrow E_{+,+}$. As a consequence, the GD behavior
is now only dependent on the effective Fermi velocity $v^{ast}$ which vanishes
at the MA. As a consequence, the small peak of the GD, emerging at the MA at
relatively high temperature, disappears.

$\begin{array}[]{ccc}\includegraphics[width=65.04366pt]{diff-
FD-E1-E2-kbT.eps}\includegraphics[width=65.04366pt]{diff-
FD-E3-E4-kbT.eps}\includegraphics[width=65.04366pt]{diff-FD-MA-zoom-
kbT.eps}\end{array}$

Figure S.4: Statistical weight $\Delta f(E)$ corresponding to the transitions
between (a) $E_{-,+}\rightarrow E_{+,+}$ and (b) $E_{-,-}\rightarrow E_{+,-}$
at different temperatures. The dots in (a) ((b) and (c)) represent the energy
$E_{+,+}$ ($E_{+,-}$) at the MA and the arrows mark the limit of the band
$E_{+,+}$ ($E_{+,-}$) at the indicated twist angle. In (c), $\Delta f(E)$ is
shown around the MA for the transition between $E_{-,-}\rightarrow E_{+,-}$.
Calculations are done for the SOC $\lambda_{I}=3\;\mathrm{meV}$,
$\lambda_{R}=4\;\mathrm{meV}$ Alex  and in the undoped TBG ($\mu=0$).

Figure S.5 shows the behavior of the normalized GD correction
$\delta\alpha_{G}/\alpha_{G}^{0}$ as function of the chemical potential $\mu$
at $k_{B}T=25\,\mathrm{meV}$ and for the FMR energy
$\hbar\Omega=0.06\,\mathrm{meV}$. The decrease of $\delta\alpha_{G}/$ is a
consequence of the thermal weight.

Figure S.5: Normalized GD correction $\delta\alpha_{G}/\alpha_{G}^{0}$ as
function of the chemical potential $\mu$ at $k_{B}T=25\,\mathrm{meV}$. The
upper limit of is $mu_{c}=\hbar v_{F}k_{c}$ corresponding to the momentum
cutoff $k_{c}=\frac{q_{0}}{2}$. Calculations are done for the SOC
$\lambda_{I}=3\;\mathrm{meV}$, $\lambda_{R}=4\;\mathrm{meV}$ Alex ,
$k_{B}T=25\,\mathrm{meV}$ and for the FMR energy
$\hbar\Omega=0.06\,\mathrm{meV}$.