# Superconductivity from electronic interactions and spin-orbit enhancement in
bilayer and trilayer graphene.

Alejandro Jimeno-Pozo<EMAIL_ADDRESS>$Imdea\ Nanoscience,\
Faraday\ 9,\ 28049\ Madrid,\ Spain$ Héctor Sainz-Cruz  $Imdea\ Nanoscience,\
Faraday\ 9,\ 28049\ Madrid,\ Spain$ Tommaso Cea  $Imdea\ Nanoscience,\
Faraday\ 9,\ 28049\ Madrid,\ Spain$ $Department\ of\ Physical\ and\ Chemical\
Sciences,\ University\ of\ L^{\prime}Aquila,\ via\ Vetoio,\ Coppito,\ 67100\
L^{\prime}Aquila,\ Italy$ Pierre A. Pantaleón  $Imdea\ Nanoscience,\ Faraday\
9,\ 28049\ Madrid,\ Spain$ Francisco Guinea  $Imdea\ Nanoscience,\ Faraday\
9,\ 28049\ Madrid,\ Spain$ $Donostia\ International\ Physics\ Center,\ Paseo\
Manuel\ de\ Lardizabal\ 4,\ 20018\ San\ Sebastian,\ Spain$ $Ikerbasque\
Foundation,Maria\ de\ Haro\ 3,\ 48013\ Bilbao,\ Spain$

(28th August 2024)

###### Abstract

We discuss a Kohn-Luttinger-like mechanism for superconductivity in Bernal
bilayer graphene and rhombohedral trilayer graphene. Working within the
continuum model description, we find that the screened long-range Coulomb
interaction alone gives rise to superconductivity with critical temperatures
that agree with experiments. We observe that the order parameter changes sign
between valleys, which implies that both materials are valley-singlet, spin-
triplet superconductors. Adding Ising spin-orbit coupling leads to a
significant enhancement in the critical temperature, also in line with
experiment, and the superconducting order parameter shows locking between the
spin and valley degrees of freedom.

Introduction. Recent experiments report cascades of correlated phases in
Bernal bilayer graphene (BBG) [1, 2, 3]. One of them is spin-polarized
superconductivity [1] with a critical temperature of $T_{c}\approx 26$ mK for
fillings near the hole-doped van Hove singularity (vHs) when an out-of-plane
electric field together with an in-plane magnetic field are applied on the
system. Moreover, the authors of Ref. [4] have assembled a heterostructure in
which a transition metal dichalcogenide (TMD) and BBG enter in synergy, making
superconductivity appear over a broader range of electron filling and magnetic
field, even without the latter, and enabling a striking increment in critical
temperature to $T_{c}\approx$ 260 mK, an effect that is attributed to Ising
spin-orbit coupling (SOC). The stability and structural simplicity of BBG are
clear advantages for experimental reproducibility, a major obstacle in Moiré
materials [5], in which all samples differ due to angle disorder [6, 7] and
strains [8]. A previous experiment showed that rhombohedral trilayer graphene
(RTG) is a superconductor as well [9] and several theories about it have been
proposed [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]. A key piece in the
puzzle of graphene superconductors is that they all display cascades of
flavour-symmetry-breaking phase transitions [1, 21, 22, 23, 2, 3], perhaps
hinting at a common origin of superconductivity.

In these materials, and in BBG in particular, the proximity of
superconductivity to flavour-polarized metallic phases and its Pauli limit
violation point to an unconventional spin-triplet pairing mediated by
electrons. In Ref. [24], a model is discussed in which short-range momentum-
independent interactions and proximity to symmetry broken phases induce spin-
triplet $f$-wave pairing. The authors of Ref. [25] propose a superconductivity
from repulsion mechanism, in which the Coulomb interaction dressed by soft
quantum-critical modes drives pairing. In their model, the magnetic field
induces spin imbalance, which makes the interaction acquire a dependence on
frequency or on soft-mode momenta, leading to valley-singlet, spin-triplet
s-wave or valley-triplet spin-triplet p-wave pairings, respectively. In Refs.
[26, 27], a pairing mechanism mediated by acoustic phonons is investigated,
which is compatible with spin-singlet s-wave and spin-triplet f-wave pairings.
However, the authors assume that the Coulomb interaction is detrimental for
superconductivity. In contrast, we find that the Coulomb interaction alone
enables superconductivity, see also [25].

Here, we present a framework in which superconductivity in BBG emerges only
from the long-range Coulomb interaction, arguably the simplest explanation.
Due to screening by particle-hole pairs, this repulsive interaction between
electrons becomes attractive and leads to pairing. We use a diagrammatic
technique similar to the Kohn-Luttinger approach [28], which we have already
applied to twisted bilayer and trilayer graphene [29, 30], as well as to RTG
[15], thus allowing for a direct comparison of superconductivity in these
systems and in BBG. We find superconductivity with critical temperature
comparable to the experimental one [1], near the vHs. Including Ising SOC
leads to significant increments in the critical temperature, as seen
experimentally [4]. We observe that the order parameter (OP) changes sign
within each valley. Adding a short-range Hubbard repulsion, we conclude that
the OP also changes sign between valleys, showing that BBG is a spin-triplet
superconductor. The rest of the paper is organized as follows: first we
describe the continuum model and the Kohn-Luttinger-like framework for
superconductivity. Then we present and discuss the results of the critical
temperature, the superconducting OP and the positive effect of adding SOC.

Figure 1: (a) Crystal lattice and hoppings of BBG. (b) Brillouin zone of BBG
and RTG. We perform the calculation of superconductivity considering electron
states in hexagonal grids around the Dirac points, with ultraviolet momentum
cutoff $k_{\Lambda}=0.025K_{D}$ for BBG (and $0.035K_{D}$ for RTG), with
$K_{D}=4\pi/3a$ the Dirac point modulus, which contain the Fermi surface and
where the continuum model is a good approximation. Red and blue hexagons are
plotted with $k_{\Lambda}=0.2K_{D}$ for clarity. (c, d) Continuum model
valence band and DOS of BBG and RTG with electric-field-induced gaps of $98$
and $74$ meV, respectively [31]. Dashed red lines mark the energies of the
vHs.

The continuum model. Bernal bilayer graphene (BBG) refers to two stacked
graphene layers so that atoms belonging to sublattice A of layer 1 lie over
the atoms in sublattice B of layer 2. Similarly to monolayer graphene, BBG is
a semi-metal in which the low energy bands touch at the Dirac points, but with
parabolic instead of linear dispersion [32, 33]. This band touching makes the
charge susceptibility and other susceptibilities diverge [34], which leads to
symmetry-broken phases [35, 36, 37, 38, 39, 40, 41, 42, 43], as seen in
experiments [44, 45, 46, 47, 48, 49, 50]. A perpendicular electric field opens
a gap and the bands acquire a ‘Mexican hat’ profile [51, 52], as shown in Fig.
1(c,d), making the density of states diverge logarithmically near the band
edges, causing vHs and setting the stage for new correlated phenomena [53, 54,
55, 56, 1, 2, 3, 4].

In the continuum approximation the hamiltonian of BBG can be expressed in the
basis $\\{\Psi_{A_{1}},\Psi_{B_{1}},\Psi_{A_{2}},\Psi_{B_{2}}\\}$ as [52],

$\mathcal{H}_{\xi}=\matrixquantity(V/2&v_{0}\pi^{\dagger}&-v_{4}\pi^{\dagger}&v_{3}\pi\\\
v_{0}\pi&V/2+\Delta^{\prime}&\gamma_{1}&-v_{4}\pi^{\dagger}\\\
-v_{4}\pi&\gamma_{1}&-V/2+\Delta^{\prime}&v_{0}\pi^{\dagger}\\\
v_{3}\pi^{\dagger}&-v_{4}\pi&v_{0}\pi&-V/2),$ (1)

where $\pi=\xi p_{x}+\mathrm{i}p_{y}$, $\xi=\pm 1$ is the valley index and
$\mathbf{p}=\mathbf{k}-\mathbf{K}^{\xi}$ is the momentum with respect to the
$\mathbf{K}^{\xi}$ Dirac point, $v_{i}=\frac{\sqrt{3}a\gamma_{i}}{2\hbar}$ are
effective velocities, $a$ is the lattice constant and $\gamma_{i}$ the
different hoppings in the system, shown in Fig. 1(a). $V$ corresponds to an
applied electric field and $\Delta^{\prime}$ is the energy difference between
dimer ($A_{1}$, $B_{2}$) and non-dimer ($A_{2}$, $B_{1}$) sites. In Ref. [57]
we list the values [52, 58], and describe the hamiltonian of RTG.

Kohn-Luttinger-like superconductivity. We propose that electronic interactions
alone are enough to induce superconductivity in BBG and RTG and that the
pairing glue is the screened long-range Coulomb potential $V_{scr}$. We employ
a similar approach to the Kohn-Luttinger theory [28], following an analogous
method as in [15]. We use the Random Phase Approximation to take into account
the screened direct interaction to infinite order, while neglecting the
contribution from the exchange interaction. The multiplicity of the direct
diagrams equals the number of flavours, $N_{f}=4$, so the approximation can be
considered an expansion in powers of $1/N_{f}$. Under these assumptions the
screened interaction is given by,

$V_{scr}(\mathbf{q})=\frac{V_{C}(\mathbf{q})}{1-\Pi(\mathbf{q})V_{C}(\mathbf{q})},$
(2)

where $V_{C}(\mathbf{q})=\frac{2\pi
e^{2}}{\epsilon|\mathbf{q}|}\tanh{\left(d|\mathbf{q}|\right)}$ is the bare
Coulomb potential, with $e$ the electron charge, $\epsilon=4$ the dielectric
constant associated hBN encapsulation and $d=40$ nm the distance to the
metallic gate. $\Pi(\mathbf{q})$ corresponds to the zero-frequency limit of
the charge susceptibility,

$\Pi(\mathbf{q})=\frac{N_{f}}{\Omega}\sum_{\mathbf{k},m,n}\frac{f(\epsilon_{n,\mathbf{k}})-f(\epsilon_{m,\mathbf{k}+\mathbf{q}})}{\epsilon_{n,\mathbf{k}}-\epsilon_{m,\mathbf{k}+\mathbf{q}}}\absolutevalue{\innerproduct{\Psi_{m,\mathbf{k}+\mathbf{q}}}{\Psi_{n,\mathbf{k}}}}^{2},$
(3)

where $m,n$ are band indexes, $f$ is the Fermi-Dirac distribution and
$\epsilon=E-\mu$ with $\mu$ being the Fermi energy, E are the eigenvalues,
$\Psi$ the eigenvectors of (1) and $\Omega$ is the area of the system.

Figure 2: Superconducting critical temperature due to only the screened long-
range Coulomb interaction (blue), and including also the short-range
interaction (green), and DOS (a.u., dashed black) versus Fermi energy in (a)
BBG with hole doping, near the vHs ($n_{e}\approx-0.87\cdot 10^{12}$ cm-2), we
observe a $T_{c}\approx 10$ mK, and (b) RTG with hole doping, near the vHs
($n_{e}\approx-0.6\cdot 10^{12}$ cm-2), $T_{c}\approx 33$ mK and it raises to
$T_{c}\approx 65$ mK after including the short-range interaction. The top
insets show the shape of the Fermi surfaces, which experience a change in
topology at the vHs due to Lifshitz transitions. Right insets: screened
Coulomb potential to $k_{\Lambda}$ at the energy for which $T_{c}$ is maximum.

Once we perform this calculation within each valley, we include a short-range,
repulsive Hubbard $U=3$ eV [59, 60], which allows electrons to exchange valley
and fixes the sign of the order parameter in each valley. We define the
dimensionless anomalous expectation values in both valleys
$\tilde{\Delta}^{+,i,j}(\mathbf{k})=\expectationvalue{c_{\mathbf{k},i,\mathbf{K}^{+},\uparrow}^{\dagger}c_{-\mathbf{k},j,\mathbf{K}^{-},\downarrow}^{\dagger}}$
and
$\tilde{\Delta}^{-,i,j}(\mathbf{k})=\expectationvalue{c_{\mathbf{k},i,\mathbf{K}^{-},\uparrow}^{\dagger}c_{-\mathbf{k},j,\mathbf{K}^{+},\downarrow}^{\dagger}}$.
Therefore the gap equation is given by,

$\begin{split}\tilde{\Delta}^{+,i,j}\left(\mathbf{k}\right)&=-\frac{K_{B}T}{\Omega}\sum_{\mathbf{q},\omega}\sum_{i^{\prime},j^{\prime}}\\\
&V_{scr}\left(\mathbf{k}-\mathbf{q}\right)G_{\mathbf{K}^{+}}^{i,i^{\prime}}\left(\mathbf{q},\mathrm{i}\omega\right)G_{\mathbf{K}^{-}}^{j,j^{\prime}}\left(-\mathbf{q},-\mathrm{i}\omega\right)\tilde{\Delta}^{+,i^{\prime},j^{\prime}}\left(\mathbf{q}\right)\\\
&+UG_{\mathbf{K}^{-}}^{i,i^{\prime}}\left(\mathbf{q},\mathrm{i}\omega\right)G_{\mathbf{K}^{+}}^{j,j^{\prime}}\left(-\mathbf{q},-\mathrm{i}\omega\right)\tilde{\Delta}^{-,j^{\prime},i^{\prime}}\left(-\mathbf{q}\right),\end{split}$
(4)

where indexes $i,j$ denote layer and sublattice degrees of freedom, $\omega$
is the Matsubara frequency and
$G_{K^{\xi}}^{i,j}\left(\mathbf{q},\mathrm{i}\omega\right)$ are electronic
Green’s functions. Eq. (4) can be understood in terms of the convolution of
the gap with a kernel,

$\tilde{\Delta}(\mathbf{k})=\matrixquantity(\tilde{\Delta}^{+}(\mathbf{k})\\\
\tilde{\Delta}^{-}(\mathbf{k}))=\sum_{\mathbf{k}^{\prime}}\matrixquantity(\Gamma^{+}(\mathbf{k},\mathbf{k}^{\prime})&\tilde{U}(\mathbf{k},\mathbf{k}^{\prime})\\\
\tilde{U}^{\dagger}(\mathbf{k},\mathbf{k}^{\prime})&\Gamma^{-}(\mathbf{k},\mathbf{k}^{\prime}))\tilde{\Delta}(\mathbf{k}^{\prime}),$
(5)

so that the solution of this gap equation is achieved when the maximum
eigenvalue of the hermitian kernel reaches a value of 1, indicating the
$T_{c}$ at which the scattering amplitude for pairs of carriers of opposite
momenta and energy within the Fermi surface has a pole. The corresponding
eigenvector is the superconducting order parameter. The long-range Coulomb
kernel is,

$\Gamma_{m,n}^{\xi}(\mathbf{k},\mathbf{k^{\prime}})=-\frac{V_{scr}(\mathbf{k}-\mathbf{k^{\prime}})}{\Omega}S^{\xi,\xi}_{m,n}(\mathbf{k},\mathbf{k^{\prime}})R^{\xi,\xi}_{m,m}(\mathbf{k},\mathbf{k})R^{\xi,\xi}_{n,n}(\mathbf{k}^{\prime},\mathbf{k}^{\prime}),$
(6)

while the short-range kernel is,

$\tilde{U}_{m,n}(\mathbf{k},\mathbf{k^{\prime}})=-\frac{U}{\Omega}S^{+,-}_{m,n}(\mathbf{k},\mathbf{k^{\prime}})R^{+,-}_{m,m}(\mathbf{k},\mathbf{k})R^{-,+}_{n,n}(\mathbf{k}^{\prime},\mathbf{k}^{\prime}),$
(7)

with
$R^{\xi,\xi^{\prime}}_{m,n}(\mathbf{k},\mathbf{k^{\prime}})=\sqrt{\left(f(-\epsilon^{\xi}_{m,\mathbf{k}})-f(\epsilon^{\xi^{\prime}}_{n,\mathbf{k}^{\prime}})\right)/\left(\epsilon^{\xi}_{m,\mathbf{k}}+\epsilon^{\xi^{\prime}}_{n,\mathbf{k}^{\prime}}\right)}$
and
$S^{\xi,\xi^{\prime}}_{m,n}(\mathbf{k},\mathbf{k^{\prime}})=\absolutevalue{\innerproduct{\Psi^{\xi}_{m,\mathbf{k}}}{\Psi^{\xi^{\prime}}_{n,\mathbf{k^{\prime}}}}}^{2}$.
Once Eq. (5) is solved, gaps with dimensions of energy can be obtained as,
$\Delta^{+}(\mathbf{k})=V_{scr}(\mathbf{k})\tilde{\Delta}^{+}(\mathbf{k})+U\tilde{\Delta}^{-}(\mathbf{k})$
and
$\Delta^{-}(\mathbf{k})=V_{scr}(\mathbf{k})\tilde{\Delta}^{-}(\mathbf{k})+U\tilde{\Delta}^{+}(\mathbf{k})$.

Results. In Figure 2 we present the superconducting critical temperature of
hole-doped BBG and RTG as a function of Fermi energy. We obtain a $T_{c}$ of
10 mK for BBG and 65 mK for RTG, which match well with those observed
experimentally [1, 9]. In both materials superconductivity survives only in
narrow energy intervals around the vHs. It is worth noting that in RTG,
including the short-range interaction doubles $T_{c}$, but it does not
increase $T_{c}$ in BBG. To obtain $T_{c}$, we use a $\mathbf{k}$-point
density of up to $\mathcal{O}(10^{7})\r{A}^{2}$ for BBG and
$\mathcal{O}(10^{6})\r{A}^{2}$ for RTG. Since states close to the Fermi
surface give the principal contribution to Eq. (5), we cut off phase space by
considering only states with $\absolutevalue{\epsilon_{n,\mathbf{k}}}\leqslant
30$ meV. These continuum model results are in good agreement with our previous
tight-binding calculation for RTG [15]. The insets in Fig. 2 show the screened
potential in momentum space, which is always repulsive, i.e. positive, and its
value increases with increasing momentum. This dependence favors non s-wave
superconductivity [61].

Moreover, inspired by the authors of Ref. [4], we study the effect of stacking
a TMD on top of BBG. We assume that the main effect is an induced Ising spin-
orbit coupling (SOC) in BBG, which breaks the equivalence between Cooper pairs
$\ket{\mathbf{K}^{+},\uparrow;\mathbf{K}^{-},\downarrow}$ and
$\ket{\mathbf{K}^{+},\downarrow;\mathbf{K}^{-},\uparrow}$, so the order
parameter cannot be defined as a spin singlet or a spin triplet. The Ising
hamiltonian which we add to (1) is just
$\mathcal{H}_{I,\xi,s}=s\xi\lambda_{I}\mathbb{I}$, where s is the spin index.
This interaction does not reshape the band structure, it just promotes one
type of Cooper pair, by lowering the energy of two spin-valley flavours and
raising that of the other two, splitting the original vHs into two.

Figure 3: Ising spin-orbit coupling enhancement of critical temperature in
hole-doped (a) BBG, for which $T_{c}$ augments by a factor of 4 and saturates
at $\lambda_{I}=4$ meV. (b) RTG, for which $T_{c}$ due to long-range
interactions (blue) augments by a factor of 2 and saturates at $\lambda_{I}=1$
meV. In contrast, Ising SOC is detrimental to short-range interactions, so
$T_{c}$ does not increase much when they are included (green).

As shown in Figure 3, Ising SOC has a very positive effect on
superconductivity in both materials: it increases $T_{c}$ by a factor of 4 in
BBG, to 40 mK, congruent with the experiment, in which a factor of 10
increment was seen [4] and in RTG it brings forth a comparable increment of
$T_{c}$ due to long-range interactions only. The two vHs lead to two
superconducting domes as a function of Fermi energy, as observed in Ref. [4],
see [57]. It is worth noting that the spin-orbit enhancements of $T_{c}$
saturate once the SOC fully polarizes the non-superconducting state into a
half-metal [57].

Our model also yields the superconducting order parameters (OPs), shown in
Figure 4. The OPs have structure within each valley: in BBG they display
hotspots and sign changes along the edges of the Fermi surface. The sign is
opposite for the inner and outer Fermi edges. In RTG, the maximum eigenvalue
of the kernel is degenerate, leading to $C_{3}$-symmetry breaking. Such
degeneracy is not typical of conventional superconductors and hints at an
exotic type of superconductivity. We observe that the OPs in both materials
change sign between valleys, forming valley-singlets. Since the overall
electron wavefunction must be antisymmetric, this implies that the pairs are
spin-triplets. The order parameter that we find in the presence of spin-orbit
coupling is similar to the spin-valley locking model in [62]. Our results are
consistent with the phenomenological model discussed in [63].

Figure 4: (a) Superconducting order parameter (OP) of BBG in valley $K^{+}$,
(b) in valley $K^{-}$, with hole doping and $T_{c}\approx 10$ mK. Within a
single valley, the OP has hotspots along the edges of the Fermi surface and
changes sign between the inner and outer edges. (c) OP of RTG in valley
$K^{+}$, (d) in valley $K^{-}$, with hole doping and $T_{c}\approx 65$ mK,
showing intensity stripes along the edges of the Fermi surface. In both
materials, the OP changes sign between valleys, which means that BBG and RTG
are valley-singlet, spin-triplet superconductors.

Discussion. BBG and RTG are 2D superconductors with great structural stability
and low disorder. These are decisive advantages to understand the physics at
play, as they enable very reproducible experiments. There is also compelling
evidence that both support spin-triplet superconducting phases and hence the
pairing is most likely unconventional and mediated by electrons. Here, we show
that a Kohn-Luttinger-like mechanism for superconductivity suffices to give
rise to superconductivity in BBG and RTG, with critical temperatures in good
agreement with experiments on both materials [1, 9, 4]. The screened long-
range Coulomb interaction is the sole responsible for pairing. At fillings
near the vHs, electron interactions strongly screen the Coulomb potential,
which glues pairs of carriers, one in each valley, and causes
superconductivity. Of course, other excitations may dress the Coulomb
interaction and contribute to superconductivity in BBG [24, 25, 26] and RTG
[10, 11, 12, 13, 14, 16, 17, 18, 19, 20].

Furthermore, our results support the proposal in Ref. [4] that Ising SOC is
what causes the enhancement of $T_{c}$ in BBG when the TMD WSe2 is stacked on
top of it. Based on the results presented here, we predict that the same idea
will boost superconductivity in RTG by a comparable amount as in BBG. Since
the four times increment we obtain in BBG underestimates the ten times
increment observed in the experiment, it is likely that our result for RTG
similarly underestimates the SOC enhancement that superconductivity will
experience in this material. Hence, we conjecture that RTG with WSe2 stacked
on top will reach $T_{c}\sim 1$ K, putting it on par with magic-angle twisted
bilayer graphene and the other twisted stacks [64, 65, 66, 67, 68, 69, 70].

The OPs reveal two key aspects of superconductivity. First, the OPs show sign
changes within each valley, which are characteristic of weak coupling
superconductivity, in which electrons interact and form pairs with high
angular momentum, such as those with p-wave or f-wave symmetry. Another
consequence of this intravalley structure is that long-range disorder will
suppress superconductivity in these materials, unlike in twisted stacks.
Second, the sign of the OPs changes between valleys, hence the pairs are spin-
triplets.

In BBG and RTG, the screened long-range Coulomb interaction suffices to reach
critical temperatures that agree with experiments, while in twisted stacks it
is necessary to include Umklapp processes involving large momentum transfers
and dressing by phonons, pointing to a rich interplay between interactions and
wavefunction complexity [29, 30]. Also, the rigidity of the bands of BBG and
RTG leads to narrow superconducing sleeves, because the phase vanishes when
the Fermi energy passes the vHs. In contrast, in twisted stacks the bands
reshape with filling and there is a pinning of the Fermi energy to the vHs
[71], allowing for wide superconducting domes.

In conclusion, we find that superconductivity emerges in BBG and RTG from the
screened long-range Coulomb interaction alone. This proposal, based on Kohn-
Luttinger theory, yields results consistent with three major experimental
advances [1, 9, 4]. We observe critical temperatures that agree with
experiments. The inclusion of Ising SOC results in prominent increments in
critical temperature for both materials, in harmony with the experimental
observation in BBG. In the absence of Ising SOC the OPs show that both
materials are valley-singlet, spin-triplet superconductors, while this
classification cannot be applied when SOC is present since it locks the spin
and valley degrees of freedom.

Acknowledgments. We are thankful to A. V. Chubukov for illuminating
discussions. We acknowledge support from the Severo Ochoa programme for
centres of excellence in R&D (Grant No. SEV-2016-0686, Ministerio de Ciencia e
Innovación, Spain); from the European Commission, within the Graphene
Flagship, Core 3, grant number 881603 and from grants NMAT2D (Comunidad de
Madrid, Spain) and SprQuMat (Ministerio de Ciencia e Innovación, Spain).

## References

  * Zhou  _et al._ [2022] H. Zhou, L. Holleis, Y. Saito, L. Cohen, W. Huynh, C. L. Patterson, F. Yang, T. Taniguchi, K. Watanabe, and A. F. Young, Isospin magnetism and spin-polarized superconductivity in bernal bilayer graphene, Science 375, 774 (2022).
  * de la Barrera  _et al._ [2022] S. C. de la Barrera, S. Aronson, Z. Zheng, K. Watanabe, T. Taniguchi, Q. Ma, P. Jarillo-Herrero, and R. Ashoori, Cascade of isospin phase transitions in bernal-stacked bilayer graphene at zero magnetic field, Nature Physics 18, 771 (2022).
  * Seiler  _et al._ [2022] A. M. Seiler, F. R. Geisenhof, F. Winterer, K. Watanabe, T. Taniguchi, T. Xu, F. Zhang, and R. T. Weitz, Quantum cascade of correlated phases in trigonally warped bilayer graphene, Nature 608, 298 (2022).
  * Zhang  _et al._ [2022a] Y. Zhang, R. Polski, A. Thomson, É. Lantagne-Hurtubise, C. Lewandowski, H. Zhou, K. Watanabe, T. Taniguchi, J. Ali, and S. Nadj-Perge, Spin-orbit enhanced superconductivity in bernal bilayer graphene, arXiv preprint arXiv:2205.05087 10.48550/arXiv.2205.05087 (2022a).
  * Lau  _et al._ [2022] C. N. Lau, M. W. Bockrath, K. F. Mak, and F. Zhang, Reproducibility in the fabrication and physics of moiré materials, Nature 602, 41 (2022).
  * Uri  _et al._ [2020] A. Uri, S. Grover, Y. Cao, J. A. Crosse, K. Bagani, D. Rodan-Legrain, Y. Myasoedov, K. Watanabe, T. Taniguchi, P. Moon, M. Koshino, P. Jarillo-Herrero, and E. Zeldov, Mapping the twist-angle disorder and landau levels in magic-angle graphene, Nature 581, 47 (2020).
  * Schäpers  _et al._ [2022] A. Schäpers, J. Sonntag, L. Valerius, B. Pestka, J. Strasdas, K. Watanabe, T. Taniguchi, L. Wirtz, M. Morgenstern, B. Beschoten, R. J. Dolleman, and C. Stampfer, Raman imaging of twist angle variations in twisted bilayer graphene at intermediate angles, 2D Materials 9, 045009 (2022).
  * Kazmierczak  _et al._ [2021] N. P. Kazmierczak, M. V. Winkle, C. Ophus, K. C. Bustillo, S. Carr, H. G. Brown, J. Ciston, T. Taniguchi, K. Watanabe, and D. K. Bediako, Strain fields in twisted bilayer graphene, Nature Materials 20, 956 (2021).
  * Zhou  _et al._ [2021a] H. Zhou, T. Xie, T. Taniguchi, K. Watanabe, and A. F. Young, Superconductivity in rhombohedral trilayer graphene, Nature 598, 434 (2021a).
  * Chou  _et al._ [2021] Y.-Z. Chou, F. Wu, J. D. Sau, and S. D. Sarma, Acoustic-phonon-mediated superconductivity in rhombohedral trilayer graphene, Physical Review Letters 127, 10.1103/physrevlett.127.187001 (2021).
  * Dai  _et al._ [2021] H. Dai, J. Hou, X. Zhang, Y. Liang, and T. Ma, Mott insulating state and d+ i d superconductivity in an abc graphene trilayer, Physical Review B 104, 10.1103/physrevb.104.035104 (2021).
  * Dong and Levitov [2021] Z. Dong and L. Levitov, Superconductivity in the vicinity of an isospin-polarized state in a cubic dirac band, arXiv preprint arXiv:2109.01133 10.48550/arxiv.2109.01133 (2021).
  * Ghazaryan  _et al._ [2021] A. Ghazaryan, T. Holder, M. Serbyn, and E. Berg, Unconventional superconductivity in systems with annular fermi surfaces: Application to rhombohedral trilayer graphene, Physical Review Letters 127, 10.1103/physrevlett.127.247001 (2021).
  * Chatterjee  _et al._ [2022] S. Chatterjee, T. Wang, E. Berg, and M. P. Zaletel, Inter-valley coherent order and isospin fluctuation mediated superconductivity in rhombohedral trilayer graphene, Nature Communications 13, 10.1038/s41467-022-33561-w (2022).
  * Cea  _et al._ [2022] T. Cea, P. A. Pantaleón, V. T. Phong, and F. Guinea, Superconductivity from repulsive interactions in rhombohedral trilayer graphene: A kohn-luttinger-like mechanism, Physical Review B 105, 10.1103/physrevb.105.075432 (2022).
  * You and Vishwanath [2022] Y.-Z. You and A. Vishwanath, Kohn-luttinger superconductivity and intervalley coherence in rhombohedral trilayer graphene, Physical Review B 105, 10.1103/physrevb.105.134524 (2022).
  * Szabó and Roy [2022] A. L. Szabó and B. Roy, Metals, fractional metals, and superconductivity in rhombohedral trilayer graphene, Physical Review B 105, 10.1103/physrevb.105.l081407 (2022).
  * Qin  _et al._ [2022] W. Qin, C. Huang, T. Wolf, N. Wei, I. Blinov, and A. H. MacDonald, Functional renormalization group study of superconductivity in rhombohedral trilayer graphene, arXiv preprint arXiv:2203.09083 10.48550/arxiv.2203.09083 (2022).
  * Dai  _et al._ [2022] H. Dai, R. Ma, X. Zhang, and T. Ma, Quantum monte carlo study of superconductivity in rhombohedral trilayer graphene under an electric field, arXiv preprint arXiv:2204.06222 10.48550/arxiv.2204.06222 (2022).
  * Lu  _et al._ [2022] D.-C. Lu, T. Wang, S. Chatterjee, and Y.-Z. You, Correlated metals and unconventional superconductivity in rhombohedral trilayer graphene: A renormalization group analysis, Physical Review B 106, 10.1103/physrevb.106.155115 (2022).
  * Zondiner  _et al._ [2020] U. Zondiner, A. Rozen, D. Rodan-Legrain, Y. Cao, R. Queiroz, T. Taniguchi, K. Watanabe, Y. Oreg, F. von Oppen, A. Stern, E. Berg, P. Jarillo-Herrero, and S. Ilani, Cascade of phase transitions and dirac revivals in magic-angle graphene, Nature 582, 203 (2020).
  * Wong  _et al._ [2020] D. Wong, K. P. Nuckolls, M. Oh, B. Lian, Y. Xie, S. Jeon, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Cascade of electronic transitions in magic-angle twisted bilayer graphene, Nature 582, 198 (2020).
  * Zhou  _et al._ [2021b] H. Zhou, T. Xie, A. Ghazaryan, T. Holder, J. R. Ehrets, E. M. Spanton, T. Taniguchi, K. Watanabe, E. Berg, M. Serbyn, and A. F. Young, Half- and quarter-metals in rhombohedral trilayer graphene, Nature 598, 429 (2021b).
  * Szabó and Roy [2022] A. L. Szabó and B. Roy, Competing orders and cascade of degeneracy lifting in doped bernal bilayer graphene, Phys. Rev. B 105, L201107 (2022).
  * Dong  _et al._ [2022] Z. Dong, A. V. Chubukov, and L. Levitov, Spin-triplet superconductivity at the onset of isospin order in biased bilayer graphene, arXiv preprint arXiv:2205.13353 10.48550/arxiv.2205.13353 (2022).
  * Chou  _et al._ [2022a] Y.-Z. Chou, F. Wu, J. D. Sau, and S. D. Sarma, Acoustic-phonon-mediated superconductivity in bernal bilayer graphene, Physical Review B 105, 10.1103/physrevb.105.l100503 (2022a).
  * Chou  _et al._ [2022b] Y.-Z. Chou, F. Wu, and S. D. Sarma, Enhanced superconductivity through virtual tunneling in bernal bilayer graphene coupled to wse $\\_2$, Physical Review B 106, 10.1103/physrevb.106.l180502 (2022b).
  * Kohn and Luttinger [1965] W. Kohn and J. M. Luttinger, New mechanism for superconductivity, Physical Review Letters 15, 524 (1965).
  * Cea and Guinea [2021] T. Cea and F. Guinea, Coulomb interaction, phonons, and superconductivity in twisted bilayer graphene, Proceedings of the National Academy of Sciences 118, 10.1073/pnas.2107874118 (2021).
  * Phong  _et al._ [2021] V. T. Phong, P. A. Pantaleón, T. Cea, and F. Guinea, Band structure and superconductivity in twisted trilayer graphene, Physical Review B 104, 10.1103/physrevb.104.l121116 (2021).
  * [31] To open these gaps one need to include electric fields of around $V=110$ meV and $V=42.5$ meV, respectively. It is worth noting that there is considerable uncertainty in the relation between the applied fields and the gaps, see Refs. [51, 72].
  * McCann and Fal’ko [2006] E. McCann and V. I. Fal’ko, Landau-level degeneracy and quantum hall effect in a graphite bilayer, Physical Review Letters 96, 10.1103/physrevlett.96.086805 (2006).
  * Novoselov  _et al._ [2006] K. S. Novoselov, E. McCann, S. V. Morozov, V. I. Fal’ko, M. I. Katsnelson, U. Zeitler, D. Jiang, F. Schedin, and A. K. Geim, Unconventional quantum hall effect and berry’s phase of 2$\pi$ in bilayer graphene, Nature Physics 2, 177 (2006).
  * Nilsson  _et al._ [2006] J. Nilsson, A. H. C. Neto, N. M. R. Peres, and F. Guinea, Electron-electron interactions and the phase diagram of a graphene bilayer, Physical Review B 73, 10.1103/physrevb.73.214418 (2006).
  * Zhang  _et al._ [2011] F. Zhang, J. Jung, G. A. Fiete, Q. Niu, and A. H. MacDonald, Spontaneous quantum hall states in chirally stacked few-layer graphene systems, Phys. Rev. Lett. 106, 156801 (2011).
  * Min  _et al._ [2008] H. Min, G. Borghi, M. Polini, and A. H. MacDonald, Pseudospin magnetism in graphene, Physical Review B 77, 10.1103/physrevb.77.041407 (2008).
  * Vafek [2010] O. Vafek, Interacting fermions on the honeycomb bilayer: From weak to strong coupling, Physical Review B 82, 10.1103/physrevb.82.205106 (2010).
  * Vafek and Yang [2010] O. Vafek and K. Yang, Many-body instability of coulomb interacting bilayer graphene: Renormalization group approach, Physical Review B 81, 10.1103/physrevb.81.041401 (2010).
  * Nandkishore and Levitov [2010] R. Nandkishore and L. Levitov, Dynamical screening and excitonic instability in bilayer graphene, Physical Review Letters 104, 10.1103/physrevlett.104.156803 (2010).
  * Lemonik  _et al._ [2012] Y. Lemonik, I. Aleiner, and V. I. Fal'ko, Competing nematic, antiferromagnetic, and spin-flux orders in the ground state of bilayer graphene, Physical Review B 85, 10.1103/physrevb.85.245451 (2012).
  * Zhang and MacDonald [2012] F. Zhang and A. H. MacDonald, Distinguishing spontaneous quantum hall states in bilayer graphene, Physical Review Letters 108, 10.1103/physrevlett.108.186804 (2012).
  * Cvetkovic  _et al._ [2012] V. Cvetkovic, R. E. Throckmorton, and O. Vafek, Electronic multicriticality in bilayer graphene, Physical Review B 86, 10.1103/physrevb.86.075467 (2012).
  * Martin  _et al._ [2010] J. Martin, B. E. Feldman, R. T. Weitz, M. T. Allen, and A. Yacoby, Local compressibility measurements of correlated states in suspended bilayer graphene, Physical Review Letters 105, 10.1103/physrevlett.105.256806 (2010).
  * Weitz  _et al._ [2010] R. T. Weitz, M. T. Allen, B. E. Feldman, J. Martin, and A. Yacoby, Broken-symmetry states in doubly gated suspended bilayer graphene, Science 330, 812 (2010).
  * Mayorov  _et al._ [2011] A. S. Mayorov, D. C. Elias, M. Mucha-Kruczynski, R. V. Gorbachev, T. Tudorovskiy, A. Zhukov, S. V. Morozov, M. I. Katsnelson, A. K. Geim, and K. S. N. and, Interaction-driven spectrum reconstruction in bilayer graphene, Science 333, 860 (2011).
  * Velasco  _et al._ [2012] J. Velasco, L. Jing, W. Bao, Y. Lee, P. Kratz, V. Aji, M. Bockrath, C. N. Lau, C. Varma, R. Stillwell, D. Smirnov, F. Zhang, J. Jung, and A. H. MacDonald, Transport spectroscopy of symmetry-broken insulating states in bilayer graphene, Nature Nanotechnology 7, 156 (2012).
  * Freitag  _et al._ [2012] F. Freitag, J. Trbovic, M. Weiss, and C. Schönenberger, Spontaneously gapped ground state in suspended bilayer graphene, Physical Review Letters 108, 10.1103/physrevlett.108.076602 (2012).
  * Ju  _et al._ [2015] L. Ju, Z. Shi, N. Nair, Y. Lv, C. Jin, J. Velasco, C. Ojeda-Aristizabal, H. A. Bechtel, M. C. Martin, A. Zettl, J. Analytis, and F. Wang, Topological valley transport at bilayer graphene domain walls, Nature 520, 650 (2015).
  * Li  _et al._ [2017] J. I. A. Li, C. Tan, S. Chen, Y. Zeng, T. Taniguchi, K. Watanabe, J. Hone, and C. R. Dean, Even-denominator fractional quantum hall states in bilayer graphene, Science 358, 648 (2017).
  * Zibrov  _et al._ [2017] A. A. Zibrov, C. Kometter, H. Zhou, E. M. Spanton, T. Taniguchi, K. Watanabe, M. P. Zaletel, and A. F. Young, Tunable interacting composite fermion phases in a half-filled bilayer-graphene landau level, Nature 549, 360 (2017).
  * Zhang  _et al._ [2009] Y. Zhang, T.-T. Tang, C. Girit, Z. Hao, M. C. Martin, A. Zettl, M. F. Crommie, Y. R. Shen, and F. Wang, Direct observation of a widely tunable bandgap in bilayer graphene, Nature 459, 820 (2009).
  * McCann and Koshino [2013] E. McCann and M. Koshino, The electronic properties of bilayer graphene, Reports on Progress in Physics 76, 056503 (2013).
  * Maher  _et al._ [2014] P. Maher, L. Wang, Y. Gao, C. Forsythe, T. Taniguchi, K. Watanabe, D. Abanin, Z. Papić, P. Cadden-Zimansky, J. Hone, P. Kim, and C. R. Dean, Tunable fractional quantum hall phases in bilayer graphene, Science 345, 61 (2014).
  * Lee  _et al._ [2014] K. Lee, B. Fallahazad, J. Xue, D. C. Dillen, K. Kim, T. Taniguchi, K. Watanabe, and E. Tutuc, Chemical potential and quantum hall ferromagnetism in bilayer graphene, Science 345, 58 (2014).
  * Kou  _et al._ [2014] A. Kou, B. E. Feldman, A. J. Levin, B. I. Halperin, K. Watanabe, T. Taniguchi, and A. Yacoby, Electron-hole asymmetric integer and fractional quantum hall effect in bilayer graphene, Science 345, 55 (2014).
  * Geisenhof  _et al._ [2021] F. R. Geisenhof, F. Winterer, A. M. Seiler, J. Lenz, T. Xu, F. Zhang, and R. T. Weitz, Quantum anomalous hall octet driven by orbital magnetism in bilayer graphene, Nature 598, 53 (2021).
  * [57] See supplementary information, which includes Refs. [73, 74].
  * [58] See Ref. [75] for a revised version of the parameters of BBG.
  * Wehling  _et al._ [2011] T. O. Wehling, E. Şaşıoğlu, C. Friedrich, A. I. Lichtenstein, M. I. Katsnelson, and S. Blügel, Strength of effective coulomb interactions in graphene and graphite, Phys. Rev. Lett. 106, 236805 (2011).
  * Schüler  _et al._ [2013] M. Schüler, M. Rösner, T. O. Wehling, A. I. Lichtenstein, and M. I. Katsnelson, Optimal hubbard models for materials with nonlocal coulomb interactions: Graphene, silicene, and benzene, Phys. Rev. Lett. 111, 036601 (2013).
  * Maiti and Chubukov [2013] S. Maiti and A. V. Chubukov, Superconductivity from repulsive interaction, in  _AIP Conference Proceedings_, Vol. 1550 (2013) p. 3.
  * Lake  _et al._ [2022] E. Lake, A. S. Patri, and T. Senthil, Pairing symmetry of twisted bilayer graphene: A phenomenological synthesis, Phys. Rev. B 106, 104506 (2022).
  * Curtis  _et al._ [2022] J. B. Curtis, N. R. Poniatowski, Y. Xie, A. Yacoby, E. Demler, and P. Narang, Stabilizing fluctuating spin-triplet superconductivity in graphene via induced spin-orbit coupling, arXiv preprint arXiv:2209.10560 10.48550/arxiv.2209.10560 (2022).
  * Cao  _et al._ [2018] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices, Nature 556, 43 (2018).
  * Yankowitz  _et al._ [2019] M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Tuning superconductivity in twisted bilayer graphene, Science 363, 1059 (2019).
  * Lu  _et al._ [2019] 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, Superconductors, orbital magnets and correlated states in magic-angle bilayer graphene, Nature 574, 653 (2019).
  * Park  _et al._ [2021] J. M. Park, Y. Cao, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Tunable strongly coupled superconductivity in magic-angle twisted trilayer graphene, Nature 590, 249 (2021).
  * Hao  _et al._ [2021] Z. Hao, A. M. Zimmerman, P. Ledwith, E. Khalaf, D. H. Najafabadi, K. Watanabe, T. Taniguchi, A. Vishwanath, and P. Kim, Electric field–tunable superconductivity in alternating-twist magic-angle trilayer graphene, Science 371, 1133 (2021).
  * Zhang  _et al._ [2022b] Y. Zhang, R. Polski, C. Lewandowski, A. Thomson, Y. Peng, Y. Choi, H. Kim, K. Watanabe, T. Taniguchi, J. Alicea, F. von Oppen, G. Refael, and S. Nadj-Perge, Promotion of superconductivity in magic-angle graphene multilayers, Science 377, 1538 (2022b).
  * Park  _et al._ [2022] J. M. Park, Y. Cao, L.-Q. Xia, S. Sun, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Robust superconductivity in magic-angle multilayer graphene family, Nature Materials 21, 877 (2022).
  * Cea  _et al._ [2019] T. Cea, N. R. Walet, and F. Guinea, Electronic band structure and pinning of fermi energy to van hove singularities in twisted bilayer graphene: A self-consistent approach, Phys. Rev. B 100, 205113 (2019).
  * Lui  _et al._ [2011] C. H. Lui, Z. Li, K. F. Mak, E. Cappelluti, and T. F. Heinz, Observation of an electrically tunable band gap in trilayer graphene, Nature Physics 7, 944 (2011).
  * Zhang  _et al._ [2010] F. Zhang, B. Sahu, H. Min, and A. H. MacDonald, Band structure of $abc$-stacked graphene trilayers, Phys. Rev. B 82, 035409 (2010).
  * Zibrov  _et al._ [2018] A. A. Zibrov, P. Rao, C. Kometter, E. M. Spanton, J. I. A. Li, C. R. Dean, T. Taniguchi, K. Watanabe, M. Serbyn, and A. F. Young, Emergent dirac gullies and gully-symmetry-breaking quantum hall states in $aba$ trilayer graphene, Phys. Rev. Lett. 121, 167601 (2018).
  * Jung and MacDonald [2014] J. Jung and A. H. MacDonald, Accurate tight-binding models for the $\pi$ bands of bilayer graphene, Phys. Rev. B 89, 035405 (2014).

Supplementary information for ‘Superconductivity from electronic interactions
and spin-orbit enhancement in bilayer and trilayer graphene’

##  I The continuum model of BBG and RTG

Bernal bilayer graphene (BBG) refers to two vertically stacked graphene layers
so that atoms belonging to the sublattice A of layer 1 lie over the atoms of
the sublattice B of layer 2, while rhombohedral trilayer graphene (RTG)
alludes to a stack of three layers in which the atoms of sublattice A of layer
3 are located over the atoms of sublattice B of layer 2 while the atoms of
sublattice A of this layer lie over the atoms of sublattice B of layer 1.

Similar to monolayer graphene, BBG and RTG are semi-metals with band touching
at the Dirac points in the low energy region, but with parabolic instead of
linear dispersion [32, 33]. An external perpendicular electric field opens a
gap in both systems and the low energy bands nearly flatten around the Dirac
points while acquiring ‘Mexican hat’ profiles [52], as can be seen in Fig.
1(c,d) in the main text. The flattening of the low energy bands leads to
logarithmic divergences in the densities of states which correspond to van
Hove singularities (vHs).

In real space the BBG and RTG unit cell is defined by vectors
$\mathbf{a}_{1}=\frac{a}{2}\matrixquantity(1,\sqrt{3})$ and
$\mathbf{a}_{2}=\frac{a}{2}\matrixquantity(1,-\sqrt{3})$, where $a=2.46\r{A}$
is the graphene lattice constant. The Brillouin zone (BZ) of both systems is
defined by the reciprocal lattice vectors
$\mathbf{b}_{1}=\frac{2\pi}{a}\matrixquantity(1,\frac{1}{\sqrt{3}})$ and
$\mathbf{b}_{2}=\frac{2\pi}{a}\matrixquantity(1,-\frac{1}{\sqrt{3}})$. In the
continuum approximation the hamiltonian of BBG can be expressed in the basis
$\\{\Psi_{A_{1}},\Psi_{B_{1}},\Psi_{A_{2}},\Psi_{B_{2}}\\}$ as follows from
[52],

$\mathcal{H}_{BBG}=\matrixquantity(V/2&v_{0}\pi^{\dagger}&-v_{4}\pi^{\dagger}&v_{3}\pi\\\
v_{0}\pi&V/2+\Delta^{\prime}&\gamma_{1}&-v_{4}\pi^{\dagger}\\\
-v_{4}\pi&\gamma_{1}&-V/2+\Delta^{\prime}&v_{0}\pi^{\dagger}\\\
v_{3}\pi^{\dagger}&-v_{4}\pi&v_{0}\pi&-V/2),$ (S1)

In the case of RTG, the hamiltonian in the basis
$\\{\Psi_{A_{1}},\Psi_{B_{1}},\Psi_{A_{2}},\Psi_{B_{2}},\Psi_{A_{3}},\Psi_{B_{3}}\\}$
takes the following form, see [73],

$\mathcal{H}_{RTG}=\matrixquantity(V+\Delta_{2}+\delta&v_{0}\pi^{\dagger}&v_{4}\pi^{\dagger}&v_{3}\pi&0&\gamma_{2}/2\\\
v_{0}\pi&V+\Delta_{2}&\gamma_{1}&v_{4}\pi^{\dagger}&0&0\\\
v_{4}\pi&\gamma_{1}&-2\Delta_{2}&v_{0}\pi^{\dagger}&v_{4}\pi^{\dagger}&v_{3}\pi\\\
v_{3}\pi^{\dagger}&v_{4}\pi&v_{0}\pi&-2\Delta_{2}&\gamma_{1}&v_{4}\pi^{\dagger}\\\
0&0&v_{4}\pi&\gamma_{1}&-V+\Delta_{2}&v_{0}\pi^{\dagger}\\\
\gamma_{2}/2&0&v_{3}\pi^{\dagger}&v_{4}\pi&v_{0}\pi&-V+\Delta_{2}+\delta),$
(S2)

with $\gamma_{0}=3.1$, eV, $\gamma_{1}=0.38$ eV, $\gamma_{2}=-0.015$ eV,
$\gamma_{3}=-0.29$ eV, $\gamma_{4}=-0.141$ eV, $\delta=-0.015$ eV and
$\Delta_{2}=-0.0023$ eV, extracted from [74, 23].

Figure S1: Comparison between tight-binding (dotted red) and continuum model
(solid blue) low energy band structure of BBG near the Dirac point for (a)
$k_{\Lambda}=0.2$, (b) $k_{\Lambda}=0.1$, (c) $k_{\Lambda}=0.05$ and (d)
$k_{\Lambda}=0.025$, in units of $K_{D}$.

In order to set an ultraviolet momentum cutoff $k_{\Lambda}$, we compare the
low energy band structure of BBG within the tight-binding and continuum
models, which can be seen in Fig. S1, Ref. [15] shows a similar comparison for
RTG. We observe that for $k_{\Lambda}\leq 0.05K_{D}$, where $K_{D}=4\pi/3a$ is
the modulus of the Dirac point momentum, both models are in good agreement.
The idea is to choose the lowest possible $k_{\Lambda}$, such that there is a
high density of $\mathbf{k}$-points near the Fermi surface, and the energy
resolution of the calculation allows for the determination of the critical
temperature even if it is low. To find a convergent $T_{c}$, an energy
resolution of $K_{B}T_{c}$ is required, where $K_{B}$ is the Boltzmann
constant. However, there is a lower bound on $k_{\Lambda}$. This is because
superconductivity comes from kernel matrix elements that involve scatterings
with momentum exchange $\mathbf{q}=\mathbf{k-k^{\prime}}$, in which
$\mathbf{k}$ and $\mathbf{k^{\prime}}$ are within the Fermi surface.
Therefore, it is necessary to set $k_{\Lambda}>max{\mid\mathbf{q}\mid}$, in
order to calculate the screened potential for all relevant momentum values.
For the displacement fields studied in the main text, $k_{\Lambda}=0.025K_{D}$
suffices for the bilayer, while the trilayer requires
$k_{\Lambda}=0.035K_{D}$. It is worth noting that the size of the Fermi
surface near the vHs, and hence the adequate $k_{\Lambda}$, increase with
displacement field.

##  II Real space screened Coulomb potential

A cut along the $x$-direction of the screened Coulomb potential in real space
is shown in Fig. S2. For this, we compute the Fourier transform of the
potential via,

$V\left(\mathbf{r}\right)=\Delta_{\mathbf{k}}^{2}\sum_{\mathbf{k}}V_{scr}\left(\mathbf{k}\right)e^{-\mathrm{i}\mathbf{k}\cdot\mathbf{r}},$
(S3)

where $\Delta_{\mathbf{k}}$ is the spacing between adjacent points in
$\mathbf{k}$-space.

Figure S2: Real space cut along the $x$-direction of the screened Coulomb
potential. (a) in BBG near the hole-doped vHs at $T=10$ mK. (b) in RTG near
the hole-doped vHs at $T=65$ mK.

As can be seen in Fig. S2, the repulsive Coulomb potential develops into an
effective attractive potential once the screening is taken into consideration
at a temperature around the $T_{c}$ of BBG and RTG. Since the continuum model
yields the correct electronic properties up to a cutoff in the reciprocal
space defined by $k_{\Lambda}$, we think that the position of the minimum of
the screened potential is overestimated as the short-range information cannot
be well captured with this model, and the magnitude is also unreliable.
However, it is good enough to qualitatively describe the appearance of an
attractive potential region which promotes superconductivity.

##  III Ising spin-orbit coupling

To account for the effects induced in the BBG by the proximity of WSe2 and
following the ideas introduced in [4], we include Ising spin-orbit coupling
(SOC) to check if it can give rise to an enhancement of the critical
temperature as measured experimentally. The Ising SOC results in an energy
imbalance between Cooper pairs
$\ket{\mathbf{K}^{+},\uparrow;\mathbf{K}^{-},\downarrow}$ and
$\ket{\mathbf{K}^{+},\downarrow;\mathbf{K}^{-},\uparrow}$, making one of them
to rise in energy by a quantity given by the SOC strength parameter,
$\lambda_{I}$, while lowering the other by the same amount. Thus, the spin-
valley-dependent hamiltonian for Ising SOC is given by,

$\mathcal{H}_{I,\xi,s}=s\xi\lambda_{I}\mathbb{I},$ (S4)

where $\mathbb{I}$ is the identity matrix. As can be inferred from Eq. (S4)
and seen in Figure S4, Ising SOC locks the spin and valley degrees of freedom
in the OP, making inapplicable the classification as spin-singlet or spin-
triplet.

In the presence of SOC, the original vHs at $\mu_{0}$, which is four-fold
degenerate, splits into two vHs which are two-fold degenerate, separated by
$2\lambda_{I}$. A practical consequence is that the susceptibility in Eq. 3 in
the main text needs to be separated in two parts, each with a prefactor of
$N_{f}=2$, which are then summed to give the total susceptibility. Regarding
the screening of the Coulomb potential, in reciprocal space, increasing the
strength of Ising SOC narrows the region around the Dirac point where the
minimum of the screened potential is concentrated, and makes it grow faster as
a function of momentum. In real space this directly translates into a more
attractive potential, see Fig. S3, which is consistent with the increment in
critical temperature shown in Fig. 3 in the main text. The saturation of the
critical temperature corresponds to the saturation of the minimum of the
screened potential in real space, as reflected in Table I.

Figure S3: (a) Screened potential in BBG in reciprocal space near the hole-doped vHs at $\mu=\mu_{0}+\lambda_{I}$, at $T_{c}$ given in Fig. 3 in the main text, for $\lambda_{I}=0,1,2,3,4,5$ meV, from left to right. (b) Real space screened potential for $\lambda_{I}=0$ meV (red), $\lambda_{I}=1$ meV (blue), $\lambda_{I}=2$ meV (green), $\lambda_{I}=3$ meV (cyan), $\lambda_{I}=4$ meV (magenta), $\lambda_{I}=5$ meV (orange). (c) Inset of (b) near the minimum of the screened potential. $\lambda_{I}$ (meV) |  $r_{min}$ (nm) |  $V_{scr}(r_{min})$ (eV)
---|---|---
0.0 | 12.1985 | -0.199664
1.0 | 12.2834 | -0.210588
2.0 | 12.3683 | -0.220475
3.0 | 12.4532 | -0.232751
4.0 | 12.5381 | -0.251644
5.0 | 12.5381 | -0.251644
Table 1: Relation between Ising SOC strength, $\lambda_{I}$, the minimum of
the screened potential and its position in real space near the hole-doped vHs
of BBG at the critical temperatures shown in Fig. 3 in the main text.

As for the results, the splitting induces two superconducting peaks as a
function of Fermi energy, as shown in Fig. S5 and seen in the experiments of
Ref. [4]. When the Fermi energy is placed at the $\mu_{0}+\lambda_{I}$ vHs,
the vHs at $\mu_{0}-\lambda_{I}$ is nearly or completely unpopulated,
depending on the value of $\lambda_{I}$, so the non-superconducting state is
nearly or completely polarized to a half-metal. In this case,
superconductivity is notably enhanced. In contrast, when $\mu$ is at the
$\mu_{0}-\lambda_{I}$ vHs, the flavours at the $\mu_{0}+\lambda_{I}$ vHs have
large Fermi surfaces. This situation leads to smaller critical temperatures
than without spin-orbit, for any reasonable value of $\lambda_{I}$, which
implies that the large Fermi surfaces away from $\mu$ have a detrimental
effect on superconductivity.

Figure S4: Order parameter of hole doped BBG with Ising SOC $\lambda_{I}=1$
meV, at $\mu=\mu_{0}+\lambda_{I}$, in (a) in valley K+ and (b) in valley K-,
(c,d) same as (a,b), at $\mu=\mu_{0}-\lambda_{I}$, where $\mu_{0}$ is the
Fermi energy at which the vHs appears in absence of Ising SOC. Figure S5:
Critical temperature in hole-doped BBG with Ising spin-orbit coupling, at both
vHs. Blue data correspond to the vHs at $\mu=\mu_{0}+\lambda_{I}$, while red
data correspond to the vHs at $\mu=\mu_{0}-\lambda_{I}$, where $\mu_{0}$ is
the Fermi energy at which the vHs appears in absence of Ising SOC.

Finally, it is also possible to consider a superconducting wavefunction
involving only Cooper pairs $\ket{K^{+}\uparrow,K^{-}\downarrow}$, which
interact with themselves via spin-valley flips. In this case, we find the same
critical temperatures as shown in Fig. 3 in the main text for the vHs at
$\mu_{0}+\lambda_{I}$.

##  IV Superconductivity in BBG with electron doping

As can be seen in Fig. 2 in the main text, the critical temperature peaks
close to the vHs in the DOS, which suggests that superconductivity could also
exist near the vHs of the conduction band. In order to confirm this we compute
the critical temperature around this vHs with the same procedure. The results
can be found in Fig. S6. We find a $T_{c}\approx 15$ mK near the vHs of the
conduction band in BBG. After including the short-range interaction, the
critical temperature increases to $T_{c}\approx 25$ mK. The OP is
qualitatively similar to the one with hole-doping.

Figure S6: (a) Critical temperature (blue) and DOS (a.u., dashed black) versus
Fermi energy in electron-doped BBG with a field-induced gap of 98 meV. The
insets show the shape of the Fermi surface. (b) Order parameter at the energy
for which $T_{c}=15$ mK, in valley $\mathbf{K}^{+}$ and (c) valley
$\mathbf{K}^{-}$. Again a sign change between valleys is visible, thus
pointing to the valley-singlet, spin-triplet nature of the superconducting
state in BBG. The maximum $T_{c}$ is found near the vHs, with $n_{e}\approx
0.385\cdot 10^{12}$ cm-2.

##  V Diagrams

Figure S7: Diagrams of the (a) long-range interaction mediated by the screened
Coulomb potential $V_{scr}$, which does not change the valley of the incoming
states. (b) Short-range interaction mediated by the momentum independent
Hubbard $U$, which makes the incoming states exchange valley.

Figure S7 shows the diagrams involved in superconductivity. We assume that
Cooper pairs are formed between electrons in different valleys. Diagram (a)
represents the long-range interaction, mediated by the screened Coulomb
potential. After scattering, each carrier remains in its original valley. The
kernel $\Gamma^{+}$ corresponds to this diagram, while $\Gamma^{-}$
corresponds to a diagram with $K^{+}$ and $K^{-}$ exchanged. The global
wavefunction is
$\ket{\Psi}\propto\ket{\mathbf{K}^{+},\uparrow;\mathbf{K}^{-},\downarrow}\pm\ket{\mathbf{K}^{+},\downarrow;\mathbf{K}^{-},\uparrow}$.
The ‘+’ combination corresponds to valley-singlet, spin-triplet pairing, while
the ‘–’ combination results in valley-triplet, spin-singlet pairing. In the
absence of short-range interactions, these two pairings are degenerate.
Diagram (b) represents the short-range interaction of the model, which is
mediated by a repulsive, momentum-independent Hubbard $U=3$ eV. After
scattering, the carriers exchange valley. This interaction can break the
degeneracy between the two pairings mentioned above.

Focusing on the superconducting OP in Fig. 4 in the main text, one might
wonder why the OP in each panel lives in the Fermi surface of only one valley,
even though one panel corresponds to $\tilde{\Delta}^{+}$ and the other to
$\tilde{\Delta}^{-}$, each of which stems from diagrams involving both
valleys. The reason for this is the application of time reversal symmetry
(TRS). The two matrix elements involved in diagram (a) can be expressed in
terms of quantities in valley $K^{+}$ only, after applying TRS. An analogous
argument applies to the panel that plots the $\tilde{\Delta}^{-}$ part of the
eigenvector. This time, we use TRS to express all quantities in the valley
$K^{-}$. The short-range interaction is a small perturbation that does not
change this picture.

##  VI The gap equation

We construct the gap equation following similar steps as if we were working
with the BCS theory but instead of setting an electron-phonon coupling as
mediator of the interaction between the pairs, we assume that the electron-
electron interaction alone leads to pairing, which essentially is the Kohn-
Luttinger methodology in [28]. This purely electronic interaction is the
screened long-range Coulomb potential that a carrier feels within the system,
which ultimately we consider to be the glue of the pairs. As there are many
carriers playing a role in electronic interactions we account for the
screening of the potential via the Random Phase Approximation (RPA).

Now we proceed to derive Eqs. (6, 7) in the main text. We start from the gap
equation (4), which includes two terms, one for the screened long-range
Coulomb potential $V_{scr}$ and another for the short-range Hubbard repulsion
$U$,

$\tilde{\Delta}^{+,i,j}\left(\mathbf{k}\right)=-\frac{K_{B}T}{\Omega}\sum_{\mathbf{q},\omega}\sum_{i^{\prime},j^{\prime}}\left[V_{scr}\left(\mathbf{k}-\mathbf{q}\right)G_{K^{+}}^{i,i^{\prime}}\left(\mathbf{q},\mathrm{i}\omega\right)G_{K^{-}}^{j,j^{\prime}}\left(-\mathbf{q},-\mathrm{i}\omega\right)\tilde{\Delta}^{+,i^{\prime},j^{\prime}}\left(\mathbf{q}\right)+UG_{K^{-}}^{i,i^{\prime}}\left(\mathbf{q},\mathrm{i}\omega\right)G_{K^{+}}^{j,j^{\prime}}\left(-\mathbf{q},-\mathrm{i}\omega\right)\tilde{\Delta}^{-,j^{\prime},i^{\prime}}\left(-\mathbf{q}\right)\right].$
(S5)

We can work separately with the long-range and short-range terms and then
simply develop the equations. We begin with the kernel of the long-range
interaction and substitute the electrons’ Green’s functions by their
definitions,
$G_{\mathbf{K}^{\xi}}^{i,i^{\prime}}\left(\mathbf{k},\pm\mathrm{i}\omega\right)=\sum_{n}\frac{\Psi^{i}_{\xi,n}\left(\mathbf{k}\right)\Psi^{i^{\prime},*}_{\xi,n}\left(\mathbf{k}\right)}{\mathrm{i}\omega\mp\epsilon_{\mathbf{K}^{\xi},n}\left(\mathbf{k}\right)}$,
leading to,

$\begin{split}\tilde{\Delta}^{+,i,j}_{n,m}\left(\mathbf{k}\right)&=-\frac{K_{B}T}{\Omega}\sum_{\mathbf{q},\omega}\sum_{i^{\prime},j^{\prime}}V_{scr}\left(\mathbf{k}-\mathbf{q}\right)G_{K^{+}}^{i,i^{\prime}}\left(\mathbf{q},\mathrm{i}\omega\right)G_{K^{-}}^{j,j^{\prime}}\left(-\mathbf{q},-\mathrm{i}\omega\right)\tilde{\Delta}^{+,i^{\prime},j^{\prime}}\left(\mathbf{q}\right),\\\
&=-\frac{K_{B}T}{\Omega}\sum_{\mathbf{q},\omega}\sum_{i^{\prime},j^{\prime}}V_{scr}\left(\mathbf{k}-\mathbf{q}\right)\frac{\Psi^{i}_{K^{+},n}\left(\mathbf{q}\right)\Psi^{i^{\prime},*}_{K^{+},n}\left(\mathbf{q}\right)}{\mathrm{i}\omega-\epsilon_{K^{+},n}\left(\mathbf{q}\right)}\frac{\Psi^{j}_{K^{-},m}\left(-\mathbf{q}\right)\Psi^{j^{\prime},*}_{K^{-},m}\left(-\mathbf{q}\right)}{\mathrm{i}\omega+\epsilon_{K^{-},m}\left(-\mathbf{q}\right)}\tilde{\Delta}^{+,i^{\prime},j^{\prime}}\left(\mathbf{q}\right),\\\
&=-\frac{K_{B}T}{\Omega}\sum_{\mathbf{q}}V_{scr}\left(\mathbf{k}-\mathbf{q}\right)\Psi^{i}_{K^{+},n}\left(\mathbf{q}\right)\Psi^{j}_{K^{-},m}\left(-\mathbf{q}\right)\sum_{\omega}\frac{1}{\mathrm{i}\omega-\epsilon_{K^{+},n}\left(\mathbf{q}\right)}\frac{1}{\mathrm{i}\omega+\epsilon_{K^{-},m}\left(-\mathbf{q}\right)}\sum_{i^{\prime},j^{\prime}}\Psi^{i^{\prime},*}_{K^{+},n}\left(\mathbf{q}\right)\Psi^{j^{\prime},*}_{K^{-},m}\left(-\mathbf{q}\right)\tilde{\Delta}^{+,i^{\prime},j^{\prime}}\left(\mathbf{q}\right),\\\
&=-\frac{K_{B}T}{\Omega}\sum_{\mathbf{q}}V_{scr}\left(\mathbf{k}-\mathbf{q}\right)\Psi^{i}_{n}\left(\mathbf{q}\right)\Psi^{j,*}_{m}\left(\mathbf{q}\right)\sum_{\omega}\frac{1}{\mathrm{i}\omega-\epsilon_{n}\left(\mathbf{q}\right)}\frac{1}{\mathrm{i}\omega+\epsilon_{m}\left(\mathbf{q}\right)}\sum_{i^{\prime},j^{\prime}}\Psi^{i^{\prime},*}_{n}\left(\mathbf{q}\right)\Psi^{j^{\prime}}_{m}\left(\mathbf{q}\right)\tilde{\Delta}^{+,i^{\prime},j^{\prime}}\left(\mathbf{q}\right),\\\
&=-\frac{1}{\Omega}\sum_{\mathbf{q}}V_{scr}\left(\mathbf{k}-\mathbf{q}\right)\Psi^{i}_{n}\left(\mathbf{q}\right)\Psi^{j,*}_{m}\left(\mathbf{q}\right)\frac{f(-\epsilon_{m}\left(\mathbf{q}\right))-f(\epsilon_{n}\left(\mathbf{q}\right))}{\epsilon_{n}\left(\mathbf{q}\right)+\epsilon_{m}\left(\mathbf{q}\right)}\tilde{\Delta}_{n,m}^{+}\left(\mathbf{q}\right),\end{split}$
(S6)

where in the last step we have performed a Matsubara sum using
$\sum_{\mathrm{i}\omega}\frac{1}{\mathrm{i}\omega-\epsilon_{1}}\frac{1}{\mathrm{i}\omega+\epsilon_{2}}=\frac{1}{K_{B}T}\frac{1}{\epsilon_{1}+\epsilon_{2}}\left[f(-\epsilon_{2})-f(\epsilon_{1})\right]$,
thus,

$\sum_{i,j}\Psi^{i,*}_{n}\left(\mathbf{k}\right)\Psi^{j}_{m}\left(\mathbf{k}\right)\tilde{\Delta}^{+,i,j}_{n,m}\left(\mathbf{k}\right)=-\frac{1}{\Omega}\sum_{\mathbf{q}}V_{scr}\left(\mathbf{k}-\mathbf{q}\right)\frac{f(-\epsilon_{m}\left(\mathbf{q}\right))-f(\epsilon_{n}\left(\mathbf{q}\right))}{\epsilon_{n}\left(\mathbf{q}\right)+\epsilon_{m}\left(\mathbf{q}\right)}\sum_{i,j}\Psi^{i,*}_{n}\left(\mathbf{k}\right)\Psi^{j}_{m}\left(\mathbf{k}\right)\Psi^{i}_{n}\left(\mathbf{q}\right)\Psi^{j,*}_{m}\left(\mathbf{q}\right)\tilde{\Delta}_{n,m}^{+}\left(\mathbf{q}\right).$
(S7)

Here we introduce the new quantity
$\tilde{\tilde{\Delta}}^{+}\left(\mathbf{p}\right)=\tilde{\Delta}^{+}\left(\mathbf{p}\right)\sqrt{\frac{f(-\epsilon_{m}\left(\mathbf{p}\right))-f(\epsilon_{n}\left(\mathbf{p}\right))}{\epsilon_{n}\left(\mathbf{p}\right)+\epsilon_{m}\left(\mathbf{p}\right)}}$
for the kernel to be hermitian, so after the substitution,

$\tilde{\tilde{\Delta}}_{n,m}^{+}\left(\mathbf{k}\right)=-\frac{1}{\Omega}\sum_{\mathbf{q}}V_{scr}\left(\mathbf{k}-\mathbf{q}\right)\sqrt{\frac{f(-\epsilon_{m}\left(\mathbf{k}\right))-f(\epsilon_{n}\left(\mathbf{k}\right))}{\epsilon_{n}\left(\mathbf{k}\right)+\epsilon_{m}\left(\mathbf{k}\right)}}\sqrt{\frac{f(-\epsilon_{m}\left(\mathbf{q}\right))-f(\epsilon_{n}\left(\mathbf{q}\right))}{\epsilon_{n}\left(\mathbf{q}\right)+\epsilon_{m}\left(\mathbf{q}\right)}}\sum_{i,j}\Psi^{i,*}_{n}\left(\mathbf{k}\right)\Psi^{j}_{m}\left(\mathbf{k}\right)\Psi^{i}_{n}\left(\mathbf{q}\right)\Psi^{j,*}_{m}\left(\mathbf{q}\right)\tilde{\Delta}^{+}_{n,m}\left(\mathbf{q}\right).$
(S8)

Now, following the same procedure we continue with the kernel of the short-
range interaction.

$\begin{split}\tilde{\Delta}^{+,i,j}_{n,m}\left(\mathbf{k}\right)&=-U\frac{K_{B}T}{\Omega}\sum_{\mathbf{q},\omega}\sum_{i^{\prime},j^{\prime}}G_{K^{-}}^{i,i^{\prime}}\left(\mathbf{q},\mathrm{i}\omega\right)G_{K^{+}}^{j,j^{\prime}}\left(-\mathbf{q},-\mathrm{i}\omega\right)\tilde{\Delta}^{-,j^{\prime},i^{\prime}}\left(-\mathbf{q}\right),\\\
&=-U\frac{K_{B}T}{\Omega}\sum_{\mathbf{q},\omega}\sum_{i^{\prime},j^{\prime}}\frac{\Psi^{i}_{K^{-},n}\left(\mathbf{q}\right)\Psi^{i^{\prime},*}_{K^{-},n}\left(\mathbf{q}\right)}{\mathrm{i}\omega-\epsilon_{K^{-},n}\left(\mathbf{q}\right)}\frac{\Psi^{j}_{K^{+},m}\left(-\mathbf{q}\right)\Psi^{j^{\prime},*}_{K^{+},m}\left(-\mathbf{q}\right)}{\mathrm{i}\omega+\epsilon_{K^{+},m}\left(-\mathbf{q}\right)}\tilde{\Delta}^{-,j^{\prime},i^{\prime}}\left(-\mathbf{q}\right),\\\
&=-U\frac{K_{B}T}{\Omega}\sum_{\mathbf{q}}\Psi^{i}_{K^{-},n}\left(\mathbf{q}\right)\Psi^{j}_{K^{+},m}\left(-\mathbf{q}\right)\sum_{\omega}\frac{1}{\mathrm{i}\omega-\epsilon_{K^{-},n}\left(\mathbf{q}\right)}\frac{1}{\mathrm{i}\omega+\epsilon_{K^{+},m}\left(-\mathbf{q}\right)}\sum_{i^{\prime},j^{\prime}}\Psi^{i^{\prime},*}_{K^{-},n}\left(\mathbf{q}\right)\Psi^{j^{\prime},*}_{K^{+},m}\left(-\mathbf{q}\right)\tilde{\Delta}^{-,j^{\prime},i^{\prime}}\left(-\mathbf{q}\right),\\\
&=-U\frac{K_{B}T}{\Omega}\sum_{\mathbf{q}}\Psi^{i,*}_{n}\left(-\mathbf{q}\right)\Psi^{j}_{m}\left(-\mathbf{q}\right)\sum_{\omega}\frac{1}{\mathrm{i}\omega-\epsilon_{n}\left(-\mathbf{q}\right)}\frac{1}{\mathrm{i}\omega+\epsilon_{m}\left(-\mathbf{q}\right)}\sum_{i^{\prime},j^{\prime}}\Psi^{i^{\prime}}_{n}\left(-\mathbf{q}\right)\Psi^{j^{\prime},*}_{m}\left(-\mathbf{q}\right)\tilde{\Delta}^{-,j^{\prime},i^{\prime}}\left(-\mathbf{q}\right),\\\
&=-U\frac{K_{B}T}{\Omega}\sum_{\mathbf{q}}\Psi^{i,*}_{n}\left(\mathbf{q}\right)\Psi^{j}_{m}\left(\mathbf{q}\right)\sum_{\omega}\frac{1}{\mathrm{i}\omega-\epsilon_{n}\left(\mathbf{q}\right)}\frac{1}{\mathrm{i}\omega+\epsilon_{m}\left(\mathbf{q}\right)}\sum_{i^{\prime},j^{\prime}}\Psi^{i^{\prime}}_{n}\left(\mathbf{q}\right)\Psi^{j^{\prime},*}_{m}\left(\mathbf{q}\right)\tilde{\Delta}^{-,j^{\prime},i^{\prime}}\left(\mathbf{q}\right),\\\
&=-U\frac{K_{B}T}{\Omega}\sum_{\mathbf{q}}\Psi^{i,*}_{n}\left(\mathbf{q}\right)\Psi^{j}_{m}\left(\mathbf{q}\right)\sum_{\omega}\frac{1}{\mathrm{i}\omega-\epsilon_{n}\left(\mathbf{q}\right)}\frac{1}{\mathrm{i}\omega+\epsilon_{m}\left(\mathbf{q}\right)}\sum_{i^{\prime},j^{\prime}}\Psi^{i^{\prime},*}_{m}\left(\mathbf{q}\right)\Psi^{j^{\prime}}_{n}\left(\mathbf{q}\right)\tilde{\Delta}^{-,i^{\prime},j^{\prime}}\left(\mathbf{q}\right),\\\
&=-U\frac{1}{\Omega}\sum_{\mathbf{q}}\Psi^{i,*}_{n}\left(\mathbf{q}\right)\Psi^{j}_{m}\left(\mathbf{q}\right)\frac{f(-\epsilon_{m}\left(\mathbf{q}\right))-f(\epsilon_{n}\left(\mathbf{q}\right))}{\epsilon_{n}\left(\mathbf{q}\right)+\epsilon_{m}\left(\mathbf{q}\right)}\tilde{\Delta}_{n,m}^{-}\left(\mathbf{q}\right),\end{split}$
(S9)

thus,

$\sum_{i,j}\Psi^{i,*}_{n}\left(\mathbf{k}\right)\Psi^{j}_{m}\left(\mathbf{k}\right)\tilde{\Delta}^{+,i,j}_{n,m}\left(\mathbf{k}\right)=-U\frac{1}{\Omega}\sum_{\mathbf{q}}\frac{f(-\epsilon_{m}\left(\mathbf{q}\right))-f(\epsilon_{n}\left(\mathbf{q}\right))}{\epsilon_{n}\left(\mathbf{q}\right)+\epsilon_{m}\left(\mathbf{q}\right)}\sum_{i,j}\Psi^{i,*}_{n}\left(\mathbf{k}\right)\Psi^{j}_{m}\left(\mathbf{k}\right)\Psi^{i,*}_{n}\left(\mathbf{q}\right)\Psi^{j}_{m}\left(\mathbf{q}\right)\tilde{\Delta}_{n,m}^{-}\left(\mathbf{q}\right),$
(S10)

and we arrive at,

$\tilde{\Delta}_{n,m}^{+}\left(\mathbf{k}\right)=-U\frac{1}{\Omega}\sum_{\mathbf{q}}\sqrt{\frac{f(-\epsilon_{m}\left(\mathbf{k}\right))-f(\epsilon_{n}\left(\mathbf{k}\right))}{\epsilon_{n}\left(\mathbf{k}\right)+\epsilon_{m}\left(\mathbf{k}\right)}}\sqrt{\frac{f(-\epsilon_{m}\left(\mathbf{q}\right))-f(\epsilon_{n}\left(\mathbf{q}\right))}{\epsilon_{n}\left(\mathbf{q}\right)+\epsilon_{m}\left(\mathbf{q}\right)}}\sum_{i,j}\Psi^{i,*}_{n}\left(\mathbf{k}\right)\Psi^{j}_{m}\left(\mathbf{k}\right)\Psi^{i,*}_{n}\left(\mathbf{q}\right)\Psi^{j}_{m}\left(\mathbf{q}\right)\tilde{\Delta}_{n,m}^{-}\left(\mathbf{q}\right).$
(S11)

Combining Eqs. (S8, S11) along with time reversal symmetry gives Eqs. (6, 7)
in the main text.

##  VII Critical temperature dependence on electric field.

Figure S8: Electric field dependence of the critical temperature near the vHs
of hole-doped (a) BBG and (b) RTG. Yellow boxes label the electronic density
at which the critical temperature has been computed for each case, in units of
$10^{12}$ cm-2

Figure S8 indicates the critical temperatures of BBG and RTG, as a function of
the external electric field. $T_{c}$ peaks at a certain field, and tends to
diminish when the field changes. These results are in qualitative agreement
with the experimental findings in Refs. [9, 1].

##  VIII The role of the Hubbard U.

Figure S9: Symmetry of the OP, as a function of the short-range interaction
strength, $U$, for the four largest eigenvalues near the hole-doped vHs of (a)
BBG at $T=10$ mK and (b) RTG at $T=85$mK. These plots are meant to illustrate
the form of the order parameter as a function of U, but note that the
eigenvalues $\lambda<1$ imply that the superconducting state has not been
reached here, and details change once it is, e.g. the degeneracy and
associated $C_{3}$-symmetry breaking of the largest eigenvalue of BBG at $U=3$
eV disappears.

We conclude with a tentative discussion on the effects of the value of $U$,
which is the least well determined parameter in our model. Realistic values
$U\sim 2-4$ eV. [59, 60]. We recall here that the screened long-range Coulomb
interaction leads to a degeneracy between valley-triplet, spin-singlet and
valley-singlet, spin-triplet pairings. The short-range Hubbard-like
interaction can break this degeneracy, selecting a preferred pairing symmetry.
For the value used in this work, $U=3$ eV, the short-range interaction can be
seen as a small perturbation. In BBG and RTG, an infinitesimally small value
of $U$ is enough to break the aforementioned degeneracy. To first order in
perturbation theory, the change in the kernel eigenvalue due to the short-
range interaction is:

$\delta\epsilon=-\sum_{\mathbf{k},\mathbf{k}^{\prime}}\tilde{\Delta}_{\mathbf{k}}\tilde{U}(\mathbf{k},\mathbf{k}^{\prime})\tilde{\Delta}_{\mathbf{k}^{\prime}},$
(S12)

where $\tilde{\Delta}_{\mathbf{k}}$ and $\tilde{\Delta}_{\mathbf{k}^{\prime}}$
are the order parameters in valley $K^{+}$ and $K^{-}$, in the absence of the
short-range interaction. From this follows the common wisdom that a repulsive,
i.e. positive, U, favours valley-singlet pairing, in which the OP changes sign
between valleys, so that $\delta\epsilon>0$, resulting in an enhanced critical
temperature. However, Eq. S12 is valid only for single pocket superconductors,
while in BBG and RTG each valley shows at least two pockets, which makes the
situation more complex, as pockets with a certain sign in one valley can
interact strongly with others with the opposite sign in the other valley, even
if the global state is valley-triplet, spin-singlet. Still, both BBG and RTG
have valley-singlet, spin-triplet pairing for any finite value of U, which is
the usual result with repulsive short-range interactions, as mentioned above.

It is illustrative to note that in the limit of large U, the short-range
interaction is not a perturbation anymore, but rather becomes the leading
interaction and overshadows the long-range Coulomb potential. When short-range
interactions are dominant, the most energetically favorable solution is the
one in which each valley as a constant sign, and there is a change of sign
between valleys. Thus, for large values of $U$ the order parameter losses the
sign changes within each valley. This phenomenon occurs in both materials,
although in the bilayer the intravalley signs only disappear completely at
$U\gtrsim 25$ eV.

##  IX Critical temperature starting from a Half Metal state.

Figure S10 shows the critical temperature as a function of Fermi energy in
both materials, when the non-superconducting state is a half metal. The
calculation is equivalent to the one in Fig. 2 in the main text, but the
factor of $N_{f}=4$ in the susceptibility is substituted by $N_{f}=2$. The
situation is qualitatively similar to the one with Ising SOC. However, the
complete disappearance of two flavours leads to a more severe change in the
susceptibility than the energy split induced by Ising SOC. It is an
approximation which gives a rough idea of what happens for very large values
of $\lambda_{I}$. The critical temperature we find here for BBG is 80 mK,
while starting from the full metal it is 10 mK, and in the experiment, 26 mK
[1]. For RTG. the critical temperatures starting from a half metal or a full
metal state are 110 mK and 65 mK, respectively, while the experiment reports
105 mK [9].

Figure S10: Superconducting critical temperature starting from a half metal
non-superconducting state (solid blue) and DOS (a.u. dashed black) versus
Fermi energy in (a) BBG with hole doping and a gap of $98$ meV; near the vHs
we observe a $T_{c}\approx 80$ mK, and in (b) RTG with hole doping and a gap
of $74$ meV; near the vHs we observe a $T_{c}\approx 110$ mK.