††thanks: These authors contributed equally to this work. ††thanks: These
authors contributed equally to this work.

# Mirror-protected Majorana zero modes in $f$-wave multilayer graphene
superconductors

Võ Tiến Phong1<EMAIL_ADDRESS>Héctor Sainz-Cruz2<EMAIL_ADDRESS>Eugene J. Mele1<EMAIL_ADDRESS>Francisco Guinea2,3
<EMAIL_ADDRESS>1Department of Physics and Astronomy, University of
Pennsylvania, Philadelphia PA 19104, U.S.A.  2IMDEA Nanoscience, Faraday 9,
28049 Madrid, Spain  3Donostia International Physics Center, Paseo Manuel de
Lardizábal 4, 20018 San Sebastián, Spain

(August 28, 2024)

###### Abstract

Inspired by recent experimental discoveries of superconductivity in chirally-
stacked and twisted multilayer graphene, we study models of $f$-wave
superconductivity on the honeycomb lattice with arbitrary numbers of layers.
These models respect a mirror symmetry that allows classification of the bands
by a mirror-projected winding number $\nu_{\pm}$. For odd numbers of layers,
the systems are topologically nontrivial with $\nu_{\pm}=\pm 1$. Along each
mirror-preserving edge in armchair nanoribbons, there are two protected
Majorana zero modes. These modes are present even if the sample is finite in
both directions, such as in rectangular and hexagonal flakes. Crucially, zero
modes can also be confined to vortex cores, which can be created by a magnetic
field or localized magnetic impurities and accessed by local scanning probes.
Finally, we apply these models to twisted bilayer and trilayer systems, which
also feature boundary-projected and vortex-confined zero modes. Since vortices
are experimentally accessible, our study suggests that superconducting
multilayer graphene systems are promising platforms to create and manipulate
Majorana zero modes.

Introduction: Due to their non-Abelian braiding statistics and immunity to
quantum decoherence, Majorana zero modes (MZMs) are highly sought after as
building blocks for topological quantum computation [1, 2, 3]. They are
believed to exist as excitations of the $\nu=5/2$ fractional quantum Hall
effect [4, 5], at the ends of a spinless $p$-wave superconducting chain [6],
or in the vortex cores of spin-triplet $p_{x}+ip_{y}$ superconductors [7, 8,
9]. However, these platforms are challenging to realize experimentally because
a full understanding of the $\nu=5/2$ state remains elusive while spin-triplet
superconductors are scarce in nature [10]. To circumvent these problems, the
modern search for MZM’s focuses primarily on proximitized systems that inherit
the desired superconducting properties from an otherwise ordinary substrate
superconductor [11, 12, 13, 14, 15]. Some of these systems have been
experimentally explored, with various groups claiming to have observed MZMs
due to the presence of zero-bias conductance peaks [16, 17, 18, 19, 20, 21,
22, 23, 24]. However, it is now clear that disorder-induced Andreev bound
states can masquerade as MZMs in conductance experiments [25, 26, 27, 28, 29,
30, 31, 32]. Therefore, quenching disorder is an important goal in the pursuit
of MZMs.

As an exceptionally low-disorder platform [33, 34], graphene is a promising
material for the realization of MZMs. Up to now, graphene-based proposals have
involved the proximity effect since intrinsic superconductivity in graphene
remained elusive for decades [35, 36, 37]. This paradigm was recently shifted
by several groundbreaking experiments showing that graphene multilayers are
robust, highly-tunable superconductors [38, 39, 40, 41, 42, 43, 44, 45, 46,
47, 48]. Moreover, there is mounting evidence that these states involve an
exotic, non $s$-wave pairing [49, 50, 51, 52, 53, 43, 44, 45]. In particular,
spin-triplet, valley-odd $f$-wave pairing is emerging as a leading candidate
for the superconducting symmetry [49, 50, 54, 55]. Inspired by these recent
developments, we study models of intrinsic $f$-wave superconductivity in
chirally-stacked multilayer graphene. We show that in odd-layer
configurations, these systems are topological mirror superconductors
characterized by a nontrivial mirror winding number. Therefore, nanoribbons
with mirror-symmetric edges must host one Majorana mirror pair per boundary.
We also perform calculations on finite-size flakes to show that these zero
modes remain, although their mirror eigenvalue and degeneracy depend
sensitively on the character of the edge terminations. Furthermore, we
calculate the spectrum of finite flakes that host vortices. We find that a
pair of zero modes are confined to each vortex core. This observation holds
significant experimental implications since vortices can be readily created
and probed by existing techniques [56, 57, 58, 59, 60, 61, 62, 63, 64].
Finally, we extend these results to twisted systems, wherein robust zero modes
are found in twisted trilayer graphene, but not twisted bilayer graphene.
Importantly, a pair of zero modes is also trapped at vortices in twisted
trilayer graphene. Our work suggests that multilayer graphene is an
experimentally-feasible, ultra-clean Majorana platform.

Figure 1: Majorana zero modes in a monolayer model. (a) Tight-binding model of
superconducting monolayer graphene. The hoppings between electrons and
electrons (and between holes and holes) are indicated in blue, while the
“hoppings” between electrons and holes are indicated in red. (b) The Brillouin
zone showing the mirror-symmetric line in red along which the winding number
is calculated. (c)-(d) Band structures of mirror-symmetric nanoribbons showing
protected edge states that cross zero energy at $k_{y}=0.$ The edge structures
are shown in the insets. Here, $\Delta=\mu=0.05t_{0},$ and the ribbons are 99a
wide.

Monolayer graphene: We begin with a simple model of spinless superconducting
monolayer graphene wherein the superconducting order parameter in real space
is described by hoppings between the electrons and holes of the same
sublattice with phase winding as indicated in Fig. 1(a). Importantly, the
order parameter changes signs between the $K$ and $K^{\prime}$ valleys,
realizing an exotic $f$-wave superconductor. We will show that at mirror-
symmetric boundaries of nanoribbons, this model supports topologically-
protected pairs of MZMs. The tight-binding Hamiltonian is described in Fig.
1(a). The Bogoliubov-de Gennes (BdG) Hamiltonian in momentum space is given by
$\hat{\mathcal{H}}=\frac{1}{2}\sum_{\mathbf{k}}\hat{\Psi}^{\dagger}(\mathbf{k})\mathcal{H}^{\text{BdG}}(\mathbf{k})\hat{\Psi}(\mathbf{k}),$
where

$\begin{split}\mathcal{H}^{\text{BdG}}(\mathbf{k})&=\tau_{z}\left[\sigma_{0}\mu+\sigma_{x}h_{r}(\mathbf{k})-\sigma_{y}h_{i}(\mathbf{k})\right]+\tau_{x}\sigma_{0}g(\mathbf{k}),\\\
h(\mathbf{k})&=h_{r}(\mathbf{k})+ih_{i}(\mathbf{k})=-t_{0}\left(1+e^{i\kappa_{1}}+e^{i\kappa_{2}}\right),\\\
g(\mathbf{k})&=\Delta\left[\sin\left(\kappa_{1}-\kappa_{2}\right)+\sin\left(-\kappa_{1}\right)+\sin\left(\kappa_{2}\right)\right],\end{split}$
(1)

$\kappa_{i}=\mathbf{k}\cdot\mathbf{a}_{i},$ $\mathbf{a}_{i}$ are the lattice
vectors, the $\tau$ and $\sigma$ Pauli matrices act on the Nambu particle-hole
and sublattice $\\{A,B\\}$ spaces respectively with the subscript $0$ denoting
the identity operator, and
$\hat{\Psi}^{\dagger}(\mathbf{k})=\left(\hat{c}^{\dagger}_{A,\mathbf{k}},\hat{c}^{\dagger}_{B,\mathbf{k}},\hat{c}_{A,-\mathbf{k}},\hat{c}_{B,-\mathbf{k}}\right).$
This Hamiltonian depends only on three real parameters: $t_{0}>0$ the hopping
between nearest neighbors, $\Delta$ the superconducting parameter, and $\mu$
the chemical potential. Assuming that $\Delta\neq 0,$ the bulk band structure
is gapped for $|\mu|<t_{0}$ and $|\mu|>3t_{0}$ 111The latter condition is
trivially true because it simply corresponds to fully-filled or fully-empty
fermionic bands.. We focus only the lightly-doped regime. Our system respects
time reversal $\mathcal{T},$ particle hole $\mathcal{P},$ and mirror
$\mathcal{M}_{y}$ symmetries:

$\begin{split}\mathcal{T}&=\tau_{z}\sigma_{0}\mathcal{K}:\quad\mathcal{T}\mathcal{H}^{\text{BdG}}(\mathbf{k})\mathcal{T}^{-1}=\mathcal{H}^{\text{BdG}}(-\mathbf{k}),\\\
\mathcal{P}&=\tau_{x}\sigma_{0}\mathcal{K}:\quad\mathcal{P}\mathcal{H}^{\text{BdG}}(\mathbf{k})\mathcal{P}^{-1}=-\mathcal{H}^{\text{BdG}}(-\mathbf{k}),\\\
\mathcal{M}_{y}&=\tau_{0}\sigma_{x}:\quad\mathcal{M}_{y}\mathcal{H}^{\text{BdG}}(\mathbf{k})\mathcal{M}_{y}^{-1}=\mathcal{H}^{\text{BdG}}(\mathcal{M}_{y}\mathbf{k}),\\\
\end{split}$ (2)

where $\mathcal{K}$ is the complex conjugation operator. Combining
$\mathcal{P}$ and $\mathcal{T}$ leads to a chiral symmetry $\mathcal{C}$

$\mathcal{C}=i\tau_{y}\sigma_{0}:\quad\mathcal{C}\mathcal{H}^{\text{BdG}}(\mathbf{k})\mathcal{C}^{-1}=-\mathcal{H}^{\text{BdG}}(\mathbf{k}),$
(3)

which requires that every state at $E(\mathbf{k})$ must have a partner at
$-E(\mathbf{k}).$ Since both $\mathcal{T}^{2}=\mathcal{P}^{2}=+1,$ our system
belongs to the topologically-trivial two-dimensional $BDI$ class [66].
Therefore, a generic termination does not guarantee the existence of boundary
states.

With mirror symmetry, we can enrich the topological classification along
mirror-symmetric lines in the Brillouin zone [67, 68, 69]. In our case, there
is only one independent mirror-symmetric line in the Brillouin zone along
$k_{y}=0$ as shown in Fig. 1(b) 222It may appear that there are two mirror-
symmetric lines: one along $k_{y}=0$ and one along $k_{y}=2\pi/\sqrt{3}a,$ but
these two are actually connected to each other by a reciprocal lattice vector.
So they are not independent.. On this line, the Hamiltonian commutes with the
mirror operator
$\left[\mathcal{M}_{y},\mathcal{H}^{\text{BdG}}(k_{x},k_{y}=0)\right]=0,$ so
we can block-diagonalize the Hamiltonian into a mirror-odd and mirror-even
sector. Within each mirror sector, we can further put the Hamiltonian into
chiral off-diagonal form since $\left[\mathcal{M}_{y},\mathcal{C}\right]=0$
[66, 71, 72],

$\mathcal{H}^{\text{BdG}}=\begin{pmatrix}\mathcal{H}_{-}&0\\\
0&\mathcal{H}_{+}\end{pmatrix},\quad\text{where}\quad\mathcal{H}_{\pm}=\begin{pmatrix}0&\mathcal{D}_{\pm}\\\
\mathcal{D}_{\pm}^{\dagger}&0\end{pmatrix}.$ (4)

The mirror-projected Hamiltonians $\mathcal{H}_{\pm}$ are characterized by
winding numbers $\nu_{\pm}=\frac{1}{2\pi i}\oint
dk_{x}\text{Tr}\left[\tilde{\mathcal{D}}_{\pm}^{\dagger}(k_{x})\partial_{k_{x}}\tilde{\mathcal{D}}_{\pm}(k_{x})\right],$
where $\tilde{\mathcal{D}}_{\pm}(k_{x})$ is the obtained from
$\mathcal{D}_{\pm}(k_{x})$ via singular value decomposition [71, 72]. For our
model, we have $\mathcal{D}_{\pm}(k_{x})=\mu\pm h(k_{x})+ig(k_{x})$ and
$\nu_{\pm}=\pm\text{sign}\left(\Delta/t_{0}\right).$ This unity winding number
predicts that when a mirror-symmetric edge is cut parallel to the
$y$-direction, there exists one topologically-protected boundary mode per
mirror sector at $k_{y}=0.$ Importantly, this mode must reside exactly at zero
energy due to chiral symmetry. More generally, the odd parity of $\nu_{\pm}$
guarantees the existence of at least one exact zero mode at $k_{y}=0$. We
illustrate these zero modes for two mirror-symmetric boundaries in Fig.
1(c)-(d). In Fig. 1(c), the termination is pristine armchair, while in Fig.
1(d), the termination is jagged armchair that includes both armchair and
zigzag characters. In both cases, since $\mathcal{M}_{y}$ is preserved, we
find two zero modes, one from each mirror sector.

The Majorana zero modes can be obtained analytically from a continuum theory
for armchair edges. Assuming that $|\Delta|,|\mu|\ll t_{0},$ the relevant
physics is described by a Dirac theory in the original sublattice basis as

$\mathcal{H}^{\text{BdG}}(\mathbf{r})=\tau_{z}\left[\mu\sigma_{0}-i\hbar
v_{F}\nabla_{\mathbf{r}}\cdot\left(\xi\sigma_{x},\sigma_{y}\right)\right]-\xi\tilde{\Delta}\tau_{x}\sigma_{0},$
(5)

where $\hbar v_{F}=\sqrt{3}t_{0}a/2,$ $\tilde{\Delta}=3\sqrt{3}\Delta/2,$ and
$\xi=\pm$ denotes valley. Putting the edge at $x=0$ and extending into
$x\rightarrow-\infty,$ the (unnormalized but normalizable) zero-energy
solutions at $k_{y}=0$ that satisfy the armchair boundary conditions are

$\psi_{\xi,m_{y}}=e^{x/\ell-
iq_{x}x}\xi\begin{pmatrix}1,&m_{y},&-im_{y}\text{sign}\Delta,&-i\text{sign}\Delta\end{pmatrix}^{T},$
(6)

where the decay length is $\ell=\hbar v_{F}/|\tilde{\Delta}|,$ the wavelength
is $q_{x}=m_{y}\xi\mu/\hbar v_{F},$ and $m_{y}=\pm 1$ is the mirror eigenvalue
of the mode. The valley-antisymmetric nature of the superconducting gap is
crucial to the existence of the zero modes in Eq. (6) because had the gap been
endowed with the same sign in both valleys, these modes would not have been
normalizeable. Thus, the presence of armchair-confined Majorana states can
serve as a diagnostic of valley-odd, spin-triplet superconductivity in
graphene [73]. Away from $k_{y}=0,$ the mirror-odd and even sectors hybridize,
lifting the energy away from zero. To find the dispersion, we can rewrite Eq.
(5) at finite small $k_{y}$ in the basis of Eq. (6)

$\mathcal{H}^{\text{BdG}}(k_{y})=\hbar\tilde{v}_{F}\begin{pmatrix}0&ik_{y}\\\
-ik_{y}&0\end{pmatrix}.$ (7)

Thus, the mid-gap boundary-projected dispersion is
$\mathcal{E}=\pm\hbar\tilde{v}_{F}k_{y},$ where
$\tilde{v}_{F}=v_{F}[1+(\mu/\tilde{\Delta})^{2}]^{-1}$, and the associated
eigenfunctions are $\psi_{\xi,\mathcal{E}=\pm}=e^{ik_{y}y}\left[\mp
i\psi_{\xi,m_{y}=-1}+\psi_{\xi,m_{y}=+1}\right].$

For potential applications in topological quantum computing, it is desirable
to isolate Majorana zero modes spatially. For our model, the Majorana zero
modes always come in pairs as demanded by $\mathcal{M}_{y}$ along
translationally-invariant edges. One way to isolate a Majorana zero mode is to
break $\mathcal{M}_{y}$ symmetry by adding a spatially-dependent $\sigma_{z}$
mass of the form $V(x,y)=m(y)\tau_{z}\sigma_{z}$ where
$\lim_{y\rightarrow\pm\infty}m(y)=\pm m$ and $m$ is a constant. The mass
changes sign at $y=0.$ Then, there are two possible zero-energy solutions
$\phi_{\xi,m_{y}}=e^{m_{y}\int^{y}m(y^{\prime})dy^{\prime}/\hbar
v_{F}}\psi_{\xi,m_{y}}$ 333We continue to label the $\phi_{\xi,m_{y}}$ states
with mirror eigenvalue $m_{y}$ only as a matter of convenience. The reader
should not take this as implying that mirror symmetry is preserved. It is
indeed broken explicitly by the potential $V(x,y)$. However, only one of them
is normalizable depending on whether $m$ is positive or negative. Thus, the
effect of the mirror-breaking mass is to gap out the two boundary modes and
create a single bound state localized where the mass changes sign.

Figure 2: Isolating Majorana zero modes on finite flakes. (a) Majorana zero
mode localized on the left armchair edge of a rectangular flake with different
values of the chemical potential and different terminations of the horizontal
top and bottom edges (either zigzag or bearded). The mirror eigenvalue of the
shown state, $m_{y},$ is indicated. (b) Dependence of the zero mode on
vertical length of the flake. (c) Majorana corner state in a hexagonal flake.
There are six such states in total. In these simulations, $\Delta=0.1t_{0}$
and $|\mu|=0.2t_{0}.$ (d) One of two vortex states confined to the vortex core
near zero energy $E\approx 10^{-17}t_{0}$, with $\Delta=0.2t_{0},$
$\mu=0.1t_{0},$ and $d\approx 0.$ The irregular shape has $9,609$ atoms, and
is chosen to illustrate that the vortex states do not depend on boundary
conditions. The energies of this flake are shown in (e), with the zero modes
emphasized. The vortex is centered at a hexagon’s center. We only show the
electron density.

Another way to isolate single zero modes is by careful design of edges and
corners. For instance, in a rectangular flake that features armchair edges on
the left and the right and either zigzag or bearded edges on both the top and
the bottom, there are two zero modes that are localized on the left and right
armchair edges. These modes are odd (even) under mirror if $\mu>0$ and the
horizontal edges are zigzag (bearded) or if $\mu<0$ and the horizontal edges
are bearded (zigzag), as shown in Fig. 2(a). Therefore, the imposition of
zigzag and bearded edges mimics a large $\sigma_{z}$ mass above the sample
that goes to zero inside the sample and then changes sign below the sample.
While the width separating the two armchair edges needs to be large to prevent
mixing of the zero modes, the length of the flake along the $y$ direction does
not qualitatively affect these modes, as shown in Fig. 2(b). Furthermore, in a
hexagonal flake, we find six localized corner states at zero energy, one of
which is shown in Fig. 2(c). The hexagonal shape can be thought of as an
example of the rule that most clean, locally straight graphene edges are
described, in the continuum limit, by the zigzag boundary conditions [75].
Edge regions connecting $A-$ and $B-$ like boundaries must include armchair
like segments, where zero modes will be localized. Examples of zero modes in a
distorted hexagon are shown in Ref. [76]. It is worth emphasizing that the
chosen simulation parameters are unrealistically large to facilitate fast
numerical convergence, but should not qualitatively affect the topologically-
derived conclusions. As demonstrated, the current platform presents many
avenues to manipulate the Majorana zero modes.

Crucially, vortices can also support zero-energy modes. We calculate these
vortex states numerically on finite flakes using the tight-binding framework,
where the gap function is modified to
$\Delta(\mathbf{r}_{i},\mathbf{r}_{j})=\pm
e^{i\phi_{\mathbf{r}}}\Delta\tanh(|\mathbf{r}|/d)/2i,$
$\mathbf{r}=(\mathbf{r}_{i}+\mathbf{r}_{j})/2$ is the center-of-mass position,
$\phi_{\mathbf{r}}$ is the angular coordinate of $\mathbf{r},$ and $d$ is the
coherence length. When the vortex is centered along a
$\mathcal{M}_{y}$-invariant line that connects a hexagon’s center and a mid-
bond, we find two vortex-confined modes at numerically-exact zero energies, as
shown in Fig. 2(d). This is strong evidence that these vortex-confined zero
modes are protected by mirror symmetry. When the vortex is centered elsewhere,
such as on a carbon site, then these low-energy modes are no longer at zero
energy. However, when $d\gtrsim a,$ the precise origin of the vortex becomes
less important, and we find near-zero modes for various geometries
numerically. The coherence length is the ratio of the Fermi velocity and the
superconducting gap, which for practically-relevant small gaps, is many times
the lattice constant. We leave the topological stability of these modes,
including in the presence of disorder, to future analyses.

Rhombohedral multilayer graphene: We now generalize the previous results to
rhombohedral $N$-layer graphene stacks, which is motivated by the recently-
discovered superconductivity in Bernal bilayer graphene and rhombohedral
trilayer graphene [43, 45, 47, 48]. The Hamiltonian is now modified to

$\mathcal{H}^{\text{BdG}}(\mathbf{k})=\tau_{z}\mathcal{H}_{N}(\mathbf{k})+\tau_{x}g(\mathbf{k}),$
(8)

where

$\begin{split}\mathcal{H}_{N}(\mathbf{k})&=\mathbb{I}_{N}\mathcal{H}_{\text{intra}}(\mathbf{k})+\mathbb{U}_{N}\mathcal{H}_{\text{inter}}(\mathbf{k})+\mathbb{L}_{N}\mathcal{H}^{\dagger}_{\text{inter}}(\mathbf{k}),\\\
\mathcal{H}_{\text{intra}}(\mathbf{k})&=\begin{pmatrix}\mu&h(\mathbf{k})\\\
h^{*}(\mathbf{k})&\mu\end{pmatrix},\quad\mathcal{H}_{\text{inter}}(\mathbf{k})=\begin{pmatrix}0&0\\\
\gamma_{1}&0\end{pmatrix},\end{split}$ (9)

$\mathbb{I}_{N}$ is the $N\times N$ identity matrix, $\mathbb{U}_{N}$
($\mathbb{L}_{N}$) is the $N\times N$ matrix with ones along the diagonal
above (below) the principal diagonal, acting on layer space, and $\gamma_{1}$
is the interlayer hopping. For multiple layers, mirror symmetry is actually
$\mathcal{C}_{2x}$ symmetry that simultaneously exchanges layers and
sublattices represented by
$\mathcal{M}_{y}=\tau_{0}\bar{\mathbb{I}}_{N}\sigma_{x},$ where
$\bar{\mathbb{I}}_{N}$ is the $N\times N$ anti-diagonal matrix of ones. For
$N$ even, all layers are exchanged under mirror, while for $N$ odd, there is
one central layer which is mapped onto itself under mirror. Due to this mirror
symmetry, we can classify the bands by the same topological invariant as
before. For $N$ even, the system is topologically trivial with $\nu_{\pm}=0.$
Consequently, there are no Majorana zero modes pinned to $k_{y}=0$ on armchair
nanoribbons, as shown for $N=2,4$ in Fig. 3(a). Interestingly, we find that
for the current models, there are still zero modes displaced away from
$k_{y}=0.$ However, it is important to emphasize that these $k_{y}\neq 0$
modes are not topologically protected, and thus they can be removed by
perturbations which are invariant under $\mathcal{M}_{y}$ and do not close the
bulk gaps. One such perturbation is a layer-antisymmetric $\sigma_{z}$
potential energy.

On the contrary, for $N$ odd, the systems are topological, protected by
$\mathcal{M}_{y}$, with $\nu_{\pm}=\pm 1.$ Therefore, there are two
topological zero modes per armchair edge pinned at $k_{y}=0,$ in addition to
many removable accidental zero modes at nonzero momenta, as shown for $N=3,5$
in Fig. 3(a). In addition to these zero modes, we also find numerous mid-gap
states for $N>1.$ Like the non-topological zero modes, these mid-gap states
are also removable; however, they are generically present at armchair edges.
The topological zero modes can be removed by breaking $\mathcal{M}_{y}$. In
the multilayer case, one can break $\mathcal{M}_{y}$ by invalidating layer
equivalence with a perpendicular electric field, which is much more
experimentally accessible than a staggered potential. Regarding vortices, for
the trilayer case, we have also found numerically a number of low-energy
states confined to a vortex core. In fact, we find a pair of numerically-exact
zero modes when the vortex center is located along a
$\mathcal{M}_{y}$-symmetric line, as shown in Fig. 3(b). Like before, if
$d\gtrsim a,$ then the vortex center is not as important, so we find these
near-zero modes to be robust. We expect all the odd-layer configurations to
feature vortex-confined zero modes, but we have not numerically calculated
these due to prohibitive computational cost as the number of layers increases.

Figure 3: Rhombohedral multilayer graphene. (a) Band structures for armchair
nanoribbons with $N=2,3,4,5$ layers. For $N=2,4,$ there are no zero modes at
$k_{y}=0,$ but there are removable zero modes displaced to nonzero $k_{y}.$
For $N=3,5,$ there are two topological zero modes per armchair edge pinned to
$k_{y}=0.$ Here, $\Delta=\mu=0.05t_{0},$ $\gamma_{1}=0.1t_{0},$ and the
ribbons are $99a$ wide. (b) Vortex-confined zero mode, with $E\approx
10^{-15}t_{0},$ in the three layers of rhombohedral trilayer graphene with
$\mu=0.1t_{0},$ $\Delta=0.2t_{0},$ and $d\approx 0.$ The flake contains
$21,600$ atoms, and the vortex is centered at a mid-bond in the middle layer.
Figure 4: Twisted multilayer graphene. (a, b) Band structures for twisted
bilayer graphene nanoribbons, in the hole superconducting dome (filling
$n=-2.4$), with a twist angle of (a) $\theta=1.06^{\circ}$ and (b)
$1.16^{\circ}$. The width of the ribbon is $W\approx 35$ moiré periods
($L_{M}$) and the superconducting gap is $\Delta_{\text{sc}}=1$ meV. (c)
Charge map of a TBG zero mode, showing localization at the last full moiré
near an edge. (d) Band structure of a twisted trilayer graphene nanoribbon
with armchair edges in the top and bottom layers. The states in red come from
the effective monolayer sector and thus have all the charge in top and bottom
layers. Here $\theta=5^{\circ}$, $\mu=0$, $\Delta_{\text{sc}}=200$ meV, [77],
and $|T|>L_{M}$ due to a commensurability condition [76]. (e) Charge map of a
TTG zero mode at $k=0$ in each of the three layers near the left edge of the
ribbon. (f) Zero mode confined to a vortex core in twisted trilayer graphene,
with $E\approx 10^{-10}t_{0}$. Here $\theta=1.5^{\circ}$, $n=-2.4$,
$\Delta_{\text{sc}}=0.05t_{0}$, and the flake contains $323,742$ atoms.

ABA multilayer graphene: The above results can be generalized to $ABA$ stacks
as well. When the number of layers is odd, such a stack can be decomposed into
a direct sum of bilayer sectors and one monolayer sector due to the presence
of $\mathcal{M}_{z}$ symmetry [78]. The monolayer sector behaves exactly like
the monolayer toy model albeit it is written in a basis of a coherent
superposition of layers. As such, in these systems, there are also protected
Majorana zero modes in the monolayer sector at armchair boundaries and robust
vortex-confined low-energy modes. Although superconductivity has not yet been
experimentally observed in $ABA$ multilayer graphene, these observations
suggest that twisted multilayer graphene will contain the same exotic physics
owning to the decomposition into a monolayer sector and bilayer sectors.

Twisted multilayer graphene: Inspired by the aforementioned findings, we now
look for zero modes in twisted multilayer graphene, using a scaled tight-
binding model [79, 80, 81, 82, 54, 83, 84, 85, 76]. Results for twisted
bilayer and trilayer graphene appear in Fig. 4. In panels (a, b), the band
structure of TBG nanoribbons is shown at two different twist angles. For these
angles, which are close to the magic angle, the bands are qualitatively
similar to those of Bernal bilayer graphene. There are 8 zero modes at
$k_{y}\neq 0$, which are not topologically protected. These modes and the rest
of sub-gap Andreev states appear near the edges of the system, at the last
complete moiré AA regions, see Fig. 4(c), and can be observed via Scanning
Tunneling Microscopy [49, 50]. The bands depend sensitively on twist angle.
Increasing the angle decreases the momentum of the near-zero modes, compare
Fig. 4(a) and (b). When the angle is increased further, the zero modes
disappear and only Andreev states near $k_{y}=0$ remain. Few sub-gap states
survive at $\theta\gtrsim 1.2^{\circ}$ [76].

Just like in TBG, in twisted trilayer graphene (TTG) [41, 42], generic
$\mathcal{M}_{y}$-symmetric edges do not lead to topologically-protected zero
modes at $k_{y}=0$. However, a new possibility arises due to the decomposition
of TTG into TBG plus monolayer graphene, as found in Ref. [78]. In particular,
since the effective monolayer comes from the odd combination of top and bottom
layers, when these two layers have armchair edges, the system includes the 4
zero modes at $k_{y}=0$, like monolayer graphene, as shown in Fig. 4(d). The
charge of these zero modes is indeed evenly distributed in top and bottom
layers, with no charge in the middle layer, as can be seen in Fig. 4(e). Since
a decomposition including a monolayer is possible for any alternating-twist
stack with an odd number of layers [78], we expect that 4 zero modes at
$k_{y}=0$ will exist in all such multilayers. Finally, Fig. 4(f) shows one of
the two zero modes confined to a vortex core in twisted trilayer graphene. As
in the ribbon geometry, these zero modes come from the effective monolayer
sector. These promising results suggest that there is a whole family of robust
superconductors in twisted odd-layer graphene that host topologically-
protected Majorana zero modes.

Discussion and conclusion: Superconductivity with $f$-wave pairing is favored
in two-dimensional materials with Fermi pockets at opposing corners of the
Brillouin zone and strong repulsive electron-electron interactions. Defects
which scatter electrons between the different pockets lead, in the
superconducting phase, to Andreev resonances deep within the gap, and even
allow for the existence of isolated Majorana states. In particular, the
Majorana states are present when $\mathcal{M}_{y}$ symmetry is respected. This
symmetry can be broken by applying a perpendicular displacement field. In
rhombohedral stacks, superconductivity has not yet been seen without such a
field; so our model might not immediately apply to these cases [43, 45, 47,
48]. However, superconductivity has been unambiguously established in twisted
trilayer graphene even in the absence of a displacement field [41, 53, 42].
Therefore, twisted trilayer graphene is the most promising candidate to detect
these Majorana bound states, the observation of which is a simple way to
identify this exotic $f$-wave superconductivity. In these materials, edges and
corners are prototypical defects that trap Andreev states and Majorana zero
modes. Although precision fabrication of these boundaries is technically
challenging, there has been much experimental progress in achieving such a
feat [86, 87, 88, 89, 90]. On the other hand, vortices are much more
experimentally accessible because they can be created in the bulk by a
magnetic field or by magnetic impurities. The electronic states at vortex
cores have been extensively studied using local probes [56, 57, 58, 59, 60,
61]. These techniques can also be used to manipulate vortices [62, 63, 64].
Therefore, the vortex-confined zero modes in these graphene-based topological
superconductors can be created, probed, and manipulated by readily-available
technology, opening the way to new types of quantum devices.

We thank Fernando de Juan and Pierre A. Pantaleón for insightful discussions.
V. T. P. and E. J. M. are supported by the U.S. Department of Energy under
grant DE-FG02-84ER45118. IMDEA Nanociencia acknowledges support from the
“Severo Ochoa” Programme for Centres of Excellence in R&D (CEX2020-001039-S /
AEI / 10.13039/501100011033). F. G. and H. S.-C. acknowledge funding from the
European Commission, within the Graphene Flagship, Core 3, grant number 881603
and from grants NMAT2D (Comunidad de Madrid, Spain), SprQuMat (Ministerio de
Ciencia e Innovación, Spain) and financial support through the (MAD2D-CM)-MRR
MATERIALES AVANZADOS-IMDEA-NC.

## References

  * Ivanov [2001] D. A. Ivanov, Non-abelian statistics of half-quantum vortices in $\mathit{p}$-wave superconductors, Phys. Rev. Lett. 86, 268 (2001).
  * Kitaev [2003] A. Y. Kitaev, Fault-tolerant quantum computation by anyons, Annals of physics 303, 2 (2003).
  * Nayak  _et al._ [2008] C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. Das Sarma, Non-abelian anyons and topological quantum computation, Rev. Mod. Phys. 80, 1083 (2008).
  * Moore and Read [1991] G. Moore and N. Read, Nonabelions in the fractional quantum hall effect, Nuclear Physics B 360, 362 (1991).
  * Read and Green [2000] N. Read and D. Green, Paired states of fermions in two dimensions with breaking of parity and time-reversal symmetries and the fractional quantum hall effect, Phys. Rev. B 61, 10267 (2000).
  * Kitaev [2001] A. Y. Kitaev, Unpaired majorana fermions in quantum wires, Physics-Uspekhi 44, 131 (2001).
  * Stern  _et al._ [2004] A. Stern, F. von Oppen, and E. Mariani, Geometric phases and quantum entanglement as building blocks for non-abelian quasiparticle statistics, Phys. Rev. B 70, 205338 (2004).
  * Das Sarma  _et al._ [2006] S. Das Sarma, C. Nayak, and S. Tewari, Proposal to stabilize and detect half-quantum vortices in strontium ruthenate thin films: Non-abelian braiding statistics of vortices in a ${p}_{x}+i{p}_{y}$ superconductor, Phys. Rev. B 73, 220502 (2006).
  * Kraus  _et al._ [2009] Y. E. Kraus, A. Auerbach, H. A. Fertig, and S. H. Simon, Majorana fermions of a two-dimensional ${p}_{x}+i{p}_{y}$ superconductor, Phys. Rev. B 79, 134515 (2009).
  * Das Sarma [2023] S. Das Sarma, In search of majorana, Nature Physics 19, 165 (2023).
  * Fu and Kane [2008] L. Fu and C. L. Kane, Superconducting proximity effect and majorana fermions at the surface of a topological insulator, Physical Review Letters 100, 10.1103/physrevlett.100.096407 (2008).
  * Sau  _et al._ [2010] J. D. Sau, R. M. Lutchyn, S. Tewari, and S. Das Sarma, Generic new platform for topological quantum computation using semiconductor heterostructures, Phys. Rev. Lett. 104, 040502 (2010).
  * Lutchyn  _et al._ [2010] R. M. Lutchyn, J. D. Sau, and S. Das Sarma, Majorana fermions and a topological phase transition in semiconductor-superconductor heterostructures, Phys. Rev. Lett. 105, 077001 (2010).
  * Oreg  _et al._ [2010] Y. Oreg, G. Refael, and F. von Oppen, Helical liquids and majorana bound states in quantum wires, Phys. Rev. Lett. 105, 177002 (2010).
  * Alicea [2010] J. Alicea, Majorana fermions in a tunable semiconductor device, Phys. Rev. B 81, 125318 (2010).
  * Mourik  _et al._ [2012] V. Mourik, K. Zuo, S. M. Frolov, S. Plissard, E. P. Bakkers, and L. P. Kouwenhoven, Signatures of majorana fermions in hybrid superconductor-semiconductor nanowire devices, Science 336, 1003 (2012).
  * Das  _et al._ [2012] A. Das, Y. Ronen, Y. Most, Y. Oreg, M. Heiblum, and H. Shtrikman, Zero-bias peaks and splitting in an al–inas nanowire topological superconductor as a signature of majorana fermions, Nature Physics 8, 887 (2012).
  * Deng  _et al._ [2012] M. Deng, C. Yu, G. Huang, M. Larsson, P. Caroff, and H. Xu, Anomalous zero-bias conductance peak in a nb–insb nanowire–nb hybrid device, Nano letters 12, 6414 (2012).
  * Finck  _et al._ [2013] A. D. K. Finck, D. J. Van Harlingen, P. K. Mohseni, K. Jung, and X. Li, Anomalous modulation of a zero-bias peak in a hybrid nanowire-superconductor device, Phys. Rev. Lett. 110, 126406 (2013).
  * Churchill  _et al._ [2013] H. O. H. Churchill, V. Fatemi, K. Grove-Rasmussen, M. T. Deng, P. Caroff, H. Q. Xu, and C. M. Marcus, Superconductor-nanowire devices from tunneling to the multichannel regime: Zero-bias oscillations and magnetoconductance crossover, Phys. Rev. B 87, 241401 (2013).
  * Nadj-Perge  _et al._ [2014] S. Nadj-Perge, I. K. Drozdov, J. Li, H. Chen, S. Jeon, J. Seo, A. H. MacDonald, B. A. Bernevig, and A. Yazdani, Observation of majorana fermions in ferromagnetic atomic chains on a superconductor, Science 346, 602 (2014).
  * Chen  _et al._ [2017] J. Chen, P. Yu, J. Stenger, M. Hocevar, D. Car, S. R. Plissard, E. P. Bakkers, T. D. Stanescu, and S. M. Frolov, Experimental phase diagram of zero-bias conductance peaks in superconductor/semiconductor nanowire devices, Science advances 3, e1701476 (2017).
  * Vaitiekėnas  _et al._ [2021] S. Vaitiekėnas, Y. Liu, P. Krogstrup, and C. Marcus, Zero-bias peaks at zero magnetic field in ferromagnetic hybrid nanowires, Nature Physics 17, 43 (2021).
  * Wang  _et al._ [2022] Z. Wang, H. Song, D. Pan, Z. Zhang, W. Miao, R. Li, Z. Cao, G. Zhang, L. Liu, L. Wen, R. Zhuo, D. E. Liu, K. He, R. Shang, J. Zhao, and H. Zhang, Plateau regions for zero-bias peaks within 5% of the quantized conductance value $2{e}^{2}/h$, Phys. Rev. Lett. 129, 167702 (2022).
  * Liu  _et al._ [2012] J. Liu, A. C. Potter, K. T. Law, and P. A. Lee, Zero-bias peaks in the tunneling conductance of spin-orbit-coupled superconducting wires with and without majorana end-states, Phys. Rev. Lett. 109, 267002 (2012).
  * Roy  _et al._ [2013] D. Roy, N. Bondyopadhaya, and S. Tewari, Topologically trivial zero-bias conductance peak in semiconductor majorana wires from boundary effects, Phys. Rev. B 88, 020502 (2013).
  * Chen  _et al._ [2019] J. Chen, B. D. Woods, P. Yu, M. Hocevar, D. Car, S. R. Plissard, E. P. A. M. Bakkers, T. D. Stanescu, and S. M. Frolov, Ubiquitous non-majorana zero-bias conductance peaks in nanowire devices, Phys. Rev. Lett. 123, 107703 (2019).
  * Prada  _et al._ [2020] E. Prada, P. San-Jose, M. W. A. de Moor, A. Geresdi, E. J. H. Lee, J. Klinovaja, D. Loss, J. Nygård, R. Aguado, and L. P. Kouwenhoven, From andreev to majorana bound states in hybrid superconductor–semiconductor nanowires, Nature Reviews Physics 2, 575 (2020).
  * Pan and Das Sarma [2020] H. Pan and S. Das Sarma, Physical mechanisms for zero-bias conductance peaks in majorana nanowires, Phys. Rev. Res. 2, 013377 (2020).
  * Pan  _et al._ [2021] H. Pan, C.-X. Liu, M. Wimmer, and S. Das Sarma, Quantized and unquantized zero-bias tunneling conductance peaks in majorana nanowires: Conductance below and above $2{e}^{2}/h$, Phys. Rev. B 103, 214502 (2021).
  * Das Sarma and Pan [2021] S. Das Sarma and H. Pan, Disorder-induced zero-bias peaks in majorana nanowires, Phys. Rev. B 103, 195158 (2021).
  * Valentini  _et al._ [2022] M. Valentini, M. Borovkov, E. Prada, S. Martí-Sánchez, M. Botifoll, A. Hofmann, J. Arbiol, R. Aguado, P. San-Jose, and G. Katsaros, Majorana-like coulomb spectroscopy in the absence of zero-bias peaks, Nature 612, 442 (2022).
  * Peres [2010] N. M. R. Peres, Colloquium: The transport properties of graphene: An introduction, Rev. Mod. Phys. 82, 2673 (2010).
  * Das Sarma  _et al._ [2011] S. Das Sarma, S. Adam, E. H. Hwang, and E. Rossi, Electronic transport in two-dimensional graphene, Rev. Mod. Phys. 83, 407 (2011).
  * San-Jose  _et al._ [2015] P. San-Jose, J. Lado, R. Aguado, F. Guinea, and J. Fernández-Rossier, Majorana zero modes in graphene, Physical Review X 5, 10.1103/physrevx.5.041042 (2015).
  * Peñaranda  _et al._ [2023] F. Peñaranda, R. Aguado, E. Prada, and P. San-Jose, Majorana bound states in encapsulated bilayer graphene, SciPost Physics 14, 10.21468/scipostphys.14.4.075 (2023).
  * Xie  _et al._ [2023] Y.-M. Xie, E. Lantagne-Hurtubise, A. F. Young, S. Nadj-Perge, and J. Alicea, Gate-defined topological josephson junctions in bernal bilayer graphene, arXiv 10.48550/ARXIV.2304.11807 (2023).
  * 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).
  * Zhou  _et al._ [2021] H. Zhou, T. Xie, T. Taniguchi, K. Watanabe, and A. F. Young, Superconductivity in rhombohedral trilayer graphene, Nature 598, 434 (2021).
  * 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).
  * 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).
  * Zhang  _et al._ [2022] 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 (2022).
  * Zhang  _et al._ [2023] Y. Zhang, R. Polski, A. Thomson, É. Lantagne-Hurtubise, C. Lewandowski, H. Zhou, K. Watanabe, T. Taniguchi, J. Alicea, and S. Nadj-Perge, Enhanced superconductivity in spin–orbit proximitized bilayer graphene, Nature 613, 268 (2023).
  * Holleis  _et al._ [2023] L. Holleis, C. L. Patterson, Y. Zhang, H. M. Yoo, H. Zhou, T. Taniguchi, K. Watanabe, S. Nadj-Perge, and A. F. Young, Ising superconductivity and nematicity in bernal bilayer graphene with strong spin orbit coupling, Arxiv 10.48550/arXiv.2303.00742 (2023).
  * Oh  _et al._ [2021] M. Oh, K. P. Nuckolls, D. Wong, R. L. Lee, X. Liu, K. Watanabe, T. Taniguchi, and A. Yazdani, Evidence for unconventional superconductivity in twisted bilayer graphene, Nature 600, 240 (2021).
  * Kim  _et al._ [2022] H. Kim, Y. Choi, C. Lewandowski, A. Thomson, Y. Zhang, R. Polski, K. Watanabe, T. Taniguchi, J. Alicea, and S. Nadj-Perge, Evidence for unconventional superconductivity in twisted trilayer graphene, Nature 606, 494 (2022).
  * Lin  _et al._ [2022] J.-X. Lin, P. Siriviboon, H. D. Scammell, S. Liu, D. Rhodes, K. Watanabe, T. Taniguchi, J. Hone, M. S. Scheurer, and J. Li, Zero-field superconducting diode effect in small-twist-angle trilayer graphene, Nature Physics 18, 1221 (2022).
  * Cao  _et al._ [2021a] Y. Cao, J. M. Park, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Pauli-limit violation and re-entrant superconductivity in moiré graphene, Nature 595, 526 (2021a).
  * Cao  _et al._ [2021b] Y. Cao, J. M. Park, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Large pauli limit violation and reentrant superconductivity in magic-angle twisted trilayer graphene, arXiv preprint arXiv:2103.12083 (2021b).
  * Sainz-Cruz  _et al._ [2022] H. Sainz-Cruz, P. A. Pantaleón, V. T. Phong, A. Jimeno-Pozo, and F. Guinea, Junctions and superconducting symmetry in twisted bilayer graphene, arXiv 10.48550/ARXIV.2211.11389 (2022), Phys. Rev. Lett., in press.
  * Crépel  _et al._ [2022] V. Crépel, T. Cea, L. Fu, and F. Guinea, Unconventional superconductivity due to interband polarization, Phys. Rev. B 105, 094506 (2022).
  * Suderow  _et al._ [2014] H. Suderow, I. Guillamón, J. G. Rodrigo, and S. Vieira, Imaging superconducting vortex cores and lattices with a scanning tunneling microscope, Superconductor Science and Technology 27, 063001 (2014).
  * Berthod  _et al._ [2017] C. Berthod, I. Maggio-Aprile, J. Bruér, A. Erb, and C. Renner, Observation of caroli–de gennes–matricon vortex states in ${\mathrm{yba}}_{2}{\mathrm{cu}}_{3}{\mathrm{o}}_{7-\delta}$, Phys. Rev. Lett. 119, 237001 (2017).
  * Kong  _et al._ [2019] L. Kong, S. Zhu, M. Papaj, H. Chen, L. Cao, H. Isobe, Y. Xing, W. Liu, D. Wang, P. Fan, Y. Sun, S. Du, J. Schneeloch, R. Zhong, G. Gu, L. Fu, H.-J. Gao, and H. Ding, Half-integer level shift of vortex bound states in an iron-based superconductor, Nature Physics 15, 1181 (2019).
  * Chen  _et al._ [2021] X. Chen, W. Duan, X. Fan, W. Hong, K. Chen, H. Yang, S. Li, H. Luo, and H.-H. Wen, Friedel oscillations of vortex bound states under extreme quantum limit in ${\mathrm{kca}}_{2}{\mathrm{fe}}_{4}{\mathrm{as}}_{4}{\mathrm{f}}_{2}$, Phys. Rev. Lett. 126, 257002 (2021).
  * Zhang  _et al._ [2021] T. Zhang, W. Bao, C. Chen, D. Li, Z. Lu, Y. Hu, W. Yang, D. Zhao, Y. Yan, X. Dong, Q.-H. Wang, T. Zhang, and D. Feng, Observation of distinct spatial distributions of the zero and nonzero energy vortex modes in $({\mathrm{li}}_{0.84}{\mathrm{fe}}_{0.16})\mathrm{OHFeSe}$, Phys. Rev. Lett. 126, 127001 (2021).
  * Li  _et al._ [2023] C. Li, Y.-F. Zhao, A. Vera, O. Lesser, H. Yi, S. Kumari, Z. Yan, C. Dong, T. Bowen, K. Wang, H. Wang, J. L. Thompson, K. Watanabe, T. Taniguchi, D. R. Hickey, Y. Oreg, J. A. Robinson, C.-Z. Chang, and J. Zhu, Proximity-induced superconductivity in epitaxial topological insulator/graphene/gallium heterostructures, Nature Materials 22, 570 (2023).
  * Ge  _et al._ [2016] J.-Y. Ge, V. N. Gladilin, J. Tempere, C. Xue, J. T. Devreese, J. V. de Vondel, Y. Zhou, and V. V. Moshchalkov, Nanoscale assembly of superconducting vortices with scanning tunnelling microscope tip, Nature Communications 7, 10.1038/ncomms13880 (2016).
  * Polshyn  _et al._ [2019] H. Polshyn, T. Naibert, and R. Budakian, Manipulating multivortex states in superconducting structures, Nano Letters 19, 5476 (2019).
  * Naibert  _et al._ [2021] T. R. Naibert, H. Polshyn, R. Garrido-Menacho, M. Durkin, B. Wolin, V. Chua, I. Mondragon-Shem, T. Hughes, N. Mason, and R. Budakian, Imaging and controlling vortex dynamics in mesoscopic superconductor–normal-metal–superconductor arrays, Phys. Rev. B 103, 224526 (2021).
  * Note [1] The latter condition is trivially true because it simply corresponds to fully-filled or fully-empty fermionic bands.
  * Chiu  _et al._ [2016] C.-K. Chiu, J. C. Y. Teo, A. P. Schnyder, and S. Ryu, Classification of topological quantum matter with symmetries, Rev. Mod. Phys. 88, 035005 (2016).
  * Zhang  _et al._ [2013a] F. Zhang, C. L. Kane, and E. J. Mele, Topological mirror superconductivity, Phys. Rev. Lett. 111, 056403 (2013a).
  * Zhang  _et al._ [2013b] F. Zhang, C. L. Kane, and E. J. Mele, Time-reversal-invariant topological superconductivity and majorana kramers pairs, Phys. Rev. Lett. 111, 056402 (2013b).
  * Ueno  _et al._ [2013] Y. Ueno, A. Yamakage, Y. Tanaka, and M. Sato, Symmetry-protected majorana fermions in topological crystalline superconductors: Theory and application to ${\text{sr}}_{2}{\text{ruo}}_{4}$, Phys. Rev. Lett. 111, 087002 (2013).
  * Note [2] It may appear that there are two mirror-symmetric lines: one along $k_{y}=0$ and one along $k_{y}=2\pi/\sqrt{3}a,$ but these two are actually connected to each other by a reciprocal lattice vector. So they are not independent.
  * Neupert and Schindler [2018] T. Neupert and F. Schindler, Topological crystalline insulators, Topological Matter: Lectures from the Topological Matter School 2017 , 31 (2018).
  * Lin  _et al._ [2021] L. Lin, Y. Ke, and C. Lee, Real-space representation of the winding number for a one-dimensional chiral-symmetric topological insulator, Phys. Rev. B 103, 224208 (2021).
  * Kokkeler  _et al._ [2023] T. H. Kokkeler, C. Huang, F. S. Bergeret, and I. V. Tokatly, Spectroscopic signature of spin triplet odd-valley superconductivity in two-dimensional materials, arXiv 10.48550/ARXIV.2304.14849 (2023).
  * Note [3] We continue to label the $\phi_{\xi,m_{y}}$ states with mirror eigenvalue $m_{y}$ only as a matter of convenience. The reader should not take this as implying that mirror symmetry is preserved. It is indeed broken explicitly by the potential $V(x,y)$.
  * [75] This arises from the fact that the projection of the $K$ and $K^{\prime}$ points on the edge only coincide (realizing the armchair boundary conditions) for very specific angles between the edges and the lattice axes [91].
  * [76] See supplementary material.
  * [77] This twist angle and superconducting gap are not realistic for TTG and are chosen only to facilitate the convergence of the computation.
  * Khalaf  _et al._ [2019] E. Khalaf, A. J. Kruchkov, G. Tarnopolsky, and A. Vishwanath, Magic angle hierarchy in twisted graphene multilayers, Physical Review B 100, 10.1103/physrevb.100.085109 (2019).
  * Gonzalez-Arraga  _et al._ [2017] L. A. Gonzalez-Arraga, J. L. Lado, F. Guinea, and P. San-Jose, Electrically controllable magnetism in twisted bilayer graphene, Phys. Rev. Lett. 119, 107201 (2017).
  * Lin and Tománek [2018] X. Lin and D. Tománek, Minimum model for the electronic structure of twisted bilayer graphene and related structures, Phys. Rev. B 98, 081410 (2018).
  * Vahedi  _et al._ [2021] J. Vahedi, R. Peters, A. Missaoui, A. Honecker, and G. T. de Laissardière, Magnetism of magic-angle twisted bilayer graphene, SciPost Physics 11, 10.21468/scipostphys.11.4.083 (2021).
  * Sainz-Cruz  _et al._ [2021] H. Sainz-Cruz, T. Cea, P. A. Pantaleón, and F. Guinea, High transmission in twisted bilayer graphene with angle disorder, Phys. Rev. B 104, 075144 (2021).
  * Guinea and Walet [2018] F. Guinea and N. R. Walet, Electrostatic effects, band distortions, and superconductivity in twisted graphene bilayers, Proceedings of the National Academy of Sciences 115, 13174 (2018).
  * Bistritzer and MacDonald [2011] R. Bistritzer and A. H. MacDonald, Moiré bands in twisted double-layer graphene, Proceedings of the National Academy of Sciences 108, 12233 (2011).
  * Moon  _et al._ [2014] P. Moon, Y.-W. Son, and M. Koshino, Optical absorption of twisted bilayer graphene with interlayer potential asymmetry, Phys. Rev. B 90, 155427 (2014).
  * Cai  _et al._ [2010] J. Cai, P. Ruffieux, R. Jaafar, M. Bieri, T. Braun, S. Blankenburg, M. Muoth, A. P. Seitsonen, M. Saleh, X. Feng, _et al._ , Atomically precise bottom-up fabrication of graphene nanoribbons, Nature 466, 470 (2010).
  * Talirz  _et al._ [2016] L. Talirz, P. Ruffieux, and R. Fasel, On-surface synthesis of atomically precise graphene nanoribbons, Advanced materials 28, 6222 (2016).
  * Ruffieux  _et al._ [2016] P. Ruffieux, S. Wang, B. Yang, C. Sánchez-Sánchez, J. Liu, T. Dienel, L. Talirz, P. Shinde, C. A. Pignedoli, D. Passerone, _et al._ , On-surface synthesis of graphene nanoribbons with zigzag edge topology, Nature 531, 489 (2016).
  * Kitao  _et al._ [2020] T. Kitao, M. W. MacLean, K. Nakata, M. Takayanagi, M. Nagaoka, and T. Uemura, Scalable and precise synthesis of armchair-edge graphene nanoribbon in metal–organic framework, Journal of the American Chemical Society 142, 5509 (2020).
  * Houtsma  _et al._ [2021] R. K. Houtsma, J. de la Rie, and M. Stöhr, Atomically precise graphene nanoribbons: interplay of structural and electronic properties, Chemical Society Reviews 50, 6541 (2021).
  * Akhmerov and Beenakker [2008] A. R. Akhmerov and C. W. J. Beenakker, Boundary conditions for dirac fermions on a terminated honeycomb lattice, Phys. Rev. B 77, 085423 (2008).
  * Chiu  _et al._ [2013] C.-K. Chiu, H. Yao, and S. Ryu, Classification of topological insulators and superconductors in the presence of reflection symmetry, Phys. Rev. B 88, 075142 (2013).
  * Abanin  _et al._ [2006] D. A. Abanin, P. A. Lee, and L. S. Levitov, Spin-filtered edge states and quantum hall effect in graphene, Phys. Rev. Lett. 96, 176803 (2006).

Supplementary Material

##  I Monolayer graphene toy model

###  I.1 Tight-binding mean-field Hamiltonian

We consider the following spinless Hamiltonian defined on a bipartite
honeycomb lattice

$\hat{\mathcal{H}}=-t_{0}\sum_{\langle
ij\rangle}\hat{c}_{i}^{\dagger}\hat{c}_{j}+\mu\sum_{i}\hat{c}_{i}^{\dagger}\hat{c}_{i}+\frac{\Delta}{4i}\sum_{\langle\langle
ij\rangle\rangle}(-1)^{3\theta_{ij}}\hat{c}_{i}^{\dagger}\hat{c}_{j}^{\dagger}-\frac{\Delta}{4i}\sum_{\langle\langle
ij\rangle\rangle}(-1)^{3\theta_{ij}}\hat{c}_{j}\hat{c}_{i},$ (S1)

where $t_{0}>0$ is the nearest-neighbor hopping constant, $\mu$ is the
chemical potential, $\Delta$ is the magnitude of the superconducting order
parameter, and $\theta_{ij}$ is the angle between
$\mathbf{r}_{j}-\mathbf{r}_{i}$ and the $x$-axis in units of $\pi$. We note
parenthetically that the particle non-conserving terms look similar to the
Haldane next-nearest-neighbor hopping terms, except that there is no
sublattice dependence in our case. The fermionic operators satisfy
$\\{\hat{c}_{i},\hat{c}_{j}^{\dagger}\\}=\delta_{ij}.$ To obtain the momentum-
space Hamiltonian, we define

$\hat{c}^{\dagger}_{\mathbf{k}}=\frac{1}{\sqrt{\mathcal{N}}}\sum_{\mathbf{r}_{i}}e^{i\mathbf{k}\cdot\mathbf{r}_{i}}\hat{c}_{i}^{\dagger},\quad\hat{c}_{i}^{\dagger}=\frac{1}{\sqrt{\mathcal{N}}}\sum_{\mathbf{k}}e^{-i\mathbf{k}\cdot\mathbf{r}_{i}}\hat{c}^{\dagger}_{\mathbf{k}},$
(S2)

to find

$\begin{split}\hat{\mathcal{H}}=&-t_{0}\sum_{\mathbf{k}}\left(h_{\mathbf{k}}\hat{c}^{\dagger}_{\mathbf{k},A}\hat{c}_{\mathbf{k},B}+h_{\mathbf{k}}^{*}\hat{c}^{\dagger}_{\mathbf{k},B}\hat{c}_{\mathbf{k},A}\right)+\mu\sum_{\mathbf{k}}\left(\hat{c}^{\dagger}_{\mathbf{k},A}\hat{c}_{\mathbf{k},A}+\hat{c}^{\dagger}_{\mathbf{k},B}\hat{c}_{\mathbf{k},B}\right)+\\\
&+\frac{\Delta}{2}\sum_{\mathbf{k}}g_{\mathbf{k}}\left(\hat{c}^{\dagger}_{\mathbf{k},A}\hat{c}^{\dagger}_{-\mathbf{k},A}+\hat{c}^{\dagger}_{\mathbf{k},B}\hat{c}^{\dagger}_{-\mathbf{k},B}\right)+\frac{\Delta}{2}\sum_{\mathbf{k}}g_{\mathbf{k}}\left(\hat{c}_{\mathbf{k},A}\hat{c}_{-\mathbf{k},A}+\hat{c}_{\mathbf{k},B}\hat{c}_{-\mathbf{k},B}\right),\end{split}$
(S3)

where
$h_{\mathbf{k}}=1+e^{i\mathbf{k}\cdot\mathbf{a}_{1}}+e^{i\mathbf{k}\cdot\mathbf{a}_{2}}$
and
$g_{\mathbf{k}}=\sin\left(\mathbf{k}\cdot\left(\mathbf{a}_{1}-\mathbf{a}_{2}\right)\right)+\sin\left(-\mathbf{k}\cdot\mathbf{a}_{1}\right)+\sin\left(\mathbf{k}\cdot\mathbf{a}_{2}\right).$
Now, using the commutation relation
$\\{\hat{c}_{\mathbf{k},\sigma},\hat{c}_{\mathbf{k}^{\prime},\sigma^{\prime}}^{\dagger}\\}=\delta_{\mathbf{k},\mathbf{k}^{\prime}}\delta_{\sigma,\sigma^{\prime}},$
we can rewrite this Hamiltonian as

$\begin{split}\hat{\mathcal{H}}=&-\frac{t_{0}}{2}\sum_{\mathbf{k}}\left(h_{\mathbf{k}}\hat{c}^{\dagger}_{\mathbf{k},A}\hat{c}_{\mathbf{k},B}+h^{*}_{\mathbf{k}}\hat{c}^{\dagger}_{\mathbf{k},B}\hat{c}_{\mathbf{k},A}\right)+\frac{\mu}{2}\sum_{\mathbf{k}}\left(\hat{c}^{\dagger}_{\mathbf{k},A}\hat{c}_{\mathbf{k},A}+\hat{c}^{\dagger}_{\mathbf{k},B}\hat{c}_{\mathbf{k},B}\right)+\\\
&+\frac{t_{0}}{2}\sum_{\mathbf{k}}\left(h_{-\mathbf{k}}\hat{c}_{-\mathbf{k},B}\hat{c}^{\dagger}_{-\mathbf{k},A}+h^{*}_{-\mathbf{k}}\hat{c}_{-\mathbf{k},A}\hat{c}^{\dagger}_{-\mathbf{k},B}\right)-\frac{\mu}{2}\sum_{\mathbf{k}}\left(\hat{c}_{-\mathbf{k},A}\hat{c}^{\dagger}_{-\mathbf{k},A}+\hat{c}_{-\mathbf{k},B}\hat{c}^{\dagger}_{-\mathbf{k},B}\right)+\\\
&+\frac{\Delta}{2}\sum_{\mathbf{k}}g_{\mathbf{k}}\left(\hat{c}^{\dagger}_{\mathbf{k},A}\hat{c}^{\dagger}_{-\mathbf{k},A}+\hat{c}^{\dagger}_{\mathbf{k},B}\hat{c}^{\dagger}_{-\mathbf{k},B}\right)+\frac{\Delta}{2}\sum_{\mathbf{k}}g_{\mathbf{k}}\left(\hat{c}_{\mathbf{k},A}\hat{c}_{-\mathbf{k},A}+\hat{c}_{\mathbf{k},B}\hat{c}_{-\mathbf{k},B}\right)+\text{constant.}\end{split}$
(S4)

We henceforth drop the inconsequential constant term. Putting this into a
convenient matrix form, we obtain, using $h_{\mathbf{k}}=h^{*}_{-\mathbf{k}},$

$\begin{split}\hat{\mathcal{H}}&=\frac{1}{2}\sum_{\mathbf{k}}\hat{\Psi}^{\dagger}(\mathbf{k})\mathcal{H}^{\text{BdG}}(\mathbf{k})\hat{\Psi}(\mathbf{k}),\\\
\mathcal{H}^{\text{BdG}}(\mathbf{k})&=\begin{pmatrix}\mu&-t_{0}h_{\mathbf{k}}&\Delta
g_{\mathbf{k}}&0\\\ -t_{0}h^{*}_{\mathbf{k}}&\mu&0&\Delta g_{\mathbf{k}}\\\
\Delta g_{\mathbf{k}}&0&-\mu&t_{0}h_{\mathbf{k}}\\\ 0&\Delta
g_{\mathbf{k}}&t_{0}h^{*}_{\mathbf{k}}&-\mu\end{pmatrix},\\\
\hat{\Psi}^{\dagger}(\mathbf{k})&=\begin{pmatrix}\hat{c}_{\mathbf{k},A}^{\dagger}&\hat{c}_{\mathbf{k},B}^{\dagger}&\hat{c}_{-\mathbf{k},A}&\hat{c}_{-\mathbf{k},B}\end{pmatrix}.\end{split}$
(S5)

Figure S1: (a) Lattice structure, (b) phase diagram, and representative band
structures of the model Hamiltonian in Eq. (S1). Here, we use $t_{0}=1,$
$\Delta=0.2,$ and $\mu=0.1$ for the gapped phase.

Next, we study the symmetries of Hamiltonian (S1). Because of the
superconducting order parameter, this Hamiltonian breaks $\mathcal{C}_{2z}$
rotation symmetry and $\mathcal{M}_{x}$ mirror symmetry. However,
$\mathcal{C}_{3z}$ rotation symmetry and $\mathcal{M}_{y}$ mirror symmetry are
preserved. These are implemented by unitary operators:

$\begin{split}\mathcal{C}_{3z}:&\quad\hat{\mathcal{C}}_{3z}\begin{pmatrix}\hat{c}^{\dagger}_{A}\left(\mathbf{r}_{i}\right)\\\
\hat{c}_{A}\left(\mathbf{r}_{i}\right)\end{pmatrix}\hat{\mathcal{C}}_{3z}^{-1}=\begin{pmatrix}\hat{c}^{\dagger}_{A}\left(\mathcal{R}_{3z}\mathbf{r}_{i}\right)\\\
\hat{c}_{A}\left(\mathcal{R}_{3z}\mathbf{r}_{i}\right)\end{pmatrix},\\\
&\quad\hat{\mathcal{C}}_{3z}\begin{pmatrix}\hat{c}^{\dagger}_{B}\left(\mathbf{r}_{i}\right)\\\
\hat{c}_{B}\left(\mathbf{r}_{i}\right)\end{pmatrix}\hat{\mathcal{C}}_{3z}^{-1}=\begin{pmatrix}\hat{c}^{\dagger}_{B}\left(\mathcal{R}_{3z}\mathbf{r}_{i}\right)\\\
\hat{c}_{B}\left(\mathcal{R}_{3z}\mathbf{r}_{i}\right)\end{pmatrix},\end{split}$
(S6)
$\begin{split}\mathcal{M}_{y}:&\quad\hat{\mathcal{M}}_{y}\begin{pmatrix}\hat{c}_{A}^{\dagger}(\mathbf{r}_{i})\\\
\hat{c}_{A}(\mathbf{r}_{i})\end{pmatrix}\hat{\mathcal{M}}_{y}^{-1}=\begin{pmatrix}\hat{c}_{B}^{\dagger}\left(\mathcal{M}_{y}\mathbf{r}_{i}\right)\\\
\hat{c}_{B}\left(\mathcal{M}_{y}\mathbf{r}_{i}\right)\end{pmatrix},\\\
&\quad\hat{\mathcal{M}}_{y}\begin{pmatrix}\hat{c}_{B}^{\dagger}(\mathbf{r}_{i})\\\
\hat{c}_{B}(\mathbf{r}_{i})\end{pmatrix}\hat{\mathcal{M}}_{y}^{-1}=\begin{pmatrix}\hat{c}_{A}^{\dagger}\left(\mathcal{M}_{y}\mathbf{r}_{i}\right)\\\
\hat{c}_{A}\left(\mathcal{M}_{y}\mathbf{r}_{i}\right)\end{pmatrix},\end{split}$
(S7)

where we have adopted an alternative convention to $\hat{c}_{i}^{\dagger}$ to
emphasize the spatial dependence of the creation operators so that
$\hat{c}_{s}^{\dagger}(\mathbf{r})$ creates a fermion on the $s$ sublattice
located at spatial point $\mathbf{r}.$ We note that while $\mathcal{C}_{3z}$
preserves sublattice, $\mathcal{M}_{y}$ exchanges sublattices. In momentum
space, the mirror symmetry is represented by

$\begin{split}\mathcal{M}_{y}:&\quad\left(\tau_{0}\otimes\sigma_{x}\right)\mathcal{H}^{\text{BdG}}(\mathbf{k})\left(\tau_{0}\otimes\sigma_{x}\right)=\mathcal{H}^{\text{BdG}}(\mathcal{M}_{y}\mathbf{k}).\end{split}$
(S8)

Here, the $\sigma$ Pauli matrices act on the sublattice degree of freedom and
the $\tau$ Pauli matrices act on the Nambu space. One can easily verify that
the Hamiltonian does indeed satisfy these constraints. We note that
$\mathcal{M}_{y}^{2}=+1.$ In addition to these spatial symmetries, there is
also a local antiunitary particle-hole symmetry. It is most easily seen in
first-quantized form in momentum space where it takes the form

$\mathcal{P}:\quad\left(\tau_{x}\otimes\sigma_{0}\right)\left[\mathcal{H}^{\text{BdG}}(\mathbf{k})\right]^{*}\left(\tau_{x}\otimes\sigma_{0}\right)=-\mathcal{H}^{\text{BdG}}(-\mathbf{k}).$
(S9)

Additionally, we observe that there is spinless time-reversal symmetry
$\mathcal{T}$ represented by

$\mathcal{T}:\quad\left(\tau_{z}\otimes\sigma_{0}\right)\left[\mathcal{H}^{\text{BdG}}(\mathbf{k})\right]^{*}\left(\tau_{z}\otimes\sigma_{0}\right)=\mathcal{H}^{\text{BdG}}(-\mathbf{k}).$
(S10)

because $-g_{\mathbf{k}}=g_{-\mathbf{k}}$ but
$g_{\mathbf{k}}^{*}=g_{\mathbf{k}},$ necessitating an extra minus sign in the
off-diagonal elements in Nambu space. We note that $\mathcal{T}^{2}=+1.$ Of
course, the Hamiltonian also respects chiral symmetry
$\mathcal{C}=\mathcal{T}\mathcal{P}$ since $\mathcal{P}$ and $\mathcal{T}$ are
individually respected

$\mathcal{C}:\quad\left(i\tau_{y}\otimes\sigma_{0}\right)\mathcal{H}^{\text{BdG}}(\mathbf{k})\left(-i\tau_{y}\otimes\sigma_{0}\right)=-\mathcal{H}^{\text{BdG}}(\mathbf{k}).$
(S11)

We now study the band structure of Eq. (S1). The four energy bands of the BdG
matrix are given by

$E(\mathbf{k})=\pm\sqrt{\Delta^{2}g_{\mathbf{k}}^{2}+\left(\mu\pm
t_{0}|h_{\mathbf{k}}|\right)^{2}}.$ (S12)

Here,
$g_{\mathbf{k}}=2\sin\left(k_{x}a/2\right)\left[\cos\left(k_{x}a/2\right)-\cos\left(\sqrt{3}k_{y}a/2\right)\right],$
which vanishes along the lines $k_{x}=0,k_{x}=\pm\sqrt{3}k_{y}$ in the first
Brillouin zone. Because of $C_{3z}$ rotation symmetry, we focus on the
$k_{x}=0$ line. The other ones are related by symmetry. Along this line,
$|h_{k_{y}}|=\sqrt{5+4\cos(\sqrt{3}k_{y}a/2)}.$ For
$k_{y}\in\left[-2\pi/\sqrt{3}a,2\pi/\sqrt{3}a\right],$ the range of
$\cos(\sqrt{3}k_{y}a/2)$ is $[-1,1].$ So $|h_{k_{y}}|$ varies continuously in
$[1,3]$ as $k_{y}$ is sampled across the Brillouin zone. Therefore, always
assuming $\Delta\neq 0,$ we have a gap closing at $E=0$ if

$t_{0}\leq\left|\mu\right|\leq 3t_{0},$ (S13)

and we have a gapped phase if

$\left|\mu\right|<t_{0}\text{ or }\left|\mu\right|>3t_{0}.$ (S14)

Conventionally, the chemical potential will be small; so we will almost always
be in a gapped phase. To access the gapless phase, we have to dope at least to
$\mu=\pm t_{0}.$ For $|\mu|>3t_{0},$ we technically have another gapped phase
at $E=0,$ but this is kind of a trivial observation because it corresponds to
fully empty or fully filled fermionic bands. So we will not consider these
large values of $\mu$ further since they are irrelevant to our purpose. A
representative band structure is shown in Fig. S1(b).

###  I.2 Symmetry classification using the winding number

Now that the spectral properties of the band structure in parameter space are
well understood, we can classify the topology of the bulk bands below $E=0$ in
the case that the spectrum is gapped. First, since we have both
$\mathcal{T}^{2}=+1$ and $\mathcal{P}^{2}=+1,$ the system is in the two-
dimensional symmetry class $BDI,$ which is topologically trivial according to
the ten-fold classification [92]. However, we can enrich the classification by
considering mirror symmetry. Due to the presence of $\mathcal{M}_{y}$
symmetry, we can label states by mirror eigenvalues $m_{y}=\pm 1$ along
mirror-symmetric momentum lines [67, 68, 69, 92]. Notice that the mirror
eigenvalues are real because we are dealing with a spinless representation.
Along each of these lines, the system is one-dimensional and admits a
classification via the Berry phase within each mirror sector.

In our case, there is only one independent mirror-invariant line in the
Brillouin zone along $k_{y}=0.$ The line at $k_{y}=\frac{2\pi}{\sqrt{3}a}$ is
connected to the line at $k_{y}=0$ via a reciprocal lattice vector so it
cannot be considered an independent line. In fact, in order to traverse a
closed loop in $\mathbf{k}$-space along a mirror-invariant line with
$k_{y}=0,$ we need to go from $k_{x}=0$ to $k_{x}=4\pi/a,$ which goes through
both the line at $k_{y}=0$ and the line at $k_{y}=\frac{2\pi}{\sqrt{3}a},$
showing that they are indeed on the same line. This is also clear from
observation of the explicit form of the Hamiltonian

$\begin{split}\mathcal{H}^{\text{BdG}}(k_{x})&=\left(\begin{array}[]{cccc}\mu&h(k_{x})&g(k_{x})&0\\\
h(k_{x})&\mu&0&g(k_{x})\\\ g(k_{x})&0&-\mu&-h(k_{x})\\\
0&g(k_{x})&-h(k_{x})&-\mu\\\ \end{array}\right),\\\
h(k_{x})&=-t_{0}\left(2\cos\left(\frac{k_{x}a}{2}\right)+1\right),\\\
g(k_{x})&=\Delta\left(\sin(k_{x}a)-2\sin\left(\frac{k_{x}a}{2}\right)\right),\end{split}$
(S15)

which has period $4\pi/a.$ Along this line, we have

$\begin{split}\mathcal{M}_{y}\mathcal{H}^{\text{BdG}}(k_{x},k_{y}=0)\mathcal{M}_{y}^{-1}&=\mathcal{H}^{\text{BdG}}(k_{x},k_{y}=0).\end{split}$
(S16)

So we can partition the Hamiltonian into a mirror-even and mirror-odd block
Hamiltonian

$\mathcal{H}^{\text{BdG}}(k_{x})=\begin{pmatrix}\mathcal{H}_{-}(k_{x})&0\\\
0&\mathcal{H}_{+}(k_{x})\end{pmatrix},$ (S17)

where

$\begin{split}\mathcal{H}_{-}(k_{x})&=\left(\begin{array}[]{cc}h(k_{x})-\mu&g(k_{x})\\\
g(k_{x})&-h(k_{x})+\mu\\\ \end{array}\right),\\\
\mathcal{H}_{+}(k_{x})&=\left(\begin{array}[]{cc}-h(k_{x})-\mu&g(k_{x})\\\
g(k_{x})&h(k_{x})+\mu\\\ \end{array}\right).\end{split}$ (S18)

The mirror basis consists of

$\begin{split}m_{y}=\pm
1:\quad\ket{1_{\pm}}=\frac{1}{\sqrt{2}}\begin{pmatrix}0\\\ 0\\\ \pm 1\\\
1\end{pmatrix},\quad\ket{2_{\pm}}=\frac{1}{\sqrt{2}}\begin{pmatrix}\pm 1\\\
1\\\ 0\\\ 0\end{pmatrix}.\end{split}$ (S19)

The subscript refers to the mirror sector (i.e. the mirror eigenvalue of that
state), and the numbers $1,2$ just label the states within each mirror sector.
Now because $[\mathcal{M}_{y},\mathcal{C}]=0,$ we use $\mathcal{C}$ also to
simultaneously put the mirror-diagonal Hamiltonian into chiral form. In the
mirror basis, the chiral operator is

$\begin{split}\mathcal{C}=\begin{pmatrix}\mathcal{C}_{-}&0\\\
0&\mathcal{C}_{+}\end{pmatrix},\quad\mathcal{C}_{-}=\mathcal{C}_{+}=\left(\begin{array}[]{cc}0&-1\\\
1&0\\\ \end{array}\right).\end{split}$ (S20)

For instance, in the following basis within the mirror-odd sector

$m_{y}=-1:\quad\ket{-,+i}=\frac{1}{\sqrt{2}}\begin{pmatrix}i\\\
1\end{pmatrix},\quad\ket{-,-i}=\frac{1}{\sqrt{2}}\begin{pmatrix}-i\\\
1\end{pmatrix},$ (S21)

the mirror-odd Hamiltonian can be put into off-diagonal chiral form

$\begin{split}\mathcal{H}_{-}(k_{x})&=\begin{pmatrix}0&\mathcal{D}_{-}(k_{x})\\\
\mathcal{D}^{\dagger}_{-}(k_{x})&0\end{pmatrix},\\\
\mathcal{D}_{-}(k_{x})&=-2i\Delta\sin\left(\frac{k_{x}a}{2}\right)+i\Delta\sin(k_{x}a)+2t_{0}\cos\left(\frac{k_{x}a}{2}\right)+\mu+t_{0}.\end{split}$
(S22)

In a similar manner, we can also rewrite the mirror-even Hamiltonian in the
following basis

$m_{y}=+1:\quad\ket{+,+i}=\frac{1}{\sqrt{2}}\begin{pmatrix}i\\\
1\end{pmatrix},\quad\ket{+,-i}=\frac{1}{\sqrt{2}}\begin{pmatrix}-i\\\
1\end{pmatrix},$ (S23)

so that it too is in chiral form

$\begin{split}\mathcal{H}_{+}(k_{x})&=\begin{pmatrix}0&\mathcal{D}_{+}(k_{x})\\\
\mathcal{D}^{\dagger}_{+}(k_{x})&0\end{pmatrix},\\\
\mathcal{D}_{+}(k_{x})&=-2i\Delta\sin\left(\frac{k_{x}a}{2}\right)+i\Delta\sin(k_{x}a)-2t_{0}\cos\left(\frac{k_{x}a}{2}\right)+\mu-
t_{0}.\end{split}$ (S24)

Next, we “flatten” the Hamiltonian so that all the positive energy eigenvalues
are $+1$ and all the negative energy eigenvalues are $-1.$ This is
accomplished using the singular value decomposition [72]

$\mathcal{D}_{\pm}(k_{x})=\mathcal{U}_{\pm}(k_{x})\Lambda_{\pm}(k_{x})\mathcal{V}^{\dagger}_{\pm}(k_{x})$
(S25)

The topological invariant is the winding number

$\begin{split}\nu_{\pm}(k_{y}=0)&=\frac{1}{2\pi
i}\int_{0}^{4\pi/a}dk_{x}\text{Tr}\left(\tilde{\mathcal{D}}^{\dagger}_{\pm}(k_{x})\partial_{k_{x}}\tilde{\mathcal{D}}_{\pm}(k_{x})\right),\\\
\tilde{\mathcal{D}}_{\pm}(k_{x})&=\mathcal{U}_{\pm}(k_{x})\mathcal{V}^{\dagger}_{\pm}(k_{x}).\end{split}$
(S26)

Putting this together, we find

$\nu_{-}(k_{y}=0)=-\text{sign}\left(\frac{\Delta}{t_{0}}\right)=-\nu_{+}(k_{y}=0).$
(S27)

This means that when a mirror-respecting edge is made, at $k_{y}=0,$ there is
exactly one edge mode belonging to the mirror-odd sector and one edge mode
belonging to the mirror-even sector. Then, by chiral symmetry, these modes
must be exactly at zero energy since if there was a mode at positive energy in
one mirror sector, there must be another one at negative energy in the same
mirror sector since chiral symmetry commutes with mirror symmetry.

###  I.3 Dirac continuum description of armchair edge states

The edge states predicted by the topological classification can be obtained
directly in the continuum limit. Let us assume that $|\Delta|,|\mu|\ll|t_{0}|$
so that the relevant physics lies at the Dirac cones of monolayer graphene.
Then, we can expand the BdG Hamiltonian around $K$ and $K^{\prime}$ to find

$\begin{split}\mathcal{H}^{\text{BdG}}(\mathbf{k})&=\begin{pmatrix}\mathcal{H}_{K}(\mathbf{k})&0\\\
0&\mathcal{H}_{K^{\prime}}(\mathbf{k})\end{pmatrix},\\\
\mathcal{H}_{K}(\mathbf{k})&=\begin{pmatrix}\mu&\hbar
v_{F}\left(k_{x}-ik_{y}\right)&-\frac{3\sqrt{3}}{2}\Delta&0\\\ \hbar
v_{F}\left(k_{x}+ik_{y}\right)&\mu&0&-\frac{3\sqrt{3}}{2}\Delta\\\
-\frac{3\sqrt{3}}{2}\Delta&0&-\mu&-\hbar v_{F}\left(k_{x}-ik_{y}\right)\\\
0&-\frac{3\sqrt{3}}{2}\Delta&-\hbar
v_{F}\left(k_{x}+ik_{y}\right)&-\mu\end{pmatrix},\\\
\mathcal{H}_{K^{\prime}}(\mathbf{k})&=\begin{pmatrix}\mu&\hbar
v_{F}\left(-k_{x}-ik_{y}\right)&+\frac{3\sqrt{3}}{2}\Delta&0\\\ \hbar
v_{F}\left(-k_{x}+ik_{y}\right)&\mu&0&+\frac{3\sqrt{3}}{2}\Delta\\\
+\frac{3\sqrt{3}}{2}\Delta&0&-\mu&-\hbar v_{F}\left(-k_{x}-ik_{y}\right)\\\
0&+\frac{3\sqrt{3}}{2}\Delta&-\hbar
v_{F}\left(-k_{x}+ik_{y}\right)&-\mu\end{pmatrix}.\end{split}$ (S28)

We are interested in solving the wavefunction in real space where
$k_{\mu}=-i\partial_{\mu}.$ The basis state is labeled as

$\psi(\mathbf{r})=\begin{pmatrix}f_{A,K}(\mathbf{r})&f_{B,K}(\mathbf{r})&g_{A,K}(\mathbf{r})&g_{B,K}(\mathbf{r})&f_{A,K^{\prime}}(\mathbf{r})&f_{B,K^{\prime}}(\mathbf{r})&g_{A,K^{\prime}}(\mathbf{r})&g_{B,K^{\prime}}(\mathbf{r})\end{pmatrix}^{T}.$
(S29)

We are interested in solving the wavefunction at an armchair edge which we
locate at $x=0$ and the sample extending in the $x\rightarrow-\infty$
direction. We take the boundary condition to be

$f_{\sigma,K}(x=0)=f_{\sigma,K^{\prime}}(x=0),\quad
g_{\sigma,K}(x=0)=g_{\sigma,K^{\prime}}(x=0).$ (S30)

Actually, the correct physically-derived boundary condition that requires both
sublattices to vanish at the interface should be

$f_{\sigma,K}(x=0)=-f_{\sigma,K^{\prime}}(x=0),\quad
g_{\sigma,K}(x=0)=-g_{\sigma,K^{\prime}}(x=0).$ (S31)

However, Eq. (S30) and Eq. (S31) are related to each other by a simple
transformation of the basis functions. So we will continue to use Eq. (S30).
If necessary, one can multiply the spinor in the $K^{\prime}$ valley by $-1.$
With the boundary condtion of Eq. (S30), instead of solving the wavefunction
in the half-plane, we can extend the domain to the entire plane where the
wavefunction on the left $(x<0)$ belongs to the $K$ valley while the
wavefunction on the right belongs to the $K^{\prime}$ valley. This amounts to
mapping $x\mapsto-x$ on the right half-plane. The boundary condition is
satisfied automatically if the wavefunction is continuous at $x=0.$ This trick
was used in Ref. [93]. However, $\Delta$ is now spatially dependent

$\zeta(x)=\begin{cases}-\frac{3\sqrt{3}}{2}\Delta&x\leq 0\\\
+\frac{3\sqrt{3}}{2}\Delta&x>0\end{cases}.$ (S32)

First, let us focus on the zero-energy modes that occur at $k_{y}=0.$ Then
going into the mirror basis, we obtain

$\mathcal{H}(\mathbf{r})=\left(\begin{array}[]{cccc}-\mu-i\hbar
v_{F}\partial_{x}&\zeta(x)&0&0\\\ \zeta(x)&\mu+i\hbar v_{F}\partial_{x}&0&0\\\
0&0&-\mu+i\hbar v_{F}\partial_{x}&\zeta(x)\\\ 0&0&\zeta(x)&\mu-i\hbar
v_{F}\partial_{x}\\\ \end{array}\right)$ (S33)

The upper block belongs to the mirror-odd sector, and the bottom block belongs
to the mirror-even sector. The problem now is essentially just solving a
$2\times 2$ Hamiltonian. There are two solutions: one belonging to each mirror
sector. If $\Delta>0,$ then the normalizable, but non-normalized, solutions
are

$\psi_{m_{y}=-1}(\mathbf{r})=\exp\left(\frac{-\zeta(x)+i\mu}{\hbar
v_{F}}x\right)\begin{pmatrix}i\\\ 1\\\ 0\\\
0\end{pmatrix},\quad\psi_{m_{y}=+1}(\mathbf{r})=\exp\left(\frac{-\zeta(x)-i\mu}{\hbar
v_{F}}x\right)\begin{pmatrix}0\\\ 0\\\ -i\\\ 1\end{pmatrix}.$ (S34)

If $\Delta<0,$ then the normalizable, but non-normalized, solutions are

$\psi_{m_{y}=-1}(\mathbf{r})=\exp\left(\frac{\zeta(x)+i\mu}{\hbar
v_{F}}x\right)\begin{pmatrix}-i\\\ 1\\\ 0\\\
0\end{pmatrix},\quad\psi_{m_{y}=+1}(\mathbf{r})=\exp\left(\frac{\zeta(x)-i\mu}{\hbar
v_{F}}x\right)\begin{pmatrix}0\\\ 0\\\ i\\\ 1\end{pmatrix}.$ (S35)

We can actually find the states away from zero energy which are connected to
these zero-energy states. They are

$\begin{split}\mathcal{E}(k_{y})=\pm\hbar
v_{F}k_{y}\frac{1}{1+\left(\frac{2\mu}{3\sqrt{3}\Delta}\right)^{2}}:\quad\psi_{\mathcal{E}=\pm}(\mathbf{r})=e^{ik_{y}y}\left[\pm
i\psi_{m_{y}=-1}(\mathbf{r})+\psi_{m_{y}=+1}(\mathbf{r})\right].\end{split}$
(S36)

Next, we consider how to spatially isolate a single Majorana mode. To do this,
we add a perturbation along the $y$ direction that breaks mirror symmetry. One
such perturbation is a $\sigma_{z}$ mass (in the original sublattice basis)
that changes sign

$\lim_{y\rightarrow\pm\infty}m(y)=\pm m,$ (S37)

where $2m$ denotes the energy difference between the $A$ and $B$ sublattices.
Strictly speaking, we should only add the perturbation at the edge because we
want to preserve the topological classification of the bulk. In this case, we
should multiply the mass term by a Dirac delta function $\delta(x)$ to limit
the spatial extent of the perturbation. However, it is analytically simpler to
let the perturbation be uniform in $x.$ Since we only care about the boundary
states, the effect of the mass deep inside the bulk should be minimal on the
edge states. The Hamiltonian is modified to

$\mathcal{H}(\mathbf{r})=\left(\begin{array}[]{cccc}-i\hbar
v_{F}\partial_{x}-\mu&\zeta(x)&-\hbar v_{F}\partial_{y}+m(y)&0\\\
\zeta(x)&\mu+i\hbar v_{F}\partial_{x}&0&\hbar v_{F}\partial_{y}-m(y)\\\ +\hbar
v_{F}\partial_{y}+m(y)&0&+i\hbar v_{F}\partial_{x}-\mu&\zeta(x)\\\ 0&-\hbar
v_{F}\partial_{y}-m(y)&\zeta(x)&-i\hbar v_{F}\partial_{x}+\mu\\\
\end{array}\right).$ (S38)

The two possible zero-energy solutions are (for $\Delta>0$)

$\begin{split}\psi_{1}(\mathbf{r})&=\exp\left(\frac{-\zeta(x)+i\mu}{\hbar
v_{F}}x\right)\exp\left(-\frac{1}{\hbar
v_{F}}\int_{0}^{y}m(y^{\prime})dy^{\prime}\right)\begin{pmatrix}i\\\ 1\\\ 0\\\
0\end{pmatrix},\\\ \psi_{2}(\mathbf{r})&=\exp\left(\frac{-\zeta(x)-i\mu}{\hbar
v_{F}}x\right)\exp\left(+\frac{1}{\hbar
v_{F}}\int_{0}^{y}m(y^{\prime})dy^{\prime}\right)\begin{pmatrix}0\\\ 0\\\
-i\\\ 1\end{pmatrix}.\end{split}$ (S39)

Only one of these solutions is normalizable depending on whether $m$ is
positive or negative. The solutions for $\Delta<0$ are

$\begin{split}\psi_{1}(\mathbf{r})&=\exp\left(\frac{\zeta(x)+i\mu}{\hbar
v_{F}}x\right)\exp\left(-\frac{1}{\hbar
v_{F}}\int_{0}^{y}m(y^{\prime})dy^{\prime}\right)\begin{pmatrix}-i\\\ 1\\\
0\\\ 0\end{pmatrix},\\\
\psi_{2}(\mathbf{r})&=\exp\left(\frac{\zeta(x)-i\mu}{\hbar
v_{F}}x\right)\exp\left(+\frac{1}{\hbar
v_{F}}\int_{0}^{y}m(y^{\prime})dy^{\prime}\right)\begin{pmatrix}0\\\ 0\\\ i\\\
1\end{pmatrix}.\end{split}$ (S40)

###  I.4 Numerical simulation of finite flakes

Figure S2: Charge density of zero modes in rectangular flakes. For all
simulations, $\Delta=0.1t_{0},$ the lengths are approximately $20a,$ and the
widths are approximately $100a-110a.$ Each plot shows a zero-energy eigenstate
localized on the left (the complementary right-localized state does exist, but
is not shown). The mirror eigenvalue of the displayed state is shown on top of
each plot. We find that as the chemical potential switches sign, the mirror
eigenvalue of the zero modes also changes. Here, we only show the charge
density of the electron sector. The two upper plots feature zigzag top and
bottom edges, and the two lower plots feature bearded top and bottom edges.
All energies are approximately $10^{-16}t_{0}.$

In this section, we numerically diagonalize the Hamiltonian for finite flakes.
We first study rectangular flakes with armchair edges along the $y$-axis and
zigzag or bearded edges along the $x$-axis. We intentionally make the flakes
respect $\mathcal{M}_{y}$ symmetry so that we can label states by their mirror
eigenvalue. Each such flake features only two zero-energy states (these states
are not exactly at zero energy due to finite-size mixing effects, but they
have very small energies $\lesssim 10^{-15}$) instead of four zero-energy
states expected in the infinite nanoribbon case. This difference in the number
of modes is probably due to the imposition of nontrivial boundary conditions
at bearded or zigzag edges. Depending on whether the top and bottom edges are
zigzag or bearded and whether the chemical potential is positive or negative,
the edge states have either positive or negative mirror eigenvalue, as shown
in Fig. S2. In other words, the presence of zigzag or bearded edges selects
out one of the two zero modes per edge with definite mirror character, and
completely suppresses the other zero mode. This behavior is very similar to
the effect of a $\sigma_{z}$ mass. Therefore, under appropriate conditions,
one can heuristically think of the imposition of zigzag and bearded edges as
having a varying $\sigma_{z}$ mass that is arbitrarily large (with the
appropriate signs) outside the sample and zero inside the sample. We have only
inferred this connection, and have not rigorously shown it to be true, but it
seems to explain the number of modes as well as their mirror character
obtained numerically.

To separate the two zero modes spatially on the left and the right of the
sample, we need to make the width along the $x$-direction large enough to
prevent wavefunction overlap. However, the length along the $y$-direction does
not qualitatively affect the presence of the zero modes. In fact, we can make
the length microscopic, and the zero modes would still be there. This is true
for both zigzag and bearded edges, as shown in Fig. S3. We also simulate a
regular hexagonal flake. In this case, we find six zero modes, each localized
at a corner of the flake. An example is shown in Fig. S4. The same conclusion
holds when the hexagonal flake is not regular, as shown in Fig. S5. All of
these finite-size calculations show that the present platform offers a myriad
of ways to manipulate the Majorana zero modes by design of appropriate edges
and corners.

Figure S3: Zero modes as a function of flake size. For all simulations,
$\Delta=0.1t_{0},$ $\mu=0.2t_{0},$ and the widths are approximately $160a.$ As
we vary the length of the flakes, we observe that the zero modes remain and
are localized to the edge. The upper plots feature zigzag edges at the top and
bottom, while the lower plots feature bearded edges at the top and bottom. All
energies are approximately $10^{-15}t_{0}.$ Figure S4: Corner zero mode in a
hexagonal flake. A zero mode is trapped at a corner of a hexagonal flake.
There are six such zero modes in total. Here, $\mu=0.2t_{0}$ and
$\Delta=0.1t_{0}.$ The energy is about $10^{-12}t_{0}.$ Figure S5: Corner
modes in a distorted hexagonal flake. The flake is not a regular hexagon, but
instead have non-congruent sides. The energies are $\lesssim 10^{-9}t_{0}.$
There are six near zero-modes located at the corners of the flake. Here,
$\mu=0.2t_{0}$ and $\Delta=0.1t_{0}.$

##  II Rhombohedral multilayer graphene

In this section, we generalize the previous results to multilayer systems.

###  II.1 Bernal bilayer graphene

The Hamiltonian for Bernal bilayer graphene is just that of two copies of
monolayer graphene with an interlayer coupling $\gamma_{1}=0.1t_{0}$

$\mathcal{H}^{\text{BdG}}(\mathbf{k})=\begin{pmatrix}\mu&-t_{0}h_{\mathbf{k}}&0&0&\Delta
g_{\mathbf{k}}&0&0&0\\\ -t_{0}h^{*}_{\mathbf{k}}&\mu&\gamma_{1}&0&0&\Delta
g_{\mathbf{k}}&0&0\\\ 0&\gamma_{1}&\mu&-t_{0}h_{\mathbf{k}}&0&0&\Delta
g_{\mathbf{k}}&0\\\ 0&0&-t_{0}h^{*}_{\mathbf{k}}&\mu&0&0&0&\Delta
g_{\mathbf{k}}\\\ \Delta g_{\mathbf{k}}&0&0&0&-\mu&t_{0}h_{\mathbf{k}}&0&0\\\
0&\Delta g_{\mathbf{k}}&0&0&t_{0}h^{*}_{\mathbf{k}}&-\mu&-\gamma_{1}&0\\\
0&0&\Delta g_{\mathbf{k}}&0&0&-\gamma_{1}&-\mu&t_{0}h_{\mathbf{k}}\\\
0&0&0&\Delta g_{\mathbf{k}}&0&0&t_{0}h^{*}_{\mathbf{k}}&-\mu\end{pmatrix}.$
(S41)

Particle-hole, time reversal, and chiral symmetries have the same
representation as before

$\begin{split}\mathcal{P}:&\quad\left(\tau_{x}\otimes\lambda_{0}\otimes\sigma_{0}\right)\left[\mathcal{H}^{\text{BdG}}(\mathbf{k})\right]^{*}\left(\tau_{x}\otimes\lambda_{0}\otimes\sigma_{0}\right)=-\mathcal{H}^{\text{BdG}}(-\mathbf{k}),\\\
\mathcal{T}:&\quad\left(\tau_{z}\otimes\lambda_{0}\otimes\sigma_{0}\right)\left[\mathcal{H}^{\text{BdG}}(\mathbf{k})\right]^{*}\left(\tau_{z}\otimes\lambda_{0}\otimes\sigma_{0}\right)=\mathcal{H}^{\text{BdG}}(-\mathbf{k}),\\\
\mathcal{C}:&\quad\left(i\tau_{y}\otimes\lambda_{0}\otimes\sigma_{0}\right)\mathcal{H}^{\text{BdG}}(\mathbf{k})\left(-i\tau_{y}\otimes\lambda_{0}\otimes\sigma_{0}\right)=-\mathcal{H}^{\text{BdG}}(\mathbf{k}),\end{split}$
(S42)

where the $\lambda$ Pauli matrices act on the layer degree of freedom.
However, mirror symmetry needs a slight modification. Actually, in the context
of Bernal bilayer graphene, mirror symmetry is a $\mathcal{C}_{2x}$ rotation
symmetry that exchanges the two layers as well as sublattices. Despite that,
we continue to use $\mathcal{M}_{y}$ when we mean $\mathcal{C}_{2x}.$ The
symmetry operator is

$\mathcal{M}_{y}:\quad\left(\tau_{0}\otimes\lambda_{x}\otimes\sigma_{x}\right)\mathcal{H}^{\text{BdG}}(\mathbf{k})\left(\tau_{0}\otimes\lambda_{x}\otimes\sigma_{x}\right)=\mathcal{H}^{\text{BdG}}(\mathcal{M}_{y}\mathbf{k}).$
(S43)

Along this line, we can again block-diagonalize the Hamiltonian

$\begin{split}\mathcal{H}^{\text{BdG}}(k_{x})&=\begin{pmatrix}\mathcal{H}_{-}(k_{x})&0\\\
0&\mathcal{H}_{+}(k_{x})\end{pmatrix},\end{split}$ (S44)

where

$\begin{split}\mathcal{H}_{-}(k_{x})&=\left(\begin{array}[]{cccc}-\mu&-h(k_{x})&g(k_{x})&0\\\
-h(k_{x})&\gamma_{1}-\mu&0&g(k_{x})\\\ g(k_{x})&0&\mu&h(k_{x})\\\
0&g(k_{x})&h(k_{x})&\mu-\gamma_{1}\\\ \end{array}\right),\\\
\mathcal{H}_{+}(k_{x})&=\left(\begin{array}[]{cccc}-\mu&-h(k_{x})&g(k_{x})&0\\\
-h(k_{x})&-\gamma_{1}-\mu&0&g(k_{x})\\\ g(k_{x})&0&\mu&h(k_{x})\\\
0&g(k_{x})&h(k_{x})&\gamma_{1}+\mu\\\ \end{array}\right).\end{split}$ (S45)

Putting into chiral form, we obtain

$\begin{split}\mathcal{H}_{\pm}(k_{x})&=\begin{pmatrix}0&\mathcal{D}_{\pm}(k_{x})\\\
\mathcal{D}^{\dagger}_{\pm}(k_{x})&0\end{pmatrix},\\\
\mathcal{D}_{\pm}(k_{x})&=\left(\begin{array}[]{cc}\pm\gamma_{1}+ig(k_{x})+\mu&h(k_{x})\\\
h(k_{x})&\mu+ig(k_{x})\\\ \end{array}\right).\end{split}$ (S46)

Calculating the winding number, we find that $\nu_{\pm}=0$ for much of
parameter space. So, this system is topologically trivial. This is confirmed
by band structure calculation on an armchair nanoribbon that respects
$\mathcal{M}_{y}$ symmetry, shown in Fig. S6(a). There, we see no zero modes
at $k_{y}=0.$ However, interestingly, we do see zero modes displaced away from
$k_{y}=0$ for the present model. These are not topological though, as is clear
from the fact that they are connected to the same bulk bands. To confirm this,
we add a layer-antisymmetric $\sigma_{z}$ mass of the form

$V=m\tau_{z}\otimes\lambda_{z}\otimes\sigma_{z}.$ (S47)

This potential energy respects $\mathcal{M}_{y}$ because it simultaneously
flips sign under layer and sublattice exchange. So it commutes with
$\mathcal{M}_{y}.$ By increasing $m,$ we can remove the zero modes at
$k_{y}\neq 0,$ as shown in Fig. S6(b)-(d), without closing the bulk gap,
illustrating that these modes are accidental.

Figure S6: Band structure of superconducting Bernal bilayer graphene. (a)
There are no zero modes at $k_{y}=0,$ but there are zero modes away from
$k_{y}=0.$ (b)-(d) By adding a mass term of the form Eq. (S47), these non-
topological zero modes can be eliminated without breaking $\mathcal{M}_{y}$ or
closing the bulk gap. Here, $\mu=\Delta=0.05t_{0}.$

###  II.2 Rhombohedral trilayer graphene

The Hamiltonian for superconducting rhombohedral trilayer graphene is

$\mathcal{H}^{\text{BdG}}(\mathbf{k})=\begin{pmatrix}\mathcal{H}_{3}(\mathbf{k})&\mathcal{H}_{\text{sc}}(\mathbf{k})\\\
\mathcal{H}_{\text{sc}}^{\dagger}(\mathbf{k})&-\mathcal{H}_{3}(\mathbf{k})\end{pmatrix},$
(S48)

where

$\begin{split}\mathcal{H}_{3}(\mathbf{k})&=\begin{pmatrix}1&0&0\\\ 0&1&0\\\
0&0&1\end{pmatrix}\otimes\mathcal{H}_{\text{intra}}(\mathbf{k})+\begin{pmatrix}0&1&0\\\
0&0&1\\\
0&0&0\end{pmatrix}\otimes\mathcal{H}_{\text{inter}}(\mathbf{k})+\begin{pmatrix}0&0&0\\\
1&0&0\\\
0&1&0\end{pmatrix}\otimes\mathcal{H}^{\dagger}_{\text{inter}}(\mathbf{k}),\\\
\mathcal{H}_{\text{intra}}(\mathbf{k})&=\begin{pmatrix}\mu&-t_{0}h_{\mathbf{k}}\\\
-t_{0}h^{*}_{\mathbf{k}}&\mu\end{pmatrix},\\\
\mathcal{H}_{\text{inter}}(\mathbf{k})&=\begin{pmatrix}0&0\\\
\gamma_{1}&0\end{pmatrix},\\\
\mathcal{H}_{\text{sc}}(\mathbf{k})&=\begin{pmatrix}1&0&0\\\ 0&1&0\\\
0&0&1\end{pmatrix}\otimes\begin{pmatrix}\Delta g_{\mathbf{k}}&0\\\ 0&\Delta
g_{\mathbf{k}}\end{pmatrix}.\end{split}$ (S49)

The mirror operator is represented by

$\mathcal{M}_{y}=\tau_{0}\otimes\begin{pmatrix}0&0&1\\\ 0&1&0\\\
1&0&0\end{pmatrix}\otimes\sigma_{x},$ (S50)

which exchanges the top and bottom layers, but maps the middle layer to
itself. We also add a mirror-preserving potential

$V=m\tau_{z}\otimes\begin{pmatrix}1&0&0\\\ 0&0&0\\\
0&0&-1\end{pmatrix}\otimes\sigma_{z}.$ (S51)

The winding number can be calculated as before. We find that $\nu_{\pm}=\pm 1$
for this system. So this system is topologically nontrivial similar to the
monolayer system. We confirm this with finite calculations on an armchair
nanoribbon. We find two Majorana zero modes per boundary pinned to $k_{y}=0,$
as shown in Fig. S7(a). There are other zero modes at $k_{y}\neq 0,$ but these
are not topological. They can be eliminated by $\mathcal{M}_{y}$-symmetric
potentials without closing the bulk gap, as illustrated in Fig. S7(b)-(d).

Figure S7: Band structure of superconducting rhombohedral trilayer graphene.
(a) There are two zero modes per edge at $k_{y}=0;$ there are also zero modes
away from $k_{y}=0.$ (b)-(d) By adding a mass term of the form Eq. (S51),
these non-topological zero modes can be eliminated without breaking
$\mathcal{M}_{y}$ or closing the bulk gap. However, the topological zero modes
at $k_{y}=0$ are not gapped since they are robust. Here,
$\mu=\Delta=0.05t_{0}.$

In trilayer graphene and its multilayer generalization, the layer equivalence
can be broken simply by adding a perpendicular electric field. The
perturbation is described by

$V=V_{0}\tau_{z}\otimes\begin{pmatrix}1&0&0\\\ 0&0&0\\\
0&0&-1\end{pmatrix}\otimes\sigma_{0}.$ (S52)

This electric field hybridizes the topological zero modes and gap them out, as
shown in Fig. S8.

Figure S8: Electric-field-induced hybridization of zero modes in rhombohedral
trilayer graphene. With the addition of a perturbation in the form of Eq.
(S52), the zero modes at $k_{y}=0$ are gapped out because of broken mirror
symmetry. Here, $\mu=\Delta=0.1t_{0}$ and $V_{0}=0.3t_{0}.$

###  II.3 Rhombohedral $N$-layer generalization

The generalization to $N>1$ layers is straightforward. The Hamiltonian is

$\mathcal{H}^{\text{BdG}}(\mathbf{k})=\begin{pmatrix}\mathcal{H}_{N}(\mathbf{k})&\mathcal{H}_{\text{sc}}(\mathbf{k})\\\
\mathcal{H}_{\text{sc}}^{\dagger}(\mathbf{k})&-\mathcal{H}_{N}(\mathbf{k})\end{pmatrix},$
(S53)

where

$\begin{split}\mathcal{H}_{N}(\mathbf{k})&=\mathbb{I}_{N\times
N}\otimes\mathcal{H}_{\text{intra}}(\mathbf{k})+\mathbb{U}_{N\times
N}\otimes\mathcal{H}_{\text{inter}}(\mathbf{k})+\mathbb{L}_{N\times
N}\otimes\mathcal{H}^{\dagger}_{\text{inter}}(\mathbf{k}),\\\
\mathcal{H}_{\text{intra}}(\mathbf{k})&=\begin{pmatrix}\mu&-t_{0}h_{\mathbf{k}}\\\
-t_{0}h^{*}_{\mathbf{k}}&\mu\end{pmatrix},\\\
\mathcal{H}_{\text{inter}}(\mathbf{k})&=\begin{pmatrix}0&0\\\
\gamma_{1}&0\end{pmatrix},\\\
\mathcal{H}_{\text{sc}}(\mathbf{k})&=\mathbb{I}_{N\times
N}\otimes\begin{pmatrix}\Delta g_{\mathbf{k}}&0\\\ 0&\Delta
g_{\mathbf{k}}\end{pmatrix}.\end{split}$ (S54)

The mirror operator is represented by

$\mathcal{M}_{y}=\tau_{0}\otimes\bar{\mathbb{I}}_{N\times
N}\otimes\sigma_{x}.$ (S55)

Here, $\mathbb{I}_{N\times N}$ is the $N\times N$ identity matrix,
$\mathbb{U}_{N\times N}$ ($\mathbb{L}_{N\times N}$) is the sparse $N\times N$
matrix with ones along the upper (lower) diagonal right above (below) the
principal diagonal, and $\bar{\mathbb{I}}_{N\times N}$ is the anti-diagonal
matrix of ones. We note that for $N$ even, all layers are flipped under mirror
symmetry, while for $N$ odd, there is one central layer which is mapped to
itself under mirror symmetry. Carrying out the same topological classification
exactly as before, we find that the even-$N$ stacks are topologically trivial,
while the odd-$N$ stacks are topologically nontrivial protected by
$\mathcal{M}_{y}$ symmetry with winding numbers $\nu_{\pm}=\pm 1.$ We have
checked this fact numerically for $N\leq 50.$ Though we have not proven this
trend to be true for all $N,$ it seems likely that it is true for all values
of $N.$ The armchair band structures for a few values of $N$ are shown in Fig.
S9.

Figure S9: Band structure of armchair nanoribbons of multilayer stacks. We
observe that for stacks with an odd number of layers, there are zero modes at
$k_{y}=0,$ while for stacks with an even number of layers, there are not.

##  III Vortex-confined zero modes

In this section, we study zero modes confined to vortex cores numerically in
the tight-binding formalism. To do this, we only change the next-nearest-
neighbor superconducting hopping parameters. Without a vortex, the hopping
$\Delta(\mathbf{r}_{i},\mathbf{r}_{j})=\pm\Delta/2i$ is a constant for
$|\mathbf{r}_{i}-\mathbf{r}_{j}|=a$ and zero otherwise. With a vortex located
at $\mathbf{r}=\mathbf{0},$ we simply modify this hopping to include a spatial
dependence

$\Delta(\mathbf{r}_{i},\mathbf{r}_{j})=\pm\frac{\Delta}{2i}e^{\pm
i\phi_{\mathbf{r}}}\tanh\left(\frac{|\mathbf{r}|}{d}\right),$ (S56)

where $d$ is the distance over which the superconducting gap vanishes,
$\mathbf{r}=(\mathbf{r}_{i}+\mathbf{r}_{j})/2$ is the center-of-mass location,
and $\phi_{\mathbf{r}}$ is the angle of that position. As required, the
superconducting order parameter vanishes exactly at the vortex
$\mathbf{r}=\mathbf{0}.$ For monolayer graphene, we find that there are two
zero modes trapped at each vortex, as shown in Fig. S10, when the center of
the vortex is at the center of a hexagon or a mid-bond or another location
along a $\mathcal{M}_{y}$-symmetric line. When the position of the vortex
center is at a random location, then these low-energy modes are not exactly at
zero. However, when $d\gtrsim a,$ then the precise center of the vortex is not
as important, so the low-energy states are practically at zero energy for
various combinations of parameters. We have checked that the existence of the
vortex low-energy pairs is not sensitive to the location of the center
$\mathbf{r}=\mathbf{0}$ when $d\gtrsim a.$ In Fig. S10, we intentionally
choose shapes with irregular boundaries to show that the vortex states are not
sensitive to boundary conditions when the flakes are large.

Figure S10: Vortex states in monolayer graphene. Each flake has two modes near
zero energy trapped at the vortex, which is located at the center. The flakes
are intentionally irregularly shaped to illustrate that the vortex states are
not sensitive to boundary conditions. Here, $\mathbf{r}=\mathbf{0}$ is located
at the center of a hexagon. We use $\Delta=\mu=0.1t_{0}.$

##  IV $ABA$ multilayer graphene

In this section, we study $ABA$ multilayer graphene, which is the most stable
multilayer stacking order of graphite. Unfortunately, superconductivity has
not been experimentally observed in such systems. So the results in this
section might not be as experimentally relevant for the time being.
Nonetheless, from a theory point of view, these systems also host robust
Majorana zero modes due to mirror symmetries if they have $f$-wave
superconductivity.

###  IV.1 $ABA$ trilayer graphene

The Hamiltonian for $ABA$ trilayer graphene looks almost exactly like that for
rhombohedral trilayer graphene with one important exception in the interlayer
hopping

$\mathcal{H}^{\text{BdG}}(\mathbf{k})=\begin{pmatrix}\mathcal{H}_{3}(\mathbf{k})&\mathcal{H}_{\text{sc}}(\mathbf{k})\\\
\mathcal{H}_{\text{sc}}^{\dagger}(\mathbf{k})&-\mathcal{H}_{3}(\mathbf{k})\end{pmatrix},$
(S57)

where

$\begin{split}\mathcal{H}_{3}(\mathbf{k})&=\begin{pmatrix}1&0&0\\\ 0&1&0\\\
0&0&1\end{pmatrix}\otimes\mathcal{H}_{\text{intra}}(\mathbf{k})+\begin{pmatrix}0&1&0\\\
0&0&0\\\
0&1&0\end{pmatrix}\otimes\mathcal{H}_{\text{inter}}(\mathbf{k})+\begin{pmatrix}0&0&0\\\
1&0&1\\\
0&0&0\end{pmatrix}\otimes\mathcal{H}^{\dagger}_{\text{inter}}(\mathbf{k}),\\\
\mathcal{H}_{\text{intra}}(\mathbf{k})&=\begin{pmatrix}\mu&-t_{0}h_{\mathbf{k}}\\\
-t_{0}h^{*}_{\mathbf{k}}&\mu\end{pmatrix},\\\
\mathcal{H}_{\text{inter}}(\mathbf{k})&=\begin{pmatrix}0&0\\\
\gamma_{1}&0\end{pmatrix},\\\
\mathcal{H}_{\text{sc}}(\mathbf{k})&=\begin{pmatrix}1&0&0\\\ 0&1&0\\\
0&0&1\end{pmatrix}\otimes\begin{pmatrix}\Delta g_{\mathbf{k}}&0\\\ 0&\Delta
g_{\mathbf{k}}\end{pmatrix}.\end{split}$ (S58)

This Hamiltonian importantly respects a mirror symmetry in the $z$ direction
represented by

$\mathcal{M}_{z}=\tau_{0}\otimes\begin{pmatrix}0&0&1\\\ 0&1&0\\\
1&0&0\end{pmatrix}\otimes\sigma_{0}.$ (S59)

We have $\left[\mathcal{H}^{\text{BdG}}(\mathbf{k}),\mathcal{M}_{z}\right]=0.$
This symmetry is diagonal in $\mathbf{k}.$ So it holds for any $\mathbf{k}.$
Thus, we can block-diagonalize to write the Hamiltonian in the basis of
eigenstates of $\mathcal{M}_{z}.$ In this new basis, the Hamiltonian

$\mathcal{H}_{3}(\mathbf{k})=\begin{pmatrix}1&0&0\\\ 0&1&0\\\
0&0&1\end{pmatrix}\otimes\mathcal{H}_{\text{intra}}(\mathbf{k})+\begin{pmatrix}0&0&0\\\
0&0&\sqrt{2}\\\
0&0&0\end{pmatrix}\otimes\mathcal{H}_{\text{inter}}(\mathbf{k})+\begin{pmatrix}0&0&0\\\
0&0&0\\\
0&\sqrt{2}&0\end{pmatrix}\otimes\mathcal{H}^{\dagger}_{\text{inter}}(\mathbf{k}).$
(S60)

This decomposition is similar to that done in Ref. [78]. As one can see, in
this new basis, the system looks like one decoupled effective monolayer system
and one decoupled effective bilayer graphene. The monolayer
($\mathcal{M}_{z}$-odd) sector is invariant under $\mathcal{M}_{y}.$ We note
that the whole Hilbert space is not invariant under $\mathcal{M}_{y},$ only
the mirror-odd sector is. So we can again classify the monolayer bands by the
mirror-projected winding number as before. This predicts that everything we
have learned in the monolayer model still applies to the mirror-odd sector of
$ABA$ trilayer graphene. In particular, in armchair nanoribbons, we find a
pair of topological zero modes at $k_{y}=0$ as shown in Fig. S11(a). We also
find two zero modes confined to a vortex, as shown in Fig. S11(b). We note
that the zero mode is localized only in the top and bottom layers because it
comes from the monolayer sector, which in this case, is just an odd
combination of the top and bottom layers.

Figure S11: $ABA$ trilayer graphene.  (a) Band structure of an armchair
nanoribbon where the red bands are topological and come from the monolayer
sector of the Hamiltonian. (b) A zero mode confined to a vortex with energy
$E\approx 10^{-16}t_{0}$. The superconducting gap here is quite large
$\Delta=0.3t_{0}$ to facilitate fast computation.

###  IV.2 $ABA$ odd-layer graphene

We generalize the previous result to $ABA$ odd-layer graphene. The Hamiltonian
is

$\mathcal{H}^{\text{BdG}}(\mathbf{k})=\begin{pmatrix}\mathcal{H}_{N}(\mathbf{k})&\mathcal{H}_{\text{sc}}(\mathbf{k})\\\
\mathcal{H}_{\text{sc}}^{\dagger}(\mathbf{k})&-\mathcal{H}_{N}(\mathbf{k})\end{pmatrix},$
(S61)

where

$\begin{split}\mathcal{H}_{N}(\mathbf{k})&=\mathbb{I}_{N\times
N}\otimes\mathcal{H}_{\text{intra}}(\mathbf{k})+\begin{pmatrix}0&1&0&0&0&0&0&\dots\\\
0&0&0&0&0&0&0&\dots\\\ 0&1&0&1&0&0&0&\dots\\\ 0&0&0&0&0&0&0&\dots\\\
0&0&0&1&0&1&0&\dots\\\ 0&0&0&0&0&0&0&\dots\\\ 0&0&0&0&0&1&0&\dots\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\ddots\\\
\end{pmatrix}\otimes\mathcal{H}_{\text{inter}}(\mathbf{k})+\begin{pmatrix}0&0&0&0&0&0&0&\dots\\\
1&0&1&0&0&0&0&\dots\\\ 0&0&0&0&0&0&0&\dots\\\ 0&0&1&0&1&0&0&\dots\\\
0&0&0&0&0&0&0&\dots\\\ 0&0&0&0&1&0&1&\dots\\\ 0&0&0&0&0&0&0&\dots\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\ddots\\\
\end{pmatrix}\otimes\mathcal{H}^{\dagger}_{\text{inter}}(\mathbf{k}),\\\
\mathcal{H}_{\text{intra}}(\mathbf{k})&=\begin{pmatrix}\mu&-t_{0}h_{\mathbf{k}}\\\
-t_{0}h^{*}_{\mathbf{k}}&\mu\end{pmatrix},\\\
\mathcal{H}_{\text{inter}}(\mathbf{k})&=\begin{pmatrix}0&0\\\
\gamma_{1}&0\end{pmatrix},\\\
\mathcal{H}_{\text{sc}}(\mathbf{k})&=\mathbb{I}_{N\times
N}\otimes\begin{pmatrix}\Delta g_{\mathbf{k}}&0\\\ 0&\Delta
g_{\mathbf{k}}\end{pmatrix}.\end{split}$ (S62)

We assume that only nearest-layer hoppings and intralayer hoppings are
allowed. The rest of the discussion depends on this physically-motivated
assumption. This Hamiltonian is invariant under $\mathcal{M}_{z}$

$\mathcal{M}_{z}=\tau_{0}\otimes\bar{\mathbb{I}}_{N\times
N}\otimes\sigma_{0},$ (S63)

where $\bar{\mathbb{I}}_{N\times N}$ is the anti-diagonal $N\times N$ matrix.
There is one unique state

$\ket{\psi_{\text{mono}}}=\sqrt{\frac{2}{N+1}}\begin{pmatrix}1&0&-1&0&1&0&-1&0&...\end{pmatrix}^{T}$
(S64)

that does not couple to the interlayer hoppings. We refer to states in this
basis as coming from the effective monolayer sector because the Hamiltonian
governing these states looks exactly like the monolayer Hamiltonian. This
sector hosts topological zero modes at armchair edges, as shown in Fig. S12.
We see from there that these modes do not exist in even-layer configurations.

Figure S12: Band structure of $ABA$ multilayer graphene. We see that for $N$
odd, there are two modes crossing $k_{y}=0$ at exactly zero energy. Here,
$\mu=\Delta=0.1t_{0}.$

##  V Twisted multilayer graphene

###  V.1 Additional results for twisted bilayer graphene

Here, we present band structures for twisted bilayer graphene (TBG)
nanoribbons at six different twist angles and for both the hole and electron
superconducting domes. As shown in Fig. S13, the sub-gap spectra depend very
sensitively on twist angle. Starting at the magic angle (1.06∘ in the model),
increments as small as $0.1-0.2^{\circ}$ eliminate the zero modes and deplete
the energy window around the center of the gap. As pointed out in the main
text, the overall dynamic is the following: near the magic angle there are 8
zero modes; as the angle is increased, the momenta of these modes diminish;
increasing the angle further eliminates these modes; and at the largest angle
of 1.31∘, all Andreev states with energy $|\epsilon|<\Delta_{\text{sc}}/4$
have disappeared. This evolution is very similar in the hole and electron
domes (compare top and bottom rows in Fig. S13). However, there is electron-
hole asymmetry, which manifests quantitatively at every angle and
qualitatively in that the zero modes in the hole dome are more robust, in the
sense that they disappear at a larger twist angle.

Figure S13: Band structures of twisted bilayer graphene nanoribbons at
different twist angles. The top rows correspond to superconductivity in the
hole dome (the gap opens at filling $n=-2.4$), and the bottom rows correspond
to the electron dome (gap at filling $n=2.4$). In all cases,
$\Delta_{\text{sc}}=1$ meV and the the width of the ribbon is $W\approx
35L_{\text{M}}$.

###  V.2 Details of the model

The TBG nanoribbon unit cell is built following the procedure of Ref. [82]. In
the notation of that paper, we build a TBG nanoribbon with chiral vectors
$(440,20)@(-440,-20)$, which has a twist angle of 4.41∘, a width of $W\approx
35L_{\text{M}}$, and $27,040$ sites in its unit cell. We have checked that the
band structures have converged with respect to ribbon width. One detail to
take into account is that in the construction of Ref. [82] the origin of
coordinates is at an atom, while in the present case it is necessary to set
the origin at the center of a graphene hexagon, for the resulting system to
have $\mathcal{M}_{y}$ symmetry.

Figure S14: Unit cell of a TBG nanoribbon with chiral vector
$(22,1)@(-22,-1)$. Joining 20 such cells horizontally gives the unit cell for
the nanoribbon with chiral vector $(440,20)@(-440,-20)$, of width $W\approx
35L_{\text{M}}$, for which we obtain the spectra in Fig. 4(a)-(c) in the main
text and Fig. S13 above.

We start from the Lin-Tománek tight-binding, non-interacting Hamiltonian
${\cal H}_{\theta}$ [80] and, for TBG, we include Hartree electron-electron
interactions through an electrostatic potential [83]. Then, we build
$\mathcal{H}^{\text{BdG}}$ with the same superconducting hoppings as before:

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

where $i,j$ run over the lattice sites and $m$ is the layer index.
${\cal{H}}_{\theta}$ includes intralayer hopping to nearest-neighbors only
$\gamma_{ij}^{mm}=t_{0}$ 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_{0}=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 [80].

The parameters in the tight binding model are scaled, so that the central
bands of TBG with twist angle $\theta$ are approximated by the central bands
of an equivalent lattice with twist angle $\lambda\theta$, with $\lambda>1$
[79, 81, 82, 54]. This scaling approximation is based on the fact that the
Dirac equation that governs each layer can be described in different ways, and
one of them is a honeycomb lattice of super-atoms, each representing a small
cluster of atoms of the original lattice. In this scaled system, the lattice
constant is magnified and the hopping reduced. In the context of twisted
bilayer graphene, the possibility of scaling manifests also in the continuum
model, which shows that the bands depend, to first order, on a dimensionless
parameter [84],

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

where $a$ is the lattice constant and $v_{F}$ is the Fermi velocity. This
parameter can be understood as a comparison between the time a carrier needs
to traverse a unit cell within a layer, and the average time between
interlayer tunneling events. Thus, a small angle $\theta$ can be simulated
with a larger one $\theta^{\prime}$ by doing the following transformations:
$t_{0}\rightarrow\frac{1}{\lambda}t_{0}$, $a\rightarrow\lambda a$,
$d\rightarrow\lambda d$, with
$\lambda=\sin(\frac{\theta^{\prime}}{2})/\sin(\frac{\theta}{2})$ [79, 81, 82,
54]. Note that the interlayer distance is scaled to keep the interlayer
hopping unchanged after scaling. This approximation reproduces well the low-
energy band structure, as shown in Refs. [81, 54].

In Eq. (S65), the Hartree term is

$V_{\text{H}}(n)=\frac{2\rho(n)}{\epsilon_{r}L_{\text{M}}}\sum^{3}_{i=1}\cos(\mathbf{G}_{i}\cdot\mathbf{r}),$
(S67)

where $\mathbf{G}_{i}$ are the reciprocal lattice vectors, $\mathbf{r}$ the
position, $L_{\text{M}}$ the moiré period, $\epsilon_{r}=4$ the dielectric
constant due to hBN encapsulation, and $\rho(n)$ a filling dependent
parameter, listed in Table 1. To obtain $\rho(n)$, we fit the band structure
of the tight-binding model to the continuum model of Ref. [85] and do a self-
consistent calculation.

Values of $\rho$ for the Hartree term
---
|  $1.06^{\circ}$ |  $1.09^{\circ}$ |  $1.13^{\circ}$ |  $1.21^{\circ}$ |  $1.31^{\circ}$
$n=-2.4$ |  $-0.768$ |  $-0.776$ |  $-0.773$ |  $-0.742$ |  $-0.601$
$n=+2.4$ |  $0.768$ |  $0.795$ |  $0.793$ |  $0.756$ |  $0.612$
Table 1: Calculated values of the constant $\rho$ as a function of twist
angle, filling and dielectric constant.

The reason why we include the Hartree term only in TBG, but not in twisted
trilayer graphene (TTG), is that in TBG, we aim at a simulation with realistic
parameters, in particular with a superconducting gap of $\Delta_{\text{sc}}=1$
meV, so it is worthwhile to include the Hartree interactions, which are a
relevant correction to the non-interacting band structure.

In contrast, in the TTG nanoribbon we discuss, cutting armchair edges in the
top and bottom layers, and a ‘chiral’ edge on the middle layer, leads to a
stringent commensurability condition: the translational modulus of the chiral
layer needs to be a multiple of the translational modulus of the armchair
layers, so that the overall translational symmetry is preserved. The
translational modulus of a nanoribbon with chiral vector
$\mathbf{C}=n\mathbf{a}_{1}+m\mathbf{a}_{2}$, with
$\mathbf{a}_{1,2}=(\sqrt{3}/2,\pm 1/2)$, is given by

$T=\frac{3}{4}(t_{1}+t_{2})^{2}+\frac{1}{4}(t_{1}-t_{2})^{2}$ (S68)

with $t_{1}=(2m+n)/N_{C}$, $t_{2}=(2n+m)/N_{C}$, where
$N_{C}=\text{gcd}\left(2m+n,2n+m\right)$ and gcd stands for the greatest
common divisor. Therefore, the commensurability of the translational modulus
of a ribbon with a generic chiral vector $(n,m)$ and an armchair ribbon with
chiral vector $(n,0)$, is only fulfilled for some pairs $(n,m)@(n,0)$. One
such pair is $(-37,0)@(33,7)$, for which $T_{(33,7)}=37T_{(-37,0)}$. As a
consequence, the unit cell of the resulting TTG ribbon is large, see Fig. S15.
This precludes a simulation of the full sub-gap spectrum with a realistic
superconducting gap of $\Delta_{\text{sc}}=1$ meV, which would also require
the ribbon to be wide enough not to mix the zero modes of the left and right
edges, resulting in a system with many moiré unit cells.

Therefore, the TTG result is meant as a proof of concept: it shows that the
decomposition of TTG into TBG plus monolayer graphene works for our model.
This is evident from the fact that TTG features the four zero modes at
$k_{y}=0$ found in monolayer graphene, and that these modes indeed come the
effective monolayer that results from the odd combination of top and bottom
layers [78], as manifest in their charge distributions, which are identical in
the top and bottom layers, and zero in the middle layer, as shown in Fig. 4(e)
in the main text. Moreover, as we have emphasized, the zero modes confined to
vortex cores in twisted trilayer graphene are much more experimentally
relevant, because edge engineering is very challenging, whereas vortices occur
naturally in the superconductor and can be probed with STM measurements.

Figure S15: Unit cell of a TTG nanoribbon with chiral vector
$(-37,0)@(33,7)@(-37,0)$. The translational vector of the middle layer is
commensurate (exactly 37 times larger) with the translational vector of the
top and bottom armchair layers. Joining 5 such cells horizontally gives the
unit cell for the TTG nanoribbon for which we obtain the spectra in Fig.
4(d)-(e) in the main text.