# Junctions and superconducting symmetry in twisted bilayer graphene

Héctor Sainz-Cruz<EMAIL_ADDRESS>$Imdea\ Nanoscience,\ Faraday\ 9,\
28015\ Madrid,\ Spain$ Pierre A. Pantaleón  $Imdea\ Nanoscience,\ Faraday\ 9,\
28015\ Madrid,\ Spain$ Võ Ti\́ecircumflexn Phong  $Department\ of\ Physics\
and\ Astronomy,\ University\ of\ Pennsylvania,\ Philadelphia\ PA\ 19104$
Alejandro Jimeno-Pozo  $Imdea\ Nanoscience,\ Faraday\ 9,\ 28015\ Madrid,\
Spain$ Francisco Guinea  $Imdea\ Nanoscience,\ Faraday\ 9,\ 28015\ Madrid,\
Spain$ $Donostia\ International\ Physics\ Center,\ Paseo\ Manuel\ de\
Lardizabal\ 4,\ 20018\ San\ Sebastian,\ Spain$

(27th August 2024)

###### Abstract

Junctions provide a wealth of information on the symmetry of the order
parameter of superconductors. We analyze junctions between a scanning
tunneling microscope (STM) tip and superconducting twisted bilayer graphene
(TBG) and TBG Josephson junctions (JJs). We compare superconducting phases
that are even or odd under valley exchange ($s$\- or $f$-wave). The critical
current in mixed ($s$\- and $f$-) JJs strongly depends on the angle between
the junction and the lattice. In STM-TBG junctions, due to Andreev reflection,
$f$-wave leads to a prominent peak in subgap conductance, as seen in
experiments.

Introduction. Graphene multilayers host a myriad of exotic correlated and
topological phases Cao  _et al._ (2018a); Polshyn  _et al._ (2019); Sharpe
_et al._ (2019); Serlin  _et al._ (2020); Chen  _et al._ (2020); Uri  _et al._
(2020); Saito  _et al._ (2020); Wong  _et al._ (2020); Zondiner  _et al._
(2020); Stepanov  _et al._ (2020); Xu  _et al._ (2020); Choi  _et al._ (2021);
Rozen  _et al._ (2021); Cao  _et al._ (2021a); Stepanov  _et al._ (2021); Oh
_et al._ (2021); Xie  _et al._ (2021); Berdyugin  _et al._ (2022); Turkel  _et
al._ (2022); Huang  _et al._ (2022); de la Barrera  _et al._ (2022); Kim  _et
al._ (2022); Seiler  _et al._ (2022). Perhaps most interesting and enigmatic
among them is superconductivity, possibly with unconventional pairing
symmetries and mechanisms, observed in alternating-twist stacks of up to five
layers Cao  _et al._ (2018b); Yankowitz  _et al._ (2019); Lu  _et al._ (2019);
Park  _et al._ (2021); Hao  _et al._ (2021); Park  _et al._ (2022); Zhang  _et
al._ (2022) and in Bernal bilayers and rhombohedral trilayers Zhou  _et al._
(2021, 2022); Zhang  _et al._ (2023); Holleis  _et al._ (2023). Crucially, the
observed superconductivity violates the Pauli limit for spin-singlet pairing
Cao  _et al._ (2021b); Park  _et al._ (2022); Zhou  _et al._ (2021, 2022);
Zhang  _et al._ (2023); Holleis  _et al._ (2023) and has been observed in
settings that break time-reversal symmetry (TRS) Lin  _et al._ (2022),
strongly suggesting a spin-triplet pairing in these materials. However, the
pairing may be a mixture of singlet and triplet Lake  _et al._ (2022), and the
exact symmetries involved ($s$-, $p$-, $d$\- and/or $f$-) are still unknown
despite intense theoretical and experimental efforts to uncover them.

Recently, several experiments have studied these unconventional
superconducting states using transport measurements: either with a scanning
tunneling microscope (STM) tip Oh  _et al._ (2021); Kim  _et al._ (2022) or
with Josephson junctions (JJs) de Vries  _et al._ (2021); Rodan-Legrain  _et
al._ (2021); Diez-Merida  _et al._ (2021); Xie  _et al._ (2023); Hu  _et al._
(2022) and Superconducting Quantum Interference Devices (SQUIDs) Portolés  _et
al._ (2022). In the former setup, by comparing the transmission between the
STM tip and the superconducting surface in the weak and strong-coupling
regimes, one can gain important insights about the symmetry of the order
parameter. For instance, the experimental observations, such as the peak in
the subgap conductance Oh  _et al._ (2021); Kim  _et al._ (2022), seem
inconsistent with $s$-, $p$\- and $d$-wave pairings Lake  _et al._ (2022);
Sukhachov  _et al._ (2023). In the latter setups, the overlap of the
superconductors' wavefunctions at the junction's link gives rise to a zero-
frequency supercurrent whose magnitude and superconducting phase-dependence
carry characteristics of the pairing symmetry Josephson (1962); Sigrist and
Ueda (1991); Golubov  _et al._ (2004).

Building on these experimental insights, we argue in this Letter that
transport measurements in junctions are ideal probes of the pairing symmetry
in twisted graphene superconductors, similar to the elucidation of $d$-wave
pairing in cuprate superconductors Tsuei  _et al._ (1994); Tsuei and Kirtley
(2000), and that existing STM data Oh  _et al._ (2021); Kim  _et al._ (2022)
are consistent with $f$-wave pairing. The Fermi surface of these graphene-
based systems contains two valleys. We consider superconducting order
parameters that are either even or odd under valley exchange, which in the
absence of spin-orbit coupling correspond to spin-singlet $s$-wave
superconductivity or spin-triplet $f$-wave superconductivity respectively. In
`mixed' Josephson junctions connecting a $s$-wave to a $f$-wave
superconductor, we observe that the critical current dramatically depends on
the angle between the junction and the graphene lattice axis. Therefore,
Josephson junctions are useful for determining whether two superconducting
phases differ in their valley exchange parity.

In the STM-superconductor junction, we find that the subgap conductance shows
a prominent zero-bias peak for $f$-wave pairing only, due to enhanced Andreev
reflection. This peak has been observed in experiments on both twisted bilayer
Oh  _et al._ (2021) and twisted trilayer graphene Kim  _et al._ (2022). This
result puts forward $f$-wave pairing as a leading candidate for the
superconducting symmetry of twisted bilayer graphene, which is also consistent
with previous theoretical models based on Coulomb-interaction-mediated Cooper
pairing Crépel  _et al._ (2022); Cea and Guinea (2021).

Figure 1: STM tip-superconducting TBG junction. (a) Sketch of the model for
the junction used to calculate its transport properties. See text for details.
(b-e): normal and Andreev reflections (top) and total conductance (bottom) for
junctions with spin-singlet $s$-wave (dashed) and spin-triplet $f$-wave
pairing (solid), in the perfect contact limit. (b) With equal Fermi velocities
in all channels,
$t_{tip}=t_{sc,K}=t_{sc,K^{\prime}}=1,\,t_{tip,K}=t_{tip,K^{\prime}}=1/\sqrt{2},\,t_{K,K^{\prime}}=0,\,\Delta_{K}=0.05,\Delta_{K^{\prime}}=\pm
0.05$. (c) With a large Fermi velocity mismatch in the normal and
superconducting channels:
$t_{tip}=10,\,t_{tip,K}=t_{tip,K^{\prime}}=10/\sqrt{2}$, others as in (b). (d)
With Fermi velocity mismatch and intervalley scattering $t_{K,K^{\prime}}=1$,
others as in (c). (e) With Fermi velocity mismatch, intervalley scattering and
spin-orbit coupling: $\Delta_{K}=0.05,\Delta_{K^{\prime}}=\pm 0.02$, others as
in (d).

STM tip-superconducting TBG junction. The model. General features of transport
in normal-superconductor junctions are described in Ref. Blonder  _et al._
(1982). The coupling between the two electrodes is given by a scattering
matrix, determined by a dimensionless transmission amplitude, $T$. The model
has been extended in Lake  _et al._ (2022); Sukhachov  _et al._ (2023). As in
Ref. Blonder  _et al._ (1982), the normal metal tip and the superconducting
electrode are described in terms of incoming and outgoing single channels. On
the superconductor, the states in the channel are defined as suitable averages
in momentum space of the quasiparticles. The momentum dependence of the gap
leads to a momentum dependence of the mixing between electron and hole-like
states in the superconductor, and it modifies the transmission of the
junction, both in the tunneling and in the contact regimes. The Blonder-
Tinkham-Klapwijk (BTK) model Blonder  _et al._ (1982) has also been extended
to strongly coupled superconductors, where the chemical potential can be below
the bottom of the band Setiawan and Hofmann (2022); Lewandowski  _et al._
(2023).

We describe the metal$-$superconductor junction as one ingoing normal channel,
which represents the tip, and two outgoing superconducting channels, which
represent the two valleys in TBG. The signs of the gaps in these channels can
be equal, describing a spin-singlet $s$-wave superconductor, or opposite,
describing a spin-triplet $f$-wave superconductor aMi . The model can also be
applied to an Ising superconductor Holleis  _et al._ (2023) in a system with
strong spin-orbit coupling, characterized by spin-valley locked Cooper pairs
of the type $|K,\uparrow;K^{\prime},\downarrow\rangle$.

The three-channel model described above is discretized as a tight-binding
model, see Fig. 1(a). The normal channel is described by nearest-neighbor
hopping $t_{tip}$, which determines its Fermi velocity and density of states.
The superconducting channels are described by two nearest neighbor hoppings,
$t_{sc,K}$ and $t_{sc,K^{\prime}}$, and two gaps, $\Delta_{K}$ and
$\Delta_{K^{\prime}}$. The coupling between the normal channel and the two
superconducting channels is described by the hoppings $t_{tip,K}$ and
$t_{tip,K^{\prime}}$. Without loss of generality, we assume that the Fermi
energy is $\epsilon_{F}=0$, so that each channel has exact electron-hole
symmetry. Finally, we consider that the tip is a local perturbation which can
induce intervalley scattering, parametrized by another hopping,
$t_{K,K^{\prime}}$.

We solve the transmission of the junction by matching incoming and outgoing
waves in the three channels. If the energy $\epsilon$ is within the
superconducting gaps, we use evanescent waves in the superconducting channels.
For each energy, there are four propagating or evanescent waves in each
channel. We assume that there is an incoming wave of electron character and
amplitude 1 in the tip channel. In the same channel, there can be one electron
and one hole outgoing channels, describing normal and Andreev reflection, with
amplitudes $R_{N}$ and $R_{A}$, respectively. In each of the two
superconducting channels there can be two decaying evanescent waves, when the
energy is within the gap, or two outgoing propagating waves. We describe the
four amplitudes as $T_{i,j}$, where $i=K,K^{\prime}$ stands for the channel,
and $j=1,2$ stands for the wavefunction within each channel. The transport
properties of the junction are determined by these six amplitudes. The
conductance of the junction is $G=1-|R_{N}|^{2}+|R_{A}|^{2}$. The matching
conditions involve the amplitudes of the wavefunctions at the three sites
which describe the junction. The equations can be found in Ref. SM .

STM tip-superconducting TBG junction. Results. When the Fermi velocities in
all channels are equal, the tip channel merges smoothly into the even
combination of the $K$ and $K^{\prime}$ channels and the junction behaves as
described by the BTK theory in the regime of perfect contact, see Fig. 1(b).
For $s$-wave pairing, and at zero voltage, Andreev scattering leads to a
conductance twice as large as a single normal channel Blonder  _et al._
(1982). For $f$-wave pairing, negative interference between the two hole
channels cancels Andreev reflection. This cancellation can be expected
whenever the order parameter has a sign change between states related by TRS
Sukhachov  _et al._ (2023). At high voltages the conductance reduces to the
conductance of a single channel in both cases.

The bandwidth and Fermi velocity in TBG are considerably smaller than in a
normal metal. This Fermi velocity mismatch induces elastic back-scattering in
the normal phase, which reduces the conductance above the gap, see Fig. 1(c).
Subgap Andreev reflection for $s$-wave superconductivity is strongly
suppressed, and it remains zero for the $f$-wave phase, for a detailed
explanation see Ref. SM . The tip can also induce a perturbation on the
superconductor, on scales comparable to the atomic spacing. Such a
perturbation will induce intervalley scattering. Fig. 1(d) shows results
obtained for an intervalley coupling comparable to the bandwidth of the
superconductor. This perturbation can be considered as disorder which does not
violate TRS. The presence of intervalley scattering does not change
significantly the conductance of the junction in an $s$-wave superconductor,
in agreement with Anderson's theorem Anderson (1959). On the other hand, it is
a pair breaking perturbation in an $f$-wave superconductor, which induces
subgap states, see Ref. SM . These states allow for subgap Andreev reflection.
As a result, the subgap conductance of the junction is strongly enhanced by
intervalley scattering in a $f$-wave superconductor, leading to a zero bias
peak, highlighted in Fig. 1(d), that has been seen in the experiments of Refs.
Oh  _et al._ (2021); Kim  _et al._ (2022).

Recent transport experiments Zhang  _et al._ (2023); Holleis  _et al._ (2023)
reveal that proximity induced spin-orbit coupling promotes the superconducting
properties of Bernal bilayer graphene. An effect of spin-orbit coupling is to
break the equivalence between the Cooper pairs
$|K,\uparrow;K^{\prime},\downarrow\rangle$ and
$|K,\downarrow;K^{\prime}\uparrow\rangle$. In the model studied here, the
spin-orbit coupling makes the two channels inequivalent. Results are shown in
Fig. 1(e).

Josephson junctions. The model. For the study of JJs, our setup consists of a
TBG crystal, in which the electrodes are superconducting and the weak link is
in a normal metal or band insulating phase, as shown in Fig. 2(c). We start
from a tight-binding, non-interacting Hamiltonian ${\cal H}_{0}$ Lin and
Tománek (2018) that includes Hartree electron-electron interactions through an
electrostatic potential Guinea and Walet (2018); Rademaker  _et al._ (2019).
The parameters in the tight binding model are scaled, such that the central
bands of a TBG with twist angle $\theta$ are approximated by the central bands
of an equivalent lattice with twist angle $\lambda\theta$, with $\lambda>1$
Gonzalez-Arraga  _et al._ (2017); Vahedi  _et al._ (2021); Sainz-Cruz  _et
al._ (2021), see Fig. 2(a).

The critical current comes from second order perturbation theory and is the
derivative of the free energy $E$ with respect to the superconducting phase
difference $\phi$:

$\mathcal{I}=\frac{e}{h}\frac{\partial E}{\partial\phi}\,.$ (1)

To obtain the energies of the TBG junction, we diagonalize the Bogoliubov-de
Gennes Hamiltonian,

$\displaystyle{\cal H}_{BdG}\ket{\Psi}=\begin{pmatrix}{\cal
H}_{0}-\epsilon_{F}&f(\Delta)\\\ f^{\dagger}(\Delta)&\epsilon_{F}-{\cal
H}_{0}\\\ \end{pmatrix}\begin{pmatrix}\Psi_{e}\\\ \Psi_{h}\\\
\end{pmatrix}=E\begin{pmatrix}\Psi_{e}\\\ \Psi_{h}\\\ \end{pmatrix}\,,$ (2)

where $\epsilon_{F}$ is the Fermi energy. Again, we compare $s$-wave pairing,
which we model with a on-site attractive Hubbard term
$f(\Delta)=-\Delta_{S}\mathbb{1}$, and $f$-wave pairing, which results from
Haldane-like hoppings Haldane (1988); Pai  that allow an electron excitation
to convert to a hole excitation via second nearest-neighbor imaginary
intralayer hoppings, see Fig. 2(b).

Figure 2: (a) Low-energy bandstructure of TBG at $\theta\approx 1.08^{\circ}$
and filling $n=-2.4$, with and without scaling. (b) Hoppings inducing $f$-wave
superconducting pairing. (c) Central part of the lattice of the TBG Josephson
junction Mis ; Gonçalves  _et al._ (2021). The electrodes are superconductors
with a phase difference of $\phi$ and the link region, with a length of four
moiré periods, is metallic or insulating. The rhombus is a unit cell of TBG.
Figure 3: TBG Josephson junctions with equal electrodes. (a) Current-phase
relations for near magic-angle junctions with different pairing symmetry:
spin-singlet $s$-wave (dashed) or spin-triplet $f$-wave (solid), for electron
and hole superconducting domes (fillings $n=\pm 2.4$), with a metallic (SNS)
or insulating link (SIS) bMi  We set the superconducting gap to 1 meV Oh  _et
al._ (2021). $\theta=1.06^{\circ}$ for SNS; $1.1^{\circ}$ for SIS $h^{+}$ and
$1.16^{\circ}$ for SIS $e^{-}$. Units: nanoampere per nanometer junction
width. (b) Critical current versus twist angle for all configurations. (c)
Andreev spectrum at $1.06^{\circ}$. Inset: charge map of an Andreev bound
state. (d) Critical current in SIS JJs compared to the difference between the
superconducting and insulating gaps, as a function of twist angle, and a
sketch of the bands in the different regions of a SIS junction.

Josephson junctions. Results. Figure 3(a) compares the current-phase relations
(CPRs) of TBG JJs with $s$\- and $f$-wave pairings in multiple configurations.
CPRs can be measured with a SQUID geometry Portolés  _et al._ (2022). The main
message of Fig. 3 is that the type of pairing, $s$-wave or $f$-wave, plays a
minor role when both electrodes are equal, compare dashed and solid lines in
Fig. 3(a-b).

In SNS JJs the CPR is skewed, due to high transmission of Andreev bound
states, which carry over $80\%$ of the current in these junctions and are
mostly localized in $AA$ stacking regions, see Fig. 3(c). In contrast, in SIS
junctions the current comes from tunnelling states, so the CPR is sinusoidal
Ambegaokar and Baratoff (1963). An exception occurs when the insulating gap in
the link is comparable to the superconducting gap, resulting in skewness and
large currents. The current in SIS junctions exponentially depends on the
similitude between both gaps, see Fig. 3(d). We note that the authors of Ref.
Portolés  _et al._ (2022) report a sinusoidal CPR in TBG, without skewness,
despite having a SNS JJ. This may be due to low transmission in the junction
Golubov  _et al._ (2004). Figure 3(b) shows the critical current for all JJs
as a function of twist angle. For a comparison to experiments, see Ref. SM .
The current in SNS JJs increases with twist angle, suggesting that larger
Fermi velocities compensate the reduced density of states. Electron-hole
asymmetry is very notable, e.g. near $\theta=1.1^{\circ}$, the current in SIS
junctions with fillings $-2.4/4/-2.4$ is over two orders of magnitude larger
than with $2.4/-4/2.4$ due to the asymmetry in the size of the gaps between
narrow bands and electron- or hole-like remote bands.

Ref. de Vries  _et al._ (2021) reports a significant length dependence of the
critical current in JJs prepared in mixed configurations, e.g. with the
electrodes doped near one superconducting dome and the link near the other.
This indicates that the superconducting pairing symmetry in the electron and
hole domes may differ. The results in Fig. 4 for mixed $f$-wave and $s$-wave
TBG JJs propose an experiment that could verify the hypothesis. For these JJs,
the critical current dramatically depends on the angle between the junction
and the lattice. A similar result in non-superconducting junctions was found
in Ref. Alvarado and Yeyati (2021). The critical current is sizeable when the
junction axis is nearly parallel to the graphene armchair direction, but close
to zero when parallel to the zigzag direction. As long as the perpendicular
momentum is conserved, the zigzag JJ suffers destructive interference of the
superconducting pockets along the green lines drawn in Fig. 4. Also, the CPRs
have a period of $\pi$, half the one of standard JJs. The origin of this
effect is the existence of two sets of energy levels, due to coupling of the
$s$-wave pocket to the two $f$-wave pockets, which have an intrinsic phase
difference of $\pi$ Sigrist and Ueda (1991); Zagoskin (1997). Furthermore, the
CPR shows a $\pi$-junction behaviour, i.e. it is first negative Bulaevskii
_et al._ (1977); Golubov  _et al._ (2004). A requisite for these phenomena is
that the triplet electrode is spin unpolarized, otherwise the current is zero
due to spin conservation. The same occurs in a one-dimensional toy model SM ;
Zazunov and Egger (2012).

Figure 4: Current-phase relation in mixed $f$-wave and $s$-wave TBG Josephson
junctions, nearly parallel to the graphene armchair, as in Fig. 2(c), or
zigzag directions. The critical current is $\sim 100$ times larger for
armchair junctions.

Discussion. We have studied the role of the superconducting order parameter in
transport through superconducting TBG junctions. We focus on $s$\- and $f$\-
wave pairing (even and odd valley combinations), as these two choices are
equally favored by long range interactions, either attractive or repulsive
Scheurer and Samajdar (2020).

We have calculated the critical current, and the current-phase relation for
different types of Josephson junctions. JJs in which both electrodes are
either $s$\- or $f$-wave superconductors show similar features (unlike the
$s$\- and $p$\- cases considered in Ref. Linder  _et al._ (2009)). On the
other hand, the critical current in mixed ($s$\- and $f$-) junctions depends
strongly on the orientation of the junction with respect to the graphene
lattice axes, with maxima for armchair junctions, and zeroes for zigzag
junctions. Hence, mixed junctions are useful for determining whether two
superconducting phases differ in their valley exchange parity. Such junctions
can exist in various setups: i) different superconducting regions in the phase
diagram of TBG show different order parameters de Vries  _et al._ (2021), ii)
the superconducting state changes locally because of the spin-orbit coupling
induced by a substrate Zhang  _et al._ (2023), iii) superconducting TBG is
combined with $s$-wave proximitized graphene Heersche  _et al._ (2007); Rui
_et al._ (2020).

For a junction between a normal STM tip and superconducting TBG, we find a
prominent peak in subgap conductance for an $f$-wave order parameter, due to
Andreev states induced by the tip, in agreement with the experiments of Refs.
Oh  _et al._ (2021); Kim  _et al._ (2022). $f$-wave is also consistent with
the U- and V-shaped densities of states measured in the weak coupling regime,
as shown in Refs. Poduval and Scheurer (2023); SM . The agreement between the
experiments Oh  _et al._ (2021); Kim  _et al._ (2022) and the results
presented here puts forward $f$-wave pairing as a leading candidate for the
pairing symmetry of twisted graphene superconductors.

We note that the results for the STM-superconductor junction apply equally
well to all graphene superconductors Park  _et al._ (2021); Hao  _et al._
(2021); Park  _et al._ (2022); Zhang  _et al._ (2022); Zhou  _et al._ (2021,
2022); Zhang  _et al._ (2023); Holleis  _et al._ (2023). Extending the
Josephson junction calculations to non-twisted graphene superconductors Zhou
_et al._ (2021, 2022); Zhang  _et al._ (2023); Holleis  _et al._ (2023) is a
promising direction for future research.

Acknowledgements We thank Tommaso Cea, Shuichi Iwakiri, Klaus Ensslin and
Fernando de Juan for fruitful discussions. This work was supported by funding
from the European Commission, within the Graphene Flagship, Core 3, grant
number 881603; the Severo Ochoa programme for centres of excellence in R&D
(CEX2020-001039-S/AEI/10.13039/501100011033) and grants SprQuMat (Ministerio
de Ciencia e Innovación, Spain) and NMAT2D (Comunidad de Madrid, Spain).
V.T.P. acknowledges support from the Department of Energy under grant DE-
FG02-84ER45118, the NSF Graduate Research Fellowships Program, and the P.D.
Soros Fellowship for New Americans.

## References

  * Cao  _et al._ (2018a) Yuan Cao, Valla Fatemi, Ahmet Demir, Shiang Fang, Spencer L. Tomarken, Jason Y. Luo, Javier D. Sanchez-Yamagishi, Kenji Watanabe, Takashi Taniguchi, Efthimios Kaxiras, Ray C. Ashoori, and Pablo Jarillo-Herrero, ``Correlated insulator behaviour at half-filling in magic-angle graphene superlattices,'' Nature 556, 80–84 (2018a).
  * Polshyn  _et al._ (2019) Hryhoriy Polshyn, Matthew Yankowitz, Shaowen Chen, Yuxuan Zhang, K. Watanabe, T. Taniguchi, Cory R. Dean, and Andrea F. Young, ``Large linear-in-temperature resistivity in twisted bilayer graphene,'' Nature Physics 15, 1011–1016 (2019).
  * Sharpe  _et al._ (2019) Aaron L. Sharpe, Eli J. Fox, Arthur W. Barnard, Joe Finney, Kenji Watanabe, Takashi Taniguchi, M. A. Kastner, and David Goldhaber-Gordon, ``Emergent ferromagnetism near three-quarters filling in twisted bilayer graphene,'' Science 365, 605–608 (2019).
  * Serlin  _et al._ (2020) M. Serlin, C. L. Tschirhart, H. Polshyn, Y. Zhang, J. Zhu, K. Watanabe, T. Taniguchi, L. Balents, and A. F. Young, ``Intrinsic quantized anomalous hall effect in a moiré heterostructure,'' Science 367, 900–903 (2020).
  * Chen  _et al._ (2020) Guorui Chen, Aaron L. Sharpe, Eli J. Fox, Ya-Hui Zhang, Shaoxin Wang, Lili Jiang, Bosai Lyu, Hongyuan Li, Kenji Watanabe, Takashi Taniguchi, Zhiwen Shi, T. Senthil, David Goldhaber-Gordon, Yuanbo Zhang, and Feng Wang, ``Tunable correlated chern insulator and ferromagnetism in a moiré superlattice,'' Nature 579, 56–61 (2020).
  * 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–52 (2020).
  * Saito  _et al._ (2020) Yu Saito, Jingyuan Ge, Kenji Watanabe, Takashi Taniguchi, and Andrea F. Young, ``Independent superconductors and correlated insulators in twisted bilayer graphene,'' Nature Physics 16, 926–930 (2020).
  * Wong  _et al._ (2020) Dillon Wong, Kevin P. Nuckolls, Myungchul Oh, Biao Lian, Yonglong Xie, Sangjun Jeon, Kenji Watanabe, Takashi Taniguchi, B. Andrei Bernevig, and Ali Yazdani, ``Cascade of electronic transitions in magic-angle twisted bilayer graphene,'' Nature 582, 198–202 (2020).
  * Zondiner  _et al._ (2020) U. Zondiner, A. Rozen, D. Rodan-Legrain, Y. Cao, R. Queiroz, T. Taniguchi, K. Watanabe, Y. Oreg, F. von Oppen, Ady Stern, E. Berg, P. Jarillo-Herrero, and S. Ilani, ``Cascade of phase transitions and dirac revivals in magic-angle graphene,'' Nature 582, 203–208 (2020).
  * Stepanov  _et al._ (2020) Petr Stepanov, Ipsita Das, Xiaobo Lu, Ali Fahimniya, Kenji Watanabe, Takashi Taniguchi, Frank H. L. Koppens, Johannes Lischner, Leonid Levitov, and Dmitri K. Efetov, ``Untying the insulating and superconducting orders in magic-angle graphene,'' Nature 583, 375–378 (2020).
  * Xu  _et al._ (2020) Yang Xu, Song Liu, Daniel A. Rhodes, Kenji Watanabe, Takashi Taniguchi, James Hone, Veit Elser, Kin Fai Mak, and Jie Shan, ``Correlated insulating states at fractional fillings of moiré superlattices,'' Nature 587, 214–218 (2020).
  * Choi  _et al._ (2021) Youngjoon Choi, Hyunjin Kim, Yang Peng, Alex Thomson, Cyprian Lewandowski, Robert Polski, Yiran Zhang, Harpreet Singh Arora, Kenji Watanabe, Takashi Taniguchi, Jason Alicea, and Stevan Nadj-Perge, ``Correlation-driven topological phases in magic-angle twisted bilayer graphene,'' Nature 589, 536–541 (2021).
  * Rozen  _et al._ (2021) Asaf Rozen, Jeong Min Park, Uri Zondiner, Yuan Cao, Daniel Rodan-Legrain, Takashi Taniguchi, Kenji Watanabe, Yuval Oreg, Ady Stern, Erez Berg, Pablo Jarillo-Herrero, and Shahal Ilani, ``Entropic evidence for a pomeranchuk effect in magic-angle graphene,'' Nature 592, 214–219 (2021).
  * Cao  _et al._ (2021a) Yuan Cao, Daniel Rodan-Legrain, Jeong Min Park, Noah F. Q. Yuan, Kenji Watanabe, Takashi Taniguchi, Rafael M. Fernandes, Liang Fu, and Pablo Jarillo-Herrero, ``Nematicity and competing orders in superconducting magic-angle graphene,'' Science 372, 264–271 (2021a).
  * Stepanov  _et al._ (2021) Petr Stepanov, Ming Xie, Takashi Taniguchi, Kenji Watanabe, Xiaobo Lu, Allan H. MacDonald, B. Andrei Bernevig, and Dmitri K. Efetov, ``Competing zero-field chern insulators in superconducting twisted bilayer graphene,'' Phys. Rev. Lett. 127, 197701 (2021).
  * Oh  _et al._ (2021) Myungchul Oh, Kevin P. Nuckolls, Dillon Wong, Ryan L. Lee, Xiaomeng Liu, Kenji Watanabe, Takashi Taniguchi, and Ali Yazdani, ``Evidence for unconventional superconductivity in twisted bilayer graphene,'' Nature 600, 240–245 (2021).
  * Xie  _et al._ (2021) Yonglong Xie, Andrew T. Pierce, Jeong Min Park, Daniel E. Parker, Eslam Khalaf, Patrick Ledwith, Yuan Cao, Seung Hwan Lee, Shaowen Chen, Patrick R. Forrester, Kenji Watanabe, Takashi Taniguchi, Ashvin Vishwanath, Pablo Jarillo-Herrero, and Amir Yacoby, ``Fractional chern insulators in magic-angle twisted bilayer graphene,'' Nature 600, 439–443 (2021).
  * Berdyugin  _et al._ (2022) Alexey I. Berdyugin, Na Xin, Haoyang Gao, Sergey Slizovskiy, Zhiyu Dong, Shubhadeep Bhattacharjee, P. Kumaravadivel, Shuigang Xu, L. A. Ponomarenko, Matthew Holwill, D. A. Bandurin, Minsoo Kim, Yang Cao, M. T. Greenaway, K. S. Novoselov, I. V. Grigorieva, K. Watanabe, T. Taniguchi, V. I. Fal'ko, L. S. Levitov, Roshan Krishna Kumar, and A. K. Geim, ``Out-of-equilibrium criticalities in graphene superlattices,'' Science 375, 430–433 (2022).
  * Turkel  _et al._ (2022) Simon Turkel, Joshua Swann, Ziyan Zhu, Maine Christos, K. Watanabe, T. Taniguchi, Subir Sachdev, Mathias S. Scheurer, Efthimios Kaxiras, Cory R. Dean, and Abhay N. Pasupathy, ``Orderly disorder in magic-angle twisted trilayer graphene,'' Science 376, 193–199 (2022).
  * Huang  _et al._ (2022) Tianye Huang, Xuecou Tu, Changqing Shen, Binjie Zheng, Junzhuan Wang, Hao Wang, Kaveh Khaliji, Sang Hyun Park, Zhiyong Liu, Teng Yang, Zhidong Zhang, Lei Shao, Xuesong Li, Tony Low, Yi Shi, and Xiaomu Wang, ``Observation of chiral and slow plasmons in twisted bilayer graphene,'' Nature 605, 63–68 (2022).
  * de la Barrera  _et al._ (2022) Sergio C. de la Barrera, Samuel Aronson, Zhiren Zheng, Kenji Watanabe, Takashi Taniguchi, Qiong Ma, Pablo Jarillo-Herrero, and Raymond Ashoori, ``Cascade of isospin phase transitions in bernal-stacked bilayer graphene at zero magnetic field,'' Nature Physics 18, 771–775 (2022).
  * Kim  _et al._ (2022) Hyunjin Kim, Youngjoon Choi, Cyprian Lewandowski, Alex Thomson, Yiran Zhang, Robert Polski, Kenji Watanabe, Takashi Taniguchi, Jason Alicea, and Stevan Nadj-Perge, ``Evidence for unconventional superconductivity in twisted trilayer graphene,'' Nature 606, 494–500 (2022).
  * Seiler  _et al._ (2022) Anna M. Seiler, Fabian R. Geisenhof, Felix Winterer, Kenji Watanabe, Takashi Taniguchi, Tianyi Xu, Fan Zhang, and R. Thomas Weitz, ``Quantum cascade of correlated phases in trigonally warped bilayer graphene,'' Nature 608, 298–302 (2022).
  * Cao  _et al._ (2018b) Yuan Cao, Valla Fatemi, Shiang Fang, Kenji Watanabe, Takashi Taniguchi, Efthimios Kaxiras, and Pablo Jarillo-Herrero, ``Unconventional superconductivity in magic-angle graphene superlattices,'' Nature 556, 43–50 (2018b).
  * Yankowitz  _et al._ (2019) Matthew Yankowitz, Shaowen Chen, Hryhoriy Polshyn, Yuxuan Zhang, K. Watanabe, T. Taniguchi, David Graf, Andrea F. Young, and Cory R. Dean, ``Tuning superconductivity in twisted bilayer graphene,'' Science 363, 1059–1064 (2019).
  * Lu  _et al._ (2019) Xiaobo Lu, Petr Stepanov, Wei Yang, Ming Xie, Mohammed Ali Aamir, Ipsita Das, Carles Urgell, Kenji Watanabe, Takashi Taniguchi, Guangyu Zhang, Adrian Bachtold, Allan H. MacDonald, and Dmitri K. Efetov, ``Superconductors, orbital magnets and correlated states in magic-angle bilayer graphene,'' Nature 574, 653–657 (2019).
  * Park  _et al._ (2021) Jeong Min Park, Yuan Cao, Kenji Watanabe, Takashi Taniguchi, and Pablo Jarillo-Herrero, ``Tunable strongly coupled superconductivity in magic-angle twisted trilayer graphene,'' Nature 590, 249–255 (2021).
  * Hao  _et al._ (2021) Zeyu Hao, A. M. Zimmerman, Patrick Ledwith, Eslam Khalaf, Danial Haie Najafabadi, Kenji Watanabe, Takashi Taniguchi, Ashvin Vishwanath, and Philip Kim, ``Electric field–tunable superconductivity in alternating-twist magic-angle trilayer graphene,'' Science 371, 1133–1138 (2021).
  * Park  _et al._ (2022) Jeong Min Park, Yuan Cao, Li-Qiao Xia, Shuwen Sun, Kenji Watanabe, Takashi Taniguchi, and Pablo Jarillo-Herrero, ``Robust superconductivity in magic-angle multilayer graphene family,'' Nature Materials 21, 877–883 (2022).
  * Zhang  _et al._ (2022) Yiran Zhang, Robert Polski, Cyprian Lewandowski, Alex Thomson, Yang Peng, Youngjoon Choi, Hyunjin Kim, Kenji Watanabe, Takashi Taniguchi, Jason Alicea, Felix von Oppen, Gil Refael, and Stevan Nadj-Perge, ``Promotion of superconductivity in magic-angle graphene multilayers,'' Science 377, 1538–1543 (2022).
  * Zhou  _et al._ (2021) Haoxin Zhou, Tian Xie, Takashi Taniguchi, Kenji Watanabe, and Andrea F. Young, ``Superconductivity in rhombohedral trilayer graphene,'' Nature 598, 434–438 (2021).
  * Zhou  _et al._ (2022) Haoxin Zhou, Ludwig Holleis, Yu Saito, Liam Cohen, William Huynh, Caitlin L. Patterson, Fangyuan Yang, Takashi Taniguchi, Kenji Watanabe, and Andrea F. Young, ``Isospin magnetism and spin-polarized superconductivity in bernal bilayer graphene,'' Science 375, 774–778 (2022).
  * Zhang  _et al._ (2023) Yiran Zhang, Robert Polski, Alex Thomson, Étienne Lantagne-Hurtubise, Cyprian Lewandowski, Haoxin Zhou, Kenji Watanabe, Takashi Taniguchi, Jason Alicea, and Stevan Nadj-Perge, ``Enhanced superconductivity in spin–orbit proximitized bilayer graphene,'' Nature 613, 268–273 (2023).
  * Holleis  _et al._ (2023) Ludwig Holleis, Caitlin L. Patterson, Yiran Zhang, Heun Mo Yoo, Haoxin Zhou, Takashi Taniguchi, Kenji Watanabe, Stevan Nadj-Perge, and Andrea F. Young, ``Ising superconductivity and nematicity in bernal bilayer graphene with strong spin orbit coupling,''  (2023), 10.48550/arXiv.2303.00742.
  * Cao  _et al._ (2021b) Yuan Cao, Jeong Min Park, Kenji Watanabe, Takashi Taniguchi, and Pablo Jarillo-Herrero, ``Pauli-limit violation and re-entrant superconductivity in moiré graphene,'' Nature 595, 526–531 (2021b).
  * Lin  _et al._ (2022) Jiang-Xiazi Lin, Phum Siriviboon, Harley D. Scammell, Song Liu, Daniel Rhodes, K. Watanabe, T. Taniguchi, James Hone, Mathias S. Scheurer, and J.I.A. Li, ``Zero-field superconducting diode effect in small-twist-angle trilayer graphene,'' Nature Physics 18, 1221–1227 (2022).
  * Lake  _et al._ (2022) Ethan Lake, Adarsh S. Patri, and T. Senthil, ``Pairing symmetry of twisted bilayer graphene: A phenomenological synthesis,'' Phys. Rev. B 106, 104506 (2022).
  * de Vries  _et al._ (2021) Folkert K. de Vries, Elías Portolés, Giulia Zheng, Takashi Taniguchi, Kenji Watanabe, Thomas Ihn, Klaus Ensslin, and Peter Rickhaus, ``Gate-defined josephson junctions in magic-angle twisted bilayer graphene,'' Nature Nanotechnology 16, 760–763 (2021).
  * Rodan-Legrain  _et al._ (2021) Daniel Rodan-Legrain, Yuan Cao, Jeong Min Park, Sergio C. de la Barrera, Mallika T. Randeria, Kenji Watanabe, Takashi Taniguchi, and Pablo Jarillo-Herrero, ``Highly tunable junctions and non-local josephson effect in magic-angle graphene tunnelling devices,'' Nature Nanotechnology 16, 769–775 (2021).
  * Diez-Merida  _et al._ (2021) J. Diez-Merida, A. Diez-Carlon, S. Y. Yang, Y. M. Xie, X. J. Gao, K. Watanabe, T. Taniguchi, X. Lu, K. T. Law, and Dmitri K. Efetov, ``Magnetic josephson junctions and superconducting diodes in magic angle twisted bilayer graphene,''  (2021), 10.48550/arXiv.2110.01067.
  * Xie  _et al._ (2023) Ying-Ming Xie, Dmitri K. Efetov, and K. T. Law, ``$\varphi_{0}$-Josephson junction in twisted bilayer graphene induced by a valley-polarized state,'' Physical Review Research 5 (2023), 10.1103/physrevresearch.5.023029.
  * Hu  _et al._ (2022) Jin-Xin Hu, Zi-Ting Sun, Ying-Ming Xie, and K. T. Law, ``Valley polarization induced josephson diode effect in twisted bilayer graphene,''  (2022), 10.48550/arXiv.2211.14846.
  * Portolés  _et al._ (2022) Elías Portolés, Shuichi Iwakiri, Giulia Zheng, Peter Rickhaus, Takashi Taniguchi, Kenji Watanabe, Thomas Ihn, Klaus Ensslin, and Folkert K. de Vries, ``A tunable monolithic SQUID in twisted bilayer graphene,'' Nature Nanotechnology 17, 1159–1164 (2022).
  * Sukhachov  _et al._ (2023) P. O. Sukhachov, Felix von Oppen, and L. I. Glazman, ``Andreev reflection in scanning tunneling spectroscopy of unconventional superconductors,'' Physical Review Letters 130 (2023), 10.1103/physrevlett.130.216002.
  * Josephson (1962) B.D. Josephson, ``Possible new effects in superconductive tunnelling,'' Physics Letters 1, 251–253 (1962).
  * Sigrist and Ueda (1991) Manfred Sigrist and Kazuo Ueda, ``Phenomenological theory of unconventional superconductivity,'' Reviews of Modern Physics 63, 239–311 (1991).
  * Golubov  _et al._ (2004) A. A. Golubov, M. Yu. Kupriyanov, and E. Il'ichev, ``The current-phase relation in josephson junctions,'' Reviews of Modern Physics 76, 411–469 (2004).
  * Tsuei  _et al._ (1994) C. C. Tsuei, J. R. Kirtley, C. C. Chi, Lock See Yu-Jahnes, A. Gupta, T. Shaw, J. Z. Sun, and M. B. Ketchen, ``Pairing symmetry and flux quantization in a tricrystal superconducting ring of $\mathrm{Y}{\mathrm{ba}}_{2}{\mathrm{cu}}_{3}{\mathrm{o}}_{7-\delta}$,'' Phys. Rev. Lett. 73, 593–596 (1994).
  * Tsuei and Kirtley (2000) C. C. Tsuei and J. R. Kirtley, ``Pairing symmetry in cuprate superconductors,'' Rev. Mod. Phys. 72, 969–1016 (2000).
  * Crépel  _et al._ (2022) Valentin Crépel, Tommaso Cea, Liang Fu, and Francisco Guinea, ``Unconventional superconductivity due to interband polarization,'' Phys. Rev. B 105, 094506 (2022).
  * Cea and Guinea (2021) Tommaso Cea and Francisco Guinea, ``Coulomb interaction, phonons, and superconductivity in twisted bilayer graphene,'' Proceedings of the National Academy of Sciences 118 (2021), 10.1073/pnas.2107874118.
  * Blonder  _et al._ (1982) G. E. Blonder, M. Tinkham, and T. M. Klapwijk, ``Transition from metallic to tunneling regimes in superconducting microconstrictions: Excess current, charge imbalance, and supercurrent conversion,'' Physical Review B 25, 4515–4532 (1982).
  * Setiawan and Hofmann (2022) F. Setiawan and Johannes Hofmann, ``Analytic approach to transport in superconducting junctions with arbitrary carrier density,'' Physical Review Research 4 (2022), 10.1103/physrevresearch.4.043087.
  * Lewandowski  _et al._ (2023) Cyprian Lewandowski, Étienne Lantagne-Hurtubise, Alex Thomson, Stevan Nadj-Perge, and Jason Alicea, ``Andreev reflection spectroscopy in strongly paired superconductors,'' Physical Review B 107 (2023), 10.1103/physrevb.107.l020502.
  * (55) In the case of the $f$-wave superconductor, there is no restriction on the value of a given component of the spin, which can be $s_{z}=0,\pm 1$.
  * (56) See supplementary information, which includes Refs. Moon  _et al._ (2014); Bistritzer and MacDonald (2011); Perfetto  _et al._ (2009); Kitaev (2001); Ishii (1970).
  * Anderson (1959) P.W. Anderson, ``Theory of dirty superconductors,'' Journal of Physics and Chemistry of Solids 11, 26–30 (1959).
  * Lin and Tománek (2018) Xianqing Lin and David Tománek, ``Minimum model for the electronic structure of twisted bilayer graphene and related structures,'' Phys. Rev. B 98, 081410 (2018).
  * Guinea and Walet (2018) Francisco Guinea and Niels R. Walet, ``Electrostatic effects, band distortions, and superconductivity in twisted graphene bilayers,'' Proceedings of the National Academy of Sciences 115, 13174–13179 (2018).
  * Rademaker  _et al._ (2019) Louk Rademaker, Dmitry A. Abanin, and Paula Mellado, ``Charge smoothening and band flattening due to hartree corrections in twisted bilayer graphene,'' Phys. Rev. B 100, 205114 (2019).
  * Gonzalez-Arraga  _et al._ (2017) Luis A. Gonzalez-Arraga, J. L. Lado, Francisco Guinea, and Pablo San-Jose, ``Electrically controllable magnetism in twisted bilayer graphene,'' Phys. Rev. Lett. 119, 107201 (2017).
  * Vahedi  _et al._ (2021) Javad Vahedi, Robert Peters, Ahmed Missaoui, Andreas Honecker, and Guy Trambly de Laissardière, ``Magnetism of magic-angle twisted bilayer graphene,'' SciPost Physics 11 (2021), 10.21468/scipostphys.11.4.083.
  * Sainz-Cruz  _et al._ (2021) Héctor Sainz-Cruz, Tommaso Cea, Pierre A. Pantaleón, and Francisco Guinea, ``High transmission in twisted bilayer graphene with angle disorder,'' Phys. Rev. B 104, 075144 (2021).
  * Haldane (1988) F. D. M. Haldane, ``Model for a quantum hall effect without landau levels: Condensed-matter realization of the "parity anomaly",'' Physical Review Letters 61, 2015–2018 (1988).
  * (65) To induce a superconducting gap of magnitude $\Delta$, the necessary amplitudes are $\Delta_{S}=\Delta/2$ for $s$-wave, $\Delta_{K}=\Delta/4$ for Kitaev and $\Delta_{F}\approx\Delta/(6\sqrt{3})$ for $f$-wave pairing.
  * (66) We use a commensurate lattice, but we caution that incommensurability has an impact on transport Gonçalves  _et al._ (2021).
  * Gonçalves  _et al._ (2021) Miguel Gonçalves, Hadi Z Olyaei, Bruno Amorim, Rubem Mondaini, Pedro Ribeiro, and Eduardo V Castro, ``Incommensurability-induced sub-ballistic narrow-band-states in twisted bilayer graphene,'' 2D Materials 9, 011001 (2021).
  * (68) In SIS junctions with the leads at $n=+2.4$ ($-2.4$), the chemical potential of the link is placed in the middle of the gap between the flat bands and the hole-like (electron-like) remote bands.
  * Ambegaokar and Baratoff (1963) Vinay Ambegaokar and Alexis Baratoff, ``Tunneling between superconductors,'' Physical Review Letters 10, 486–489 (1963).
  * Alvarado and Yeyati (2021) M. Alvarado and A. Levy Yeyati, ``Transport and spectral properties of magic-angle twisted bilayer graphene junctions based on local orbital models,'' Phys. Rev. B 104, 075406 (2021).
  * Zagoskin (1997) Alexandre M Zagoskin, ``The half-periodic josephson effect in an s-wave superconductor - normal-metal - d-wave superconductor junction,'' Journal of Physics: Condensed Matter 9, L419–L426 (1997).
  * Bulaevskii  _et al._ (1977) L N Bulaevskii, V V Kuzii, and A A Sobyanin, ``Superconducting system with weak coupling to the current in the ground state,'' JETP Lett. 25 (1977).
  * Zazunov and Egger (2012) Alex Zazunov and Reinhold Egger, ``Supercurrent blockade in josephson junctions with a majorana wire,'' Phys. Rev. B 85, 104514 (2012).
  * Scheurer and Samajdar (2020) Mathias S. Scheurer and Rhine Samajdar, ``Pairing in graphene-based moiré superlattices,'' Phys. Rev. Research 2, 033062 (2020).
  * Linder  _et al._ (2009) Jacob Linder, Annica M. Black-Schaffer, Takehito Yokoyama, Sebastian Doniach, and Asle Sudbø, ``Josephson current in graphene: Role of unconventional pairing symmetries,'' Phys. Rev. B 80, 094522 (2009).
  * Heersche  _et al._ (2007) Hubert B. Heersche, Pablo Jarillo-Herrero, Jeroen B. Oostinga, Lieven M. K. Vandersypen, and Alberto F. Morpurgo, ``Bipolar supercurrent in graphene,'' Nature 446, 56–59 (2007).
  * Rui  _et al._ (2020) Dingran Rui, Luzhao Sun, N. Kang, Hailin Peng, Zhongfan Liu, and H. Q. Xu, ``Superconductivity in an al-twisted bilayer graphene-al junction device,'' Japanese Journal of Applied Physics 59, SGGI07 (2020).
  * Poduval and Scheurer (2023) Prathyush P. Poduval and Mathias S. Scheurer, ``Vestigial singlet pairing in a fluctuating magnetic triplet superconductor: Applications to graphene moiré systems,''  (2023), 10.48550/arXiv.2301.01344.
  * Moon  _et al._ (2014) Pilkyung Moon, Young-Woo Son, and Mikito Koshino, ``Optical absorption of twisted bilayer graphene with interlayer potential asymmetry,'' Phys. Rev. B 90, 155427 (2014).
  * Bistritzer and MacDonald (2011) Rafi Bistritzer and Allan H. MacDonald, ``Moiré bands in twisted double-layer graphene,'' Proceedings of the National Academy of Sciences 108, 12233–12237 (2011).
  * Perfetto  _et al._ (2009) Enrico Perfetto, Gianluca Stefanucci, and Michele Cini, ``Equilibrium and time-dependent josephson current in one-dimensional superconducting junctions,'' Phys. Rev. B 80, 205408 (2009).
  * Kitaev (2001) A Yu Kitaev, ``Unpaired majorana fermions in quantum wires,'' Physics-Uspekhi 44, 131–136 (2001).
  * Ishii (1970) Chikara Ishii, ``Josephson currents through junctions with normal metal barriers,'' Progress of Theoretical Physics 44, 1525–1547 (1970).

## Supplementary information

## S1. STM tip-superconductor junctions

###  .1 Equations

There are six equations relating the scattering amplitudes, as each node leads
to an equation for the electron part of the wavefunction, and to another
equation for the hole part, see Fig. 1(a). These equations are:

$\displaystyle\epsilon(1+R_{N})$
$\displaystyle=t_{tip}(e_{N}^{*}+R_{N}e_{N})+t_{tip,K}\left(\frac{T_{K,1}}{N_{K,1}}+\frac{T_{K,2}}{N_{K,2}}\right)+t_{tip,K^{\prime}}\left(\frac{T_{K^{\prime},1}}{N_{K^{\prime},1}}+\frac{T_{K^{\prime},2}}{N_{K^{\prime},2}}\right)$
$\displaystyle\epsilon\left(\frac{T_{K,1}}{N_{K,1}}+\frac{T_{K,2}}{N_{K,2}}\right)$
$\displaystyle=t^{*}_{tip,K}(1+R_{N})+t_{sc,K}\left(\frac{T_{K,1}}{N_{K,1}}e_{K,1}+\frac{T_{K,2}}{N_{K,2}}e_{K,2}\right)+t_{K,K^{\prime}}\left(\frac{T_{K^{\prime},1}}{N_{K^{\prime},1}}+\frac{T_{K^{\prime},2}}{N_{K^{\prime},2}}\right)$
$\displaystyle\epsilon\left(\frac{T_{K^{\prime},1}}{N_{K^{\prime},1}}+\frac{T_{K^{\prime},2}}{N_{K^{\prime},2}}\right)$
$\displaystyle=t^{*}_{tip,K^{\prime}}(1+R_{N})+t_{sc,K^{\prime}}\left(\frac{T_{K^{\prime},1}}{N_{K^{\prime},1}}e_{K^{\prime},1}+\frac{T_{K^{\prime},2}}{N_{K^{\prime},2}}e_{K^{\prime},2}\right)+t_{K,K^{\prime}}^{*}\left(\frac{T_{K,1}}{N_{K,1}}+\frac{T_{K,2}}{N_{K,2}}\right)$
$\displaystyle\epsilon R_{A}$
$\displaystyle=-t_{tip}R_{A}e_{N}^{*}-t_{tip,K}\left(\frac{T_{K,1}A_{K,1}}{N_{K,1}}+\frac{T_{K,2}A_{K,2}}{N_{K,2}}\right)-t_{tip,K^{\prime}}\left(\frac{T_{K^{\prime},1}A_{K^{\prime},1}}{N_{K^{\prime},1}}+\frac{T_{K^{\prime},2}A_{K^{\prime},2}}{N_{K^{\prime},2}}\right)$
$\displaystyle\epsilon\left(\frac{T_{K,1}A_{K,1}}{N_{K,1}}+\frac{T_{K,2}A_{K,2}}{N_{K,2}}\right)$
$\displaystyle=-t_{tip,K}^{*}R_{A}-t_{sc,K}\left(\frac{T_{K,1}A_{K,1}}{N_{K,1}}e_{K,1}+\frac{T_{K,2}A_{K,2}}{N_{K,2}}e_{K,2}\right)$
$\displaystyle-
t_{K,K^{\prime}}\left(\frac{T_{K^{\prime},1}A_{K^{\prime},1}}{N_{K^{\prime},1}}+\frac{T_{K^{\prime},2}A_{K^{\prime},2}}{N_{K^{\prime},2}}\right)$
$\displaystyle\epsilon\left(\frac{T_{K^{\prime},1}A_{K^{\prime},1}}{N_{K^{\prime},1}}+\frac{T_{K^{\prime},2}A_{K^{\prime},2}}{N_{K^{\prime},2}}\right)$
$\displaystyle=-t_{tip,K^{\prime}}^{*}R_{A}-t_{sc,K^{\prime}}\left(\frac{T_{K^{\prime},1}A_{K^{\prime},1}}{N_{K^{\prime},1}}e_{K^{\prime},1}+\frac{T_{K^{\prime},2}A_{K^{\prime},2}}{N_{K^{\prime},2}}e_{K^{\prime},2}\right)$
$\displaystyle-
t_{K,K^{\prime}}^{*}\left(\frac{T_{K,1}A_{K,1}}{N_{K,1}}+\frac{T_{K,2}A_{K,2}}{N_{K,2}}\right)$
(S1)

where $R_{N},R_{A},T_{K,1},T_{K,2},T_{K^{\prime}_{1}},T_{K^{\prime},2}$ are
the normal reflection coefficient, the Andreev reflection coefficient
(associated to the tip electron and hole channels, respectively), and the
transmission coefficients for the electron and hole channels in each valley,
$\epsilon$ is the energy, and:

$\displaystyle e_{N}$
$\displaystyle=-\frac{\epsilon}{2t_{tip}}+i\sqrt{1-\frac{\epsilon^{2}}{4t_{tip}^{2}}}$
$\displaystyle e_{K,1}$
$\displaystyle=\left\\{\begin{array}[]{cc}-i\sqrt{1-\frac{\Delta_{K}^{2}-\epsilon^{2}}{4t_{sc,K}^{2}}}+i\sqrt{\frac{\Delta_{K}^{2}-\epsilon^{2}}{4t_{sc,K}^{2}}}&|\epsilon|\leq\Delta_{K}\\\
-i\sqrt{1+\frac{\epsilon^{2}-\Delta_{K}^{2}}{4t_{sc,K}^{2}}}+\sqrt{\frac{\epsilon^{2}-\Delta_{K}^{2}}{4t_{sc,K}^{2}}}&|\epsilon|>\Delta_{K}\end{array}\right.$
$\displaystyle e_{K,2}$
$\displaystyle=\left\\{\begin{array}[]{cc}i\sqrt{1-\frac{\Delta_{K}^{2}-\epsilon^{2}}{4t_{sc,K}^{2}}}-i\sqrt{\frac{\Delta_{K}^{2}-\epsilon^{2}}{4t_{sc,K}^{2}}}&|\epsilon|\leq\Delta_{K}\\\
i\sqrt{1-\frac{\epsilon^{2}-\Delta_{K}^{2}}{4t_{sc,K}^{2}}}+\sqrt{\frac{\epsilon^{2}-\Delta_{K}^{2}}{4t_{sc,K}^{2}}}&|\epsilon|>\Delta_{K}\end{array}\right.$
(S6) $\displaystyle A_{K,1}$
$\displaystyle=\left\\{\begin{array}[]{cc}\frac{-\epsilon-i\sqrt{\Delta_{K}^{2}-\epsilon^{2}}}{\Delta_{K}}&|\epsilon|\leq\Delta_{K}\\\
\frac{-\epsilon-\sqrt{\epsilon^{2}-\Delta_{K}^{2}}}{\Delta_{K}}&|\epsilon|>\Delta_{K}\end{array}\right.$
$\displaystyle A_{K,2}$
$\displaystyle=\left\\{\begin{array}[]{cc}\frac{-\epsilon+i\sqrt{\Delta_{K}^{2}-\epsilon^{2}}}{\Delta_{K}}&|\epsilon|\leq\Delta_{K}\\\
\frac{-\epsilon+\sqrt{\epsilon^{2}-\Delta_{K}^{2}}}{\Delta_{K}}&|\epsilon|>\Delta_{K}\end{array}\right.$
(S11) $\displaystyle N_{K,1}$ $\displaystyle=\sqrt{1+|A_{K,1}|^{2}}$
$\displaystyle N_{K,2}$ $\displaystyle=\sqrt{1+|A_{K,2}|^{2}}$ (S12)

and similar expressions for $K^{\prime}$.

The model can be extended to intravalley superconducting gaps with angular
dependence, like the cases discussed in Ref. Sukhachov  _et al._ (2023). Each
individual scattering angle can be treated as an independent channel.
Reflection and transmission coefficients need to be defined for each angle,
and the total currents will be given by integrals over all angles.

###  .2 Subgap Andreev conductance.

The enhancement of the conductance for voltages below the superconducting gap
discussed in Ref. Blonder  _et al._ (1982) arises from processes where an
incoming electron is reflected as a hole, or vice versa. In an $s$-wave
superconductor, and in the limit of perfect transmission, these processes lead
to a conductance which is twice the conductance in the normal state, as shown
in the BTK theory Blonder  _et al._ (1982). In an $s$-wave superconductor, an
electron is injected into the superconductor as a coherent sum of even
combinations of plane waves with momenta $\vec{k}$ and $-\vec{k}$, because due
to time reversal symmetry, the normal-superconductor hopping elements satisfy
$t_{tip,\vec{k}}=t_{tip,-\vec{k}}$. This even combination becomes coupled to
another even combination of hole states, which can move back into the normal
electrode. In a non $s$-wave superconductor, the superconducting gap changes
sign. If there are pairs of momenta such that
$\Delta_{\vec{k}}=-\Delta_{\vec{k^{\prime}}}$ and
$t_{tip,\vec{k}}=t_{tip,\vec{k^{\prime}}}$, the amplitudes of the injected
electron will be equal for $\vec{k}$ and $\vec{k^{\prime}}$, but, inside the
superconductor, it will be coupled to holes with amplitudes of opposite signs.
Such a hole state cannot tunnel back into the normal electrode, and subgap
Andreev conductance will be fully suppressed. Local tunneling processes, as
expected in an STM experiment, imply momentum independent hopping elements,
$t_{tip,\vec{k}}=t$. Hence, we can expect that the subgap Andreev conductance
will be suppressed when a normal tip is coupled to generic $p$\- and $d$-wave
superconductors Sukhachov  _et al._ (2023), and also in the case of the
$f$-wave, two valley superconductor considered here. The situation changes
when there is intervalley scattering, see below.

###  .3 Tip induced Andreev states.

We can understand the formation of subgap Andreev states by the intervalley
elastic scattering induced by the tip by using the model shown in Fig. 1(a).
The model reduces to a simple tight binding model:

$\displaystyle{\cal H}$ $\displaystyle={\cal H}_{sc1}+{\cal H}_{sc2}+{\cal
H}_{tip}$ $\displaystyle{\cal H}_{sc1}$
$\displaystyle=t_{sc1}\sum_{n=-\infty}^{0}\left(c_{e,n}^{\dagger}c_{e,n-1}-c_{h,n}^{\dagger}c_{h,n-1}\right)+\Delta_{sc1}\sum_{n=-\infty}^{0}c_{e,n}^{\dagger}c_{h,n}+h.c.$
$\displaystyle{\cal H}_{sc2}$
$\displaystyle=t_{sc2}\sum_{n=1}^{n=\infty}\left(c_{e,n}^{\dagger}c_{e,n+1}-c_{h,n}^{\dagger}c_{h,n+1}\right)+\Delta_{sc2}\sum_{n=1}^{n=\infty}c_{e,n}^{\dagger}c_{h,n}+h.c.$
$\displaystyle{\cal H}_{tip}$
$\displaystyle=t_{K,K^{\prime}}\left(c_{e,0}^{\dagger}c_{e,1}-c_{h,0}^{\dagger}c_{h,1}\right)+h.c.$
(S13)

The Green's function at sites $n=0$ and $n=1$ of the system can be obtained
from the Green's functions at the same sites in the absence of intervalley
coupling:

$\displaystyle\left(\begin{array}[]{cc}G_{0,0}(\omega)&G_{0,1}(\omega)\\\
G_{1,0}(\omega)&G_{1,1}(\omega)\end{array}\right)$
$\displaystyle=\left(\begin{array}[]{cc}\bar{G}_{0,0}^{-1}(\omega)&t_{K,K^{\prime}}{\cal
I}_{2}\\\ t_{K,K^{\prime}}{\cal
I}_{2}&\bar{G}_{1,1}^{-1}(\omega)\end{array}\right)^{-1}$ (S18)

where $\bar{G}_{0,0}(\omega)$ and $\bar{G}_{1,1}(\omega)$ are surface Green's
functions associated to ${\cal H}_{sc1}$ and ${\cal H}_{sc2}$, and ${\cal
I}_{2}$ is a $2\times 2$ identity matrix. By changing to a basis defined by
$c_{e,n}\pm c_{h,n}$, these matrix functions are:

$\displaystyle\bar{G}_{0,0}^{-1}(\omega)$
$\displaystyle=\left(\begin{array}[]{cc}\frac{\omega-\Delta_{sc1}}{2}+\frac{\sqrt{(\omega^{2}-\Delta_{sc1}^{2})(\omega^{2}-\Delta_{sc1}^{2}-4t_{sc1}^{2})}}{2(\omega+\Delta_{sc1})}&0\\\
0&\frac{\omega+\Delta_{sc1}}{2}+\frac{\sqrt{(\omega^{2}-\Delta_{sc1}^{2})(\omega^{2}-\Delta_{sc1}^{2}-4t_{sc1}^{2})}}{2(\omega-\Delta_{sc1})}\end{array}\right)$
(S21)

and an equivalent expression for $G_{1,1}^{-1}(\omega)$. Finally, for
$t_{sc1}=t_{sc2}=t$ and $\Delta_{sc1}=\Delta_{sc2}=\Delta$, the Andreev states
are defined by the equations:

$\displaystyle\frac{\omega\mp\Delta}{2}+\frac{\sqrt{(\omega^{2}-\Delta^{2})(\omega^{2}-\Delta^{2}-4t^{2})}}{2(\omega\pm\Delta)}$
$\displaystyle=\pm t_{K,K^{\prime}}$ (S22)

For $t_{K,K^{\prime}}\ll\Delta,t$ this equation gives Andreev states near the
edge of the superconducting gap, $\omega=\pm\Delta$, and for
$t_{K,K^{\prime}}=t$ the Andreev states move to the center of the gap,
$\omega=0$. The parameter $t$ describes a high energy cutoff of the order of
the bandwidth. For TBG near a magic angle, it is reasonable to expect that the
perturbation due to the tip in the contact regime is such that
$t_{K,K^{\prime}}\gtrsim t$, so that, in an $f$-wave superconductor, Andreev
states near the center of the gap will exist.

###  .4 Weak coupling conductance

Figure S1: Normal and Andreev reflection (top) and total conductance (bottom)
in a STM tip-superconducting TBG junction in the weak coupling regime. (a)
With Fermi velocity mismatch and intervalley scattering. The parameters used
are
$t_{tip}=10,\,t_{sc,K}=t_{sc,K^{\prime}}=1,\,t_{tip,K}=t_{tip,K^{\prime}}=1/\sqrt{2},\,t_{K,K^{\prime}}=0.2,\,\Delta_{K}=0.05,\,\Delta_{K^{\prime}}=\pm
0.05$. (b) With Fermi velocity mismatch, intervalley scattering and spin-orbit
coupling. $\Delta_{K}=0.05,\,\Delta_{K^{\prime}}=\pm 0.02$, others as in (a).
The insets show a zoom near the edge of the smaller gap.

We show in Fig. S1(a) results in the regime where the normal transmission of
the junction is small. The intervalley coupling has also been reduced. The
Andreev reflection for the $s$-wave superconductor is notably suppressed. On
the other hand, the intervalley coupling still induces subgap states near the
edges of the gap of the $f$-wave superconductor. As a result, Andreev
reflection persists, and leads to a peak in the junction conductance. Still,
the conductance above and below the gap are similar for $s$-wave and $f$-wave
superconductor. Therefore, in this weak coupling regime, it would be difficult
for an experiment to tell apart the two pairings. As shown in Fig. S1(b), the
transport characteristics of $f$-wave and $s$-wave superconductors remain
similar to each other when spin-orbit coupling is included.

## S2. Josephson junctions

###  .5 Critical current: comparison to experiments

In this section, we compare the critical currents obtained with our
microscopic simulations to the ones measured in the experiments of Ref. de
Vries  _et al._ (2021) In the main text we have used a superconducting gap of
1 meV, inferred from STM measurements Oh  _et al._ (2021). However the
transport data of Ref. de Vries  _et al._ (2021) suggests a gap of $\sim 0.1$
meV instead, so we set $\Delta=0.1$ meV in the following calculations. The
experimental twist angle is 1.06${}^{\circ}\pm 0.04^{\circ}$ and the junction
width is 1200 nm. Simulating a junction of this width is numerically
prohibitive, so we extrapolate the results for a junction of width $\sim 50$
nm. We focus on the setup with a link of length $L_{j}\approx 100$ nm. Note
that, since the gate-defined potential profile is expected to be smooth, the
effective link is most likely shorter than 100 nm.

In this setup, when the electrodes are tuned to the optimal filling for
superconductivity ($n\approx-2.4$) and the link is slightly detuned from that
($n\approx-2.4+\delta n$), the junction is in a SNS configuration and the
measured critical current is $\mathcal{I}_{c}\approx 250$ nA. When the link is
doped instead to the gap between the narrow bands and the hole-like remote
bands ($n=-4$), the junction is in a SIS configuration and the critical
current is $\mathcal{I}_{c}\approx 40$ nA.

Experiment de Vries  _et al._ (2021)
---
Twist angle |  Link length |  SNS $\mathcal{I}_{c}$ |  SIS $\mathcal{I}_{c}$
1.06${}^{\circ}\pm 0.04^{\circ}$ |  100 nm |  250 nA |  40 nA
Theory
Twist angle and pairing |  Link length |  SNS $\mathcal{I}_{c}$ |  SIS $\mathcal{I}_{c}$
1.06∘ and $s$-wave |  80 nm |  420 nA |  10-4 nA
1.06∘ and $f$-wave |  80 nm |  240 nA |  10-4 nA
1.06∘ and $s$-wave |  40 nm |  490 nA |  0.1 nA
1.12∘ and $s$-wave |  80 nm |  350 nA |  3 nA
1.12∘ and $s$-wave |  40 nm |  420 nA |  40 nA
1.12∘ and $f$-wave |  40 nm |  325 nA |  16 nA
Table 1: Critical current in TBG Josephson junctions: comparison between
experiment and theory.

In Table 1, we present the results of the calculations in the two
configurations. We consider two possibilities for the twist angle (1.06∘ and
1.12∘), as well as a longer (80 nm) and a shorter (40 nm) link. In some cases
we also compare $s$\- to $f$-wave pairing. We set the dielectric constant to
$\epsilon_{r}=7$. The results are in very good agreement with the experiments.
In particular, the critical current in the SNS configuration always matches
the order magnitude of the experimental one, regardless of details such as the
twist angle or the effective length of the link. In the case of the SIS
junction, the current depends exponentially on the length of the link and the
insulating gap between narrow and remote bands. Crucially, this insulating gap
is 8 meV for 1.06∘ and diminishes to 4 meV at 1.12∘. At 1.12∘ and with a 40 nm
link, the SIS current matches the experimental value. Therefore, a junction
with this angle and link length reproduces simultaneously the SNS and SIS
critical currents measured in Ref. de Vries  _et al._ (2021). Note that an
angle of 1.12∘ is well within the combined error bar of the experiment plus
the uncertainty in the tight binding parameters. The results suggest that the
top gate of 100 nm used in the experiment induces a rather smooth potential
profile, leading to an effective link around half that length. A complementary
or even alternative explanation is that the insulating phase has some metallic
character due to defects. $s$\- and $f$-wave lead to similar critical
currents.

###  .6 Details of the model

The junction lattice is built following the same procedure as in Ref. Sainz-
Cruz  _et al._ (2021). To avoid border transport, we impose periodic boundary
conditions from top to bottom, which leads to a folding of the Brillouin zone.
The folded bandstructure has more than two flat bands, e.g. four in Fig. 2(a).
Following the notation of Sainz-Cruz  _et al._ (2021), we build a TBG nanotube
with chiral vectors (44,2)@(-44,-2), which has a twist angle of 4.41∘ and 2704
sites in its unit cell.

We use a tight binding Hamiltonian given by Lin and Tománek (2018)

${{\cal H}_{0}=-\sum_{i\neq
j,m}\gamma_{ij}^{mm}(c^{\dagger}_{i,m}c_{j,m}+h.c.)}-\sum_{i,j,m}\gamma_{ij}^{m,m+1}(c^{\dagger}_{i,m}c_{j,m+1}+h.c.)+\sum_{i,m}V_{H}(n)c^{\dagger}_{i,m}c_{i,m}\,,$
(S23)

where $i,j$ run over the lattice sites and $m$ is the layer index. $H_{0}$
includes intralayer hopping to nearest-neighbors only
$\gamma_{ij}^{mm}=t_{\parallel}$ and interlayer hopping that decays
exponentially away from the vertical direction,
$\gamma_{ij}^{m,m+1}=t_{\perp}e^{-(\sqrt{r^{2}+d^{2}}-d)/\lambda_{\perp}}\frac{d^{2}}{r^{2}+d^{2}}$,
where $d=0.335$ nm is the distance between layers, $t_{\parallel}=3.09$ eV and
$t_{\perp}=0.39$ eV are the intralayer and interlayer hopping amplitudes and
$\lambda_{\perp}=0.027$ nm is a cutoff for the interlayer hopping Lin and
Tománek (2018).

We perform a scaling approximation, based on the fact that, within the
continuum model, the bands of TBG depend, to first order, on a dimensionless
parameter Bistritzer and MacDonald (2011),

$\alpha=\frac{at_{\perp}}{2\hbar
v_{F}\sin(\theta/2)}\propto\frac{t_{\perp}}{t_{\parallel}\theta}\,.$ (S24)

where $a$ is the lattice constant, $v_{F}$ is the Fermi velocity. Thus, a
small angle $\theta$ can be simulated with a larger one $\theta^{\prime}$,
doing the following transformations:
$t_{\parallel}\rightarrow\frac{1}{\lambda}t_{\parallel}$, $a\rightarrow\lambda
a$, $d\rightarrow\lambda d$, with
$\lambda=\sin(\frac{\theta^{\prime}}{2})/\sin(\frac{\theta}{2})$ Gonzalez-
Arraga  _et al._ (2017); Vahedi  _et al._ (2021); Sainz-Cruz  _et al._ (2021).
This approximation reproduces well the low-energy bandstructure, as shown in
Fig. 2(a) in the main text. It is worth noting that scaling leads to a rigid
blueshift of the bandstructure, of up to $\sim 20$ meV, which we have removed
in the figure. We use scaling factors $\lambda\sim$ 4.

In Eq. (S23), the Hartree term is
$V_{H}(n)=\frac{2\rho(n)}{\epsilon_{r}L_{M}}\sum^{3}_{i=1}\cos({G_{i}}\cdot{r})$,
where $G_{i}$ are the reciprocal lattice vectors, $r$ the position, $L_{M}$
the moiré period, $\epsilon_{r}$ the dielectric constant due to hBN
encapsulation and $\rho(n)$ a filling dependent parameter, listed in Table 2.
Realistic values for the dielectric constant of hBN-encapsulated twisted
bilayer graphene are usually in the range $\epsilon_{r}=4-10$, and our
original choice was $\epsilon_{r}=4$. However, this value leads to band
distortions so severe that, for some twist angles, the gaps between the
central narrow bands and the remote dispersive bands are closed. Therefore it
is not possible to form a proper SIS junction in these cases, because the
bandgap for the insulating `I' region is not available. For this reason, the
results for SIS junctions in Fig. 3 in the main text are calculated with
$\epsilon_{r}=7$, so the gaps are preserved, while the results for SNS and
mixed junctions are calculated with $\epsilon_{r}=4$. However, it is worth
noting that, in contrast to SIS JJs, the current in SNS and mixed JJs depends
only weakly on $\epsilon$. We show this explicitly in Fig. S2, which includes
CPRs and critical currents for these JJs, for $\epsilon_{r}=7$. The results
are qualitatively the same as those in Fig. 3 and 4 in the main text.

Figure S2: (a) Current-phase relation in TBG SNS and mixed JJs. (b) Critical
current versus twist angle for all configurations. Everything equal to Fig.
3(a-b) in the main text, but with dielectric constant $\epsilon_{r}=7$ instead
of 4.

To obtain $\rho(n)$, we fit the bandstructure of the tight-binding model to
the continuum model of Ref. Moon  _et al._ (2014) and do a self-consistent
calculation. Fig. S3 compares the resulting bandstructures obtained with the
continuum and tight-binding Hamiltonians including the Hartree term. The bands
are in fair agreement. In particular, note that the very narrow bands near -55
meV, which are similar in both cases, set the Fermi level, and over $95\%$ of
the critical current comes from states in a small window around the gap which
opens at the Fermi level.

Figure S3: (a) Low energy TBG bands for a tight-binding model (red) and a
continuum model, both including Hartree interactions. Filling $n=-2.4$,
$\theta=1.06^{\circ}$ and $\epsilon_{r}=4$. (b) There are four bands because
the system is technically a nanotube, in this case with a unit cell twice the
size that of TBG, so its Brillouin zone results from the folding depicted in
the figure. The green lines are the momentum values allowed by the periodic
boundary condition which closes the nanotube. For a detailed discussion, see
Ref. Sainz-Cruz  _et al._ (2021). Values of $\rho$ for the Hartree term
---
|  1.04∘ |  1.05∘ |  1.06∘ |  1.07∘ |  1.08∘ |  1.09∘ |  1.10∘ |  1.11∘ |  1.12∘ |  1.13∘ |  1.14∘ |  1.15∘ |  1.16∘
$n=-2.4$; $\epsilon=4$ |  -0.802 |  -0.795 |  -0.768 |  -0.778 |  -0.776 |  -0.776 |  -0.77 |  -0.776 |  -0.77 |  -0.773 |  -0.766 |  -0.766 |  -0.769
$n=+2.4$; $\epsilon=4$ |  0.802 |  0.812 |  0.768 |  0.76 |  0.78 |  0.795 |  0.787 |  0.779 |  0.785 |  0.793 |  0.775 |  0.785 |  0.781
$n=-2.4$; $\epsilon=7$ |  -0.962 |  -0.961 |  -0.948 |  -0.946 |  -0.942 |  -0.934 |  -0.927 |  -0.921 |  -0.911 |  -0.9 |  -0.896 |  -0.887 |  -0.876
$n=+2.4$; $\epsilon=7$ |  0.97 |  0.971 |  0.971 |  0.966 |  0.961 |  0.956 |  0.944 |  0.936 |  0.93 |  0.917 |  0.907 |  0.896 |  0.888
Table 2: Calculated values of the constant $\rho$ as a function of twist
angle, filling and dielectric constant.

To obtain the results in Fig. 3 and 4, we have exploited the fact that most of
the critical current comes from states near the superconducting gap $\Delta=1$
meV, in fact, we have observed that states in the window $[-2\Delta,2\Delta]$
carry over 95$\%$ of the current. Therefore, we have approximated Eq. 1 in the
main text as:

$\mathcal{I}\approx\frac{e}{h}\frac{\partial}{\partial\phi}\sum_{\mid\epsilon_{i}\mid<2\Delta}\epsilon_{i}$
(S25)

For the calculation of the current in the mixed junction, spin is explicitly
included in the model, so the Hamiltonian in Eq. 2 is doubled. The $f$-wave
electrode is at filling $n=2.4$ and the $s$-wave electrode at $n=-2.4$. The
interface region interpolates smoothly between these two fillings. The gaps
are similarly smoothed, see Fig. S4, in contrast to SNS and SIS junctions, for
which we use hard boundary conditions.

Figure S4: Superconducting gaps and filling at the interface between $f$-wave
and $s$-wave electrodes in the mixed junction, versus position in units of
moiré period, measured from the center of the junction.

With regards to the lenght of the system, we have found that 24 moiré unit
cells per electrode are enough to reach convergent results, except in the
mixed armchair junction, where 41 cells per electrode are needed. The low-
energy spectrum was obtain with the library ARPACK. The complexity was
approximately $\mathcal{O}(N^{2})$ with N the number of sites in the system.
To verify the algorithm was working as intended, we first reproduced some of
the results in one-dimensional chains obtained with a Green's functions
technique in Ref. Perfetto  _et al._ (2009).

###  .7 Toy model junction

In the main text we showed that the critical current in mixed
$f$-wave/$s$-wave TBG junctions has a phase periodicity of $\pi$, half of
conventional junctions. This phenomenon was also found in Ref. Zazunov and
Egger (2012) in one-dimensional mixed chains. Fig. S5 depicts a toy model
which reproduces the result: in a chain of atoms with spin-unpolarized Kitaev
($p$-wave) pairing Kitaev (2001) on one side and $s$-wave pairing on the
other, the current is $\pi$-periodic. Note that the saw-tooth profiles are a
consequence of fully ballistic transport Ishii (1970); Golubov  _et al._
(2004).

Figure S5: (a) Schematic of a mixed one-dimensional toy model junction, with
Kitaev pairing on one side and $s$-wave pairing on the other. (b) CPR of the
mixed junction, compared to SNS junctions. The parameters used are $t=1$,
$\Delta_{K}=0.1$ and $\Delta_{s}=0.2$. SNS junctions have five metallic atoms
in the link.

###  .8 Andreev spectra

Figure S6: Subgap Andreev spectra in TBG SNS junctions, as a function of the
superconducting phase difference, for $f$-wave and $s$-wave pairings and for
electron (top row) and hole (bottom) domes, at different twist angles.
$\epsilon_{r}=4$.

Figure S6 shows the subgap Andreev spectra in TBG junctions in the SNS
configuration, for $f$-wave and $s$-wave parings at different fillings and
twist angles. As stated in the main text, these states carry most of the
current in SNS junctions. $f$-wave and $s$-wave have similar spectra, except
for the quasi-flat levels that appear for $f$-wave pairing. These states,
which are localized near the edges of the sample, are precursors of Majorana
modes, which will be analyzed in a forthcoming publication. The spectrum
changes fast near the magic angle, compare 1.06∘ and 1.09∘. The fact that the
critical current increases with twist angle in these junctions is seen here as
the growth of the slope of the Andreev levels with angle. There is marked
electron-hole asymmetry. Finally, Fig. S7 shows the subgap spectrum of a mixed
$f$-wave/$s$-wave junction, which has a periodicity of $\pi$.

Figure S7: Subgap Andreev spectra in a mixed $f$-wave/$s$-wave TBG Josephson
junction. $\theta=1.06^{\circ}$. $\epsilon_{r}=4$.