$\psi_{+,-,A,\uparrow}$, $\psi_{+,+,B,\uparrow}$, $\psi_{+,-,B,\uparrow}$,
$\psi_{+,+,A,\downarrow}$, $\psi_{+,-,A,\downarrow}$,
$\psi_{+,+,B,\downarrow}$, $\psi_{+,-,B,\downarrow}$. We take the ansatz of
the zero modes as

$\phi(x,y)=\frac{1}{\sqrt{\mathcal{N}}}[g_{1}(r),g_{2}(r),g_{3}(r),g_{4}(r)]^{T}\otimes
v$ (101)

where $v$ is a two-component vector, $g_{1,2,3,4}(r)$ are functions of $r$ to
be determined and $\mathcal{N}$ is a normalization factor. For either
$v=(1,0)^{T}$ or $v=(0,1)^{T}$, we find that the zero mode satisfies
$g_{1}=g_{4}=0$ and

$-i\partial_{r}g_{2}+\Delta(r)g_{3}=0,\qquad\Delta(r)g_{2}+i\partial_{r}g_{3}=0.$
(102)

Combining the two equations, we obtain the equations

$\Delta(r)g_{2}-\partial_{r}\big{(}\frac{1}{\Delta(r)}\partial_{r}g_{2}\big{)}=0,\quad
g_{3}=\frac{i}{\Delta(r)}\partial_{r}g_{2}(r)$ (103)

and the solution

$g_{2}(r)=e^{-\int_{0}^{r}dr^{\prime}\Delta(r^{\prime})},\quad
g_{3}(r)=-ig_{2}(r).$ (104)

A vortex of opposite vorticity $\Delta(\mathbf{r})=e^{i\theta}\Delta(r)$ will
instead have $g_{2},g_{3}=0$ and $g_{1},g_{4}$ nonzero.

The combined symmetries, $C_{2z}^{\prime}TS=\xi_{z}\sigma_{x}K$ and
$S\mathcal{P}=-\xi_{x}\mu_{y}\sigma_{x}K$ are sufficient to protect the two
(fermionic) Abrikosov modes per valley and pin them to $0$ energy. To see
this, $S\mathcal{P}$ acts as a time-reversal symmetry and squares to $-1$,
guaranteeing Kramers’ degeneracy, while $C_{2z}^{\prime}TS$ acts as particle-
hole and squares to $+1$, forcing zero modes remain at $0$ energy, so long as
the vortex is located at a $C_{2z}$-symmetric point.

##  Appendix C Numerical verification of the boundary states of TBG-TSC

To study the edge states or corner states located around the domain walls in
the TBG, we bring Eq. (38) and Eq. (52) into the real-space form under the
basis
$(c_{\mathbf{r},l,+,\alpha,\uparrow},c^{\dagger}_{\mathbf{r},l^{\prime},-,\beta,\downarrow})^{T}$,

$\displaystyle\mathcal{H}^{(+,0)}(\mathbf{r})$ $\displaystyle=$ $\displaystyle
v_{F}k_{D}\xi_{z}\left(\begin{array}[]{cc}-i\partial_{\mathbf{r}}\cdot\bm{\sigma}&0\\\
0&-i\partial_{\mathbf{r}}\cdot\bm{\sigma}\\\ \end{array}\right)$ (110)
$\displaystyle+\xi_{z}\left(\begin{array}[]{cc}0&T(\mathbf{r})\\\
T^{\dagger}(\mathbf{r})&0\\\ \end{array}\right)$
$\displaystyle\mathcal{H}^{(+,1)}(\mathbf{r})$ $\displaystyle=$
$\displaystyle\Delta(\mathbf{r})\xi_{x},$ (111)

where $T(\mathbf{r})=\sum_{j=1}^{3}e^{-i\mathbf{q}_{j}\cdot\mathbf{r}}T^{j}$,
and $\sigma$, $\xi$, and the explicitly written $2\times 2$ matrices denote
sublattice, particle-hole, and the layer degrees of freedom, respectively. The
continuous position $\mathbf{r}$ is the coarse-grained graphene lattice. In
the following subsections, we will consider different geometries of dividing
the material into pairing domains with positive and negative
$\Delta(\mathbf{r})$, and discuss whether the edge states or corner states
would appear.

###  C.1 Edge states

To solve edge states in the direction perpendicular to $x$ ($y$), we take the
periodic boundary condition in the $x$ ($y$) direction with $x$ ($y$) ranging
from $-L/2$ to $L/2$. Then we assign positive and negative $\Delta(x)$ for
$|x|<L/4$ ($|y|<L/4$) and $|x|>L/4$ ($|y|>L/4$), respectively, as shown in
Fig. 6(b)(c). Two domain walls that bind zero modes are formed at $x=\pm L/4$
($y=\pm L/4$). Using periodic boundary condition greatly simplifies the
calculation as the basis states are plane waves. The matrix elements of the
real-space Hamiltonian $H(\mathbf{r})$ are calculated with respect to plane
waves $\ket{\mathbf{k},\mathbf{G},i}$, where $\mathbf{k}$ indexes the Moiré
Brillouin zone momentum, $\mathbf{G}$ the element of the $\mathbf{Q}$ lattice,
and $i$ the remaining indices of the TBG Hamiltonian, which are sublattice,
valley, Moiré valley, and particle-hole. The Hamiltonian matrix

$\displaystyle\braket{\mathbf{k},\mathbf{G},i}{H(\mathbf{r})}{\mathbf{k},\mathbf{G}^{\prime},i^{\prime}}$
(112)

is then numerically diagonalized.

Figure 6:  (a) Definitions of $\mathbf{a}_{\rm Mx}$ and $\mathbf{a}_{\rm My}$,
the Moiré lattice vectors. (b)(c) Pairing domain geometries used to solve for
edge states in the $x$ and $y$ directions, respectively. The dashed squares
denotes supercells in these geometries, and the specific sizes used in the
numerics are indicated. The orange (blue) regions denote domains with
$\Delta(\mathbf{r})=\Delta_{0}>0$ ($\Delta(\mathbf{r})=-s\Delta_{0}<0$).
Domain walls that bind zero modes are formed at the boundaries between
positive and negative pairings. $s>0$ controls the relative magnitude of the
pairing potential between the positive and negative regions, and is designed
to break accidental symmetries Eq. (113), (114). $O_{1}$ denotes an
“inversion” center of the $P$ symmetry, and $O_{2}$ denotes an “inversion”
center of the $P_{2}$ symmetry. (d)-(o) Edge bands in the $x$ direction.
(p)-(aa) Edge bands in the $y$ direction. Three different twisted angles are
investigated: $\theta=1.05^{\circ},1.20^{\circ}$, and $1.40^{\circ}$. Four
different sets of parameters are chosen: $s=1$, $\beta=0$ (1st row); $s=0.4$,
$\beta=0$ (2nd row); $s=1$, $\beta=20$meV (3rd row); $s=1$, $\beta=40$meV (4th
row). In the calculation we have chosen the parameters $w_{1}=110$meV,
$w_{0}=0.8w_{1}$, $v_{F}=5.944{\rm eV\cdot\mathring{A}}$, and
$|\Delta_{0}|=20$meV.

Using plane waves as the basis of Eq. (110), we obtain the edge bands in the
$x$ direction (i.e. when the pairing domain walls are set perpendicular to the
$x$ direction) at three twist angles $\theta=1.05^{\circ}$, $1.2^{\circ}$,
$1.4^{\circ}$, as shown in Fig. 6(d)-(f), as well as the edge bands in the $y$
direction, as shown in Fig. 6(p)-(r). It is seen that, gapless edge states not
only appear in the $x$ direction, as analyzed in the main text, but also
appear accidentally in the $y$ direction, where no known symmetry could
protect such gapless modes.

The origin of the accidental edge states in the $y$ direction may be ascribed
to an additional emergent symmetry on the domain wall, which we would denote
by $P_{2}\equiv\xi_{z}P$, where
$P=\delta_{\mathbf{Q,-Q^{\prime}}}\zeta_{\mathbf{Q}}\sigma_{0}\xi_{z}$ is the
unitary symmetry defined in Appendix A that anti-commutes with the bulk
Hamiltonian,

$\displaystyle
P\mathcal{H}^{(+,0)}(\mathbf{r})P^{-1}=-\mathcal{H}^{(+,0)}(-\mathbf{r})$
(113) $\displaystyle
P\mathcal{H}^{(+,1)}(\mathbf{r})P^{-1}=-\mathcal{H}^{(+,1)}(-\mathbf{r})$
(114)

where $\mathcal{H}^{(\eta,1)}(\mathbf{r})=\xi_{z}{H}^{(+)}(\mathbf{r})$ is the
free part, and $\mathcal{H}^{(\eta,1)}(\mathbf{r})=\Delta_{0}\xi_{x}$ is a
uniform spin-singlet pairing potential. However, at the domain wall,
$\Delta(\mathbf{r})$ no longer stays constant, but has opposite signs at the
two sides, i.e., $\Delta(\mathbf{r})=-\Delta(\mathbf{-r})$. (Here
$\mathbf{r}=0$ corresponds to a point at the domain wall.) Therefore, $P$ is
broken by $\mathcal{H}^{(+,1)}(\mathbf{r})$ by the domain wall. However, both
$\mathcal{H}^{(+,0)}(\mathbf{r})$ and $\mathcal{H}^{(+,1)}(\mathbf{r})$ anti-
commute with $P_{2}$ instead. In real space, $P_{2}$ is $P$, which inverts
position $\mathbf{r}\rightarrow-\mathbf{r}$ and switches layer, combined with
$\xi_{z}$, which introduces a relative sign of the transformation depending on
whether the operator being transformed is a particle or a hole.

Within the Dirac theory of the edge states, we can show that this emergent
symmetry $P_{2}$ can indeed protect edge states in any arbitrary direction. As
has been presented in Appendix B.1, one can show that
$P_{2}=i\xi_{0}^{\prime}\mu_{y}$ on the edge by definition, and that
$P_{2}\mathcal{H}_{\text{eff}}(k_{y})P_{2}=-\mathcal{H}_{{\rm eff}}(-k_{y})$.
Therefore, the only gap term allowed by both $S$ and $\mathcal{P}$,
$\xi_{x}^{\prime}\mu_{0}$, would not obey $P_{2}$, in the sense that
$\\{P_{2},\xi_{x}^{\prime}\mu_{0}\\}\not=0$, thus not permitted. Therefore,
although $C_{2x}$ is absent in an arbitrary direction, the emergent $P_{2}$
would still rule out other possible gap terms.

Enlightened by this, we tried several methods to break the $P_{2}$ symmetry in
order to gap the edge states in the $y$ direction. One natural way is to tune
the relative magnitude between the negative and the positive domains, so that

$\displaystyle\Delta(\mathbf{r})$ $\displaystyle=$
$\displaystyle\Delta_{0},~{}~{}~{}~{}\text{for the positive domains}$ (115)
$\displaystyle\Delta(\mathbf{r})$ $\displaystyle=$
$\displaystyle-{s}\Delta_{0},~{}~{}~{}\text{for the negative domains}$ (116)

The numerical results are shown in Fig. 6(g)-(i)(s)-(u). Only at
$\theta=1.2^{\circ}$, the edge state in the $y$ direction has a non-negligible
gap. However, the edge states in the $y$ direction at
$\theta=1.05^{\circ},1.4^{\circ}$ are still almost gapless, meaning our
perturbation does not do a very good job of breaking the accidental symmetry.
Another term that breaks $P_{2}$ is given by

$\mathcal{H}^{(+,2)}=f(\mathbf{r})\xi_{z}\left(\begin{array}[]{cc}0&\sigma_{0}\\\
\sigma_{0}&0\\\ \end{array}\right)$ (117)

where $f(\mathbf{r})$ is odd on the domain wall,

$\displaystyle f(\mathbf{r})=+\beta,~{}\text{if}~{}\Delta(\mathbf{r})>0$ (118)
$\displaystyle f(\mathbf{r})=-\beta,~{}\text{if}~{}\Delta(\mathbf{r})<0$ (119)

With different magitude of $\beta$, we plot the edge bands in Fig.
6(j)-(o)(v)-(aa). A significant gap in the $y$ direction is opened and the gap
increases with $\beta$, while the zero modes in the $x$ direction, as
expected, are stable against this perturbation. Note that the low-lying states
are always gapped from the bulk; this is expected as the zero modes are
pointlike: instead of giving a continuous set of energies like an edge, the
energies from the zero modes are discrete.

###  C.2 Corner states

Now that the accidental edge states are gapped out, we can numerically
investigate the evolution of the symmetry-protected corner modes when the
chemical potential is tuned to finite values. The largest gap in the
$y$-direction appears at $\theta=1.40^{\circ}$ with $\beta=40$meV (Fig. 6). We
hence choose $\theta=1.4^{\circ}$, $\beta=40$meV in the following numerical
calculations to demonstrate the existence of the corner states. Two different
pairing domain geometries are devised. One is the rhombic domains as shown in
Fig. 7(a), chosen to break the $C_{3z}$ symmetry, and the other is the
circular domains as shown in Fig. 8(a).

In the first geometry, when $\mu=0$, as illustrated in Fig. 7(b), there are 4
zero modes that are well separated from other states in the spectrum, and
their total spatial densities are localized at the two $C_{2x}$-invariant
corners, which shows perfect accordance with the theoretical prediction in the
main text. As $\mu$ is tuned up, as illustrated in Fig. 7(c), the zero modes
slightly “diffuse” along the edges, but are still well separated from other
states in energy. At some critical value (Fig. 7(d)), the bulk states collapse
onto the zero modes in the spectrum, and for chemical potential beyond the
critical value, zero modes disappear, and the spatial densities of the 4
lowest states start to spread into the bulk (Fig. 7(e)). The reason the zero
modes spread, and not move, is due to the straight edges of the rhombic
geometry.

Figure 7: Rhombic pairing domains. The periodic boundary condition is realized
by a supercell spanned by $18\mathbf{a}_{M1}$ and $18\mathbf{a}_{M2}$. (a) The
configuration of the domains. The dashed square indicates the spatial region
displayed in the following spatial density diagrams. (b)-(e) energy spectra
near the zero-energy and the spatial densities of the 4 lowest states. The red
dashed rhomboids indicate the the domain walls. Results for
$\mu=0,10,16,22$meV are displayed. Through all the numerics, we choose
$w_{1}=110$meV, $w_{0}=0.8w_{1}$, $v_{F}=5.944{\rm eV\cdot\mathring{A}}$,
$\theta=1.40^{\circ}$, $\beta=40$meV, and $|\Delta_{0}|=20$meV.

To see how zero modes from different $C_{3z}$ corners evolve with chemical
potential and annihilate each other, a second geometry is investigated, as
shown in Fig. 8(a). When $\mu=0$ (Fig. 8(b)), there are 12 zero modes per
valley well separated from other states in the spectrum, which corresponds to
3 copies of the states in the rhombic geometry, and their spatial
distributions are localized at the two $C_{2x}$-invariant corners and four
other equivalent corners related by $C_{3z}$. As $\mu$ is tuned up (Fig.
8(c)), while their energies are still pinned close to zero and separated from
other states, they start moving towards their adjacent corners. At some
critical value, states from adjacent corners meet (Fig. 8(d)), and
correspondingly, the gap between the zero modes and other states vanish
(<1meV). For chemical potential beyond this critical value (Fig. 8(e)), zero
modes from adjacent corners gap each other out.

The numerical results are highly consistent with the theoretical predictions
presented in the main text.

Figure 8: Circular pairing domains. The periodic boundary condition is
realized by a supercell spanned by $36\mathbf{a}_{M1}$ and
$36\mathbf{a}_{M2}$. (a) The configuration of the domains. The dashed square
indicates the spatial region displayed in the following spatial density
diagrams. (b)-(e) energy spectra near the zero-energy and the spatial
densities of the 12 lowest states. The red dashed circles indicate the the
domain walls. Results for $\mu=0,8,12,20$meV are displayed. Through all the
numerics, we choose $w_{1}=110$meV, $w_{0}=0.8w_{1}$, $v_{F}=5.944{\rm
eV\cdot\mathring{A}}$, $\theta=1.40^{\circ}$, $\beta=40$meV, and
$|\Delta_{0}|=20$meV.

##  Appendix D Other pairing symmetries

In Appendices A, B and C we have analytically and numerically shown that the
TBG Hamiltonian with a uniform spin-singlet pairing is a topological
superconductor guaranteed by the $C_{2z}T$ and $\mathcal{P}$ symmetries. Now
we show that the superconducting phase with the $B_{1}$ triplet pairing
(defined in the next paragraph) is also topological and has the similar edge
and corner states as the uniform spin-singlet pairing.

Only four intervalley terms open full gaps in the Bogoliubov bands: $A_{1}$
(spin-singlet), $A_{2}$ (spin-singlet), $B_{1}$ (spin-triplet), $B_{2}$ (spin-
triplet) [164]. $A_{1},A_{2}$ and $B_{1},B_{2}$ are even and odd under the
$C_{2z}$ operation, respectively. And, $A_{1},B_{1}$ and $A_{2},B_{2}$ are
even and odd under the $C_{2x}$ operation, respectively. The four terms are
invariant under the other crystalline symmetry operations. We find that the
$A_{1}$ and $B_{1}$ pairings are topological whereas the $A_{2}$ and $B_{2}$
pairings are trivial. In fact, the $B_{1}$ pairing is very similar to the
uniform spin-singlet ($A_{1}$) pairing: their BdG Hamiltonians in a single
valley are the same, but the sign of the pairing switches for the opposite
valley. Explicitly, the pairing forms read
$\xi_{x}\tau_{0}\mu_{0}\sigma_{0},\xi_{x}\tau_{z}\mu_{z}\sigma_{0},\xi_{x}\tau_{z}\mu_{0}\sigma_{0},\xi_{x}\tau_{0}\mu_{z}\sigma_{0}$,
for the $A_{1},A_{2},B_{1},B_{2}$ pairings, respectively.

###  D.1 The $A_{2}$ singlet pairing

In this subsection we show the triviality of the $A_{2}$ pairing. We assume
the valley-U(1) symmetry and use the same BdG basis (Eq. 73) as in the Dirac
theory of the $A_{1}$ pairing. As discussed in the main text and in Section
A.3, due to the time-reversal and spin-SU(2) symmetries, the Dirac theory of
any spin-singlet pairing takes the form

$\mathcal{H}=k_{x}\xi_{z}\tau_{z}\mu_{0}\sigma_{x}+k_{y}\xi_{z}\tau_{0}\mu_{0}\sigma_{y}+\xi_{x}\otimes\Delta,$
(120)

with $\Delta$ being an eight by eight Hermitian matrix, built of
$\tau,\mu,\sigma$ operators. Here $\tau$, $\mu$, $\sigma$ are Pauli (identity)
matrices for the valley, Moiré valley, and sublattice degrees of freedoms,
respectively. For $A_{1}$, the pairing matrix
$\Delta=\tau_{0}\mu_{0}\sigma_{0}=I$. Anti-commutation with $\xi_{y}$ forces
the pairing term to $\xi_{x}\otimes\Delta.$

Due to the valley-U(1) symmetry, $\Delta$ is diagonal in the valley index, as
Eq. 33 combines creation operators of one valley with annihilation operators
of the other. The symmetry operators in the BdG formalism are given in Eq. 74.
Among these symmetries, $S$ and $\mathcal{P}$ are constraints respectively
given by the spinful time-reversal symmetry and the charge conjugation
$\mathcal{P}_{c}$ (Section A.3), and $C_{3z}$, $C_{2x}$, $T$, $C_{2z}T$ are
crystalline symmetries. For simplicity, we only consider the pairings that are
diagonal in the Moiré valley index $\mu_{0}$, i.e., pairings preserving the
Moiré translation symmetry. The $A_{2}$ pairing is invariant under $C_{3z}$,
$C_{2z}$, $C_{2z}T$ and is odd under $C_{2x}$. For the pairing to be invariant
(with eigenvalue $1$) under $C_{3z}$
($e^{i\frac{2\pi}{3}\tau_{z}{\mu_{0}}\sigma_{z}}$), $\Delta$ must be
proportional to $\sigma_{0}$ or $\sigma_{z}$; For it to be odd under $C_{2x}$
(${\tau_{0}}\mu_{x}\sigma_{x}$), $\Delta$ can be chosen as
$\tau_{0,z}\mu_{z}\sigma_{0}$ or $\tau_{0,z}\mu_{0}\sigma_{z}$; For it to be
invariant under $C_{2z}$ ($\tau_{x}\mu_{x}\sigma_{x}$), only
$\tau_{z}\mu_{z}\sigma_{0}$ remains. The term $\tau_{z}\mu_{z}\sigma_{0}$ also
preserves the $C_{2z}T$ (${\tau_{0}\mu_{0}}\sigma_{x}K$) symmetry. This form
of the $A_{2}$ pairing is also given in Table II of Ref. [164]. Therefore, the
BdG Hamiltonian on the basis Eq. 73 reads

$\mathcal{H}=k_{x}\xi_{z}\tau_{z}\mu_{0}\sigma_{x}+k_{y}\xi_{z}\tau_{0}\mu_{0}\sigma_{y}+\Delta_{0}\xi_{x}\tau_{z}\mu_{z}\sigma_{0}.$
(121)

The $C_{2x}$ and $\mathcal{P}$ (given by the charge conjugation
$\mathcal{P}_{c}$ Section A.3) are broken by the $A_{2}$ pairing.
Nevertheless, there exists a projective $C_{2x}$ and a projective
$\mathcal{P}$

$C_{2x}^{\prime}=\xi_{z}\tau_{0}\mu_{x}\sigma_{x},\qquad\mathcal{P}^{\prime}=i\xi_{0}\tau_{0}\mu_{y}\sigma_{x}K\
,$ (122)

which correspond to $C_{2x}$ followed by a charge-U(1) rotation and
$\mathcal{P}_{c}$ followed by a charge-U(1) rotation, respectively [165]. We
will show these projective symmetries do not protect topology.

Under a charge rotation $c\rightarrow e^{i\gamma}c,c^{\dagger}\rightarrow
e^{-i\gamma}c^{\dagger}$, $\psi\rightarrow e^{i\gamma\xi_{z}}\psi$, as $\psi$
consists of both creation and annihilation operators. The charge rotation
implemented is $\gamma=\pi/2$, with additional factors of $i$ placed to
preserve the algebra of the operators $C_{2x}^{\prime}$ and
$\mathcal{P}^{\prime}$. Following the calculation in Section B.1, we solve for
the effective edge theory (in valley $\eta=+$) in the $x$ direction. First
take ansatz $\phi_{k_{y}}(x)\propto e^{ik_{y}y-\lambda|x|}u(k_{y})$ and solve
the Hamiltonian with domain wall in $\Delta=-\Delta_{0}\text{sgn}(x)$ for
$k_{y}=0$, $\Delta_{0}>0$:

$\displaystyle{\cal
H}=-i\partial_{x}\xi_{z}\tau_{z}\mu_{0}\sigma_{x}+\Delta_{0}\text{sgn}(x)\xi_{x}\tau_{z}\mu_{z}\sigma_{0}$
(123) $\displaystyle
i\lambda\text{sgn}(x)\xi_{z}\tau_{z}\mu_{0}\sigma_{x}u-\Delta_{0}\text{sgn}(x)\xi_{x}\tau_{z}\mu_{z}\sigma_{0}u=0$
(124) $\displaystyle\lambda u=\Delta_{0}\xi_{y}\mu_{z}\sigma_{x}u.$ (125)

These equations have solutions when $\lambda=\Delta_{0}$,
$\xi_{y}\mu_{z}\sigma_{x}=1$. Notice that the valley operator $\tau$ drops out
of the solution, allowing us to focus on a single sector $\tau_{z}=+1$:

$\displaystyle u_{1}$ $\displaystyle=\frac{1}{2}\begin{bmatrix}1\\\
i\end{bmatrix}\otimes\begin{bmatrix}1\\\
0\end{bmatrix}\otimes\begin{bmatrix}1\\\
1\end{bmatrix},\;u_{2}=\frac{1}{2}\begin{bmatrix}1\\\
i\end{bmatrix}\otimes\begin{bmatrix}0\\\
1\end{bmatrix}\otimes\begin{bmatrix}1\\\ -1\end{bmatrix},$ $\displaystyle
u_{3}$ $\displaystyle=\frac{1}{2}\begin{bmatrix}-1\\\
i\end{bmatrix}\otimes\begin{bmatrix}1\\\
0\end{bmatrix}\otimes\begin{bmatrix}1\\\
-1\end{bmatrix},\;u_{4}=\frac{1}{2}\begin{bmatrix}-1\\\
i\end{bmatrix}\otimes\begin{bmatrix}0\\\
1\end{bmatrix}\otimes\begin{bmatrix}1\\\ 1\end{bmatrix}.$ (126)

We find

$\mathcal{H}=k_{y}\xi_{y}^{\prime}\mu_{z}\ ,$ (127)

with the projected symmetries on edge

$\displaystyle
S=\xi_{z}^{\prime},\quad\mathcal{P}^{\prime}=-i\xi_{y}^{\prime}\mu_{x}K,\quad
C_{2x}^{\prime}=\xi_{y}^{\prime}\mu_{y}.$ (128)

The gap term $\xi_{y}^{\prime}\mu_{y}$ is allowed by the symmetries and hence
there is no protected edge state. On the edge in the $y$ direction, we can
obtain the same effective Hamiltonian with the projected symmetries

$S=\xi_{z}^{\prime},\quad\mathcal{P}^{\prime}=i\xi_{y}^{\prime}\mu_{x}K,\quad
C_{2x}^{\prime}C_{2z}T=\xi_{0}^{\prime}\mu_{x}K.$ (129)

The action of the symmetries on the edge in the $y$ direction is found by
rotating the wavefunctions
$u_{i}^{\prime}=e^{-i\frac{\pi}{4}\sigma_{z}}u_{i}$. The prefactor
$e^{-i\frac{\pi}{4}\sigma_{z}}$ commutes with both $S,\mathcal{P}^{\prime}$,
so no change occurs to those symmetries, but $C_{2y}T$ is modified. The gap
term $\xi_{x}^{\prime}\mu_{0}$ is allowed by the symmetries and hence there is
no protected edge state. For edges in generic directions the only symmetries
on the edge are $S$ and $\mathcal{P}^{\prime}$, which, as shown in the $x$ and
$y$ directions, cannot forbid gap terms. Therefore, there is no gapless edge
in any direction.

There is also no corner state. Unlike the $A_{1}$ pairing, where
$[S,\mathcal{P}]=0$, for the $A_{2}$ pairing we have
$\\{S,\mathcal{P}^{\prime}\\}=0$. Thus, even there existed zero energy doublet
at a corner, the two states in the doublet must have opposite chiral
eigenvalues and hence could be gapped out. To be specific, we can choose the
$S$ and $\mathcal{P}^{\prime}$ operators in the doublet space as $\sigma_{z}$
and $i\sigma_{y}K$, respectively. Then both $\sigma_{x}$ and $\sigma_{y}$ gap
terms are allowed to appear.

From the absence of gapless boundary mode, we conclude the triviality of the
$A_{2}$ pairing phase.

###  D.2 The $B_{1}$ and $B_{2}$ spin triplet pairings

We assume the valley-U(1) symmetry and the spinful time-reversal symmetry in
the triplet pairing. The triplet pairing must break the spin-SU(2) symmetry,
as the $S_{z}=\pm 1,0$ sectors will rotate into each other. For example,
rotating spins by $S_{z}$ will introduce phases into the $S_{z}=\pm 1$
sectors. Here we assume a spin-U(1) remains. Without loss of generality, we
assume the rotation axis of the spin-U(1) in the $z$ direction. Due to the
valley-U(1) and spin-U(1), the operator
$c^{\dagger}_{\mathbf{k},\eta,l,\alpha,\uparrow}$ can only couple to
$c^{\dagger}_{-\mathbf{k},-\eta,l,\alpha,\downarrow}$ in the pairing (we are
focusing on the $S_{z}=0$ sector of the triplet pairing). Hence we can still
use the BdG basis Eq. 73 for the triplet pairings. Due to the spinful time-
reversal symmetry, the BdG Hamiltonian still has the chiral symmetry
$\\{\xi_{y},\mathcal{H}^{(\eta)}\\}=0$ (Eq. 43). We can write the Dirac BdG
Hamiltonian as

$\mathcal{H}=k_{x}\xi_{z}\tau_{z}\mu_{0}\sigma_{x}+k_{y}\xi_{z}\tau_{0}\mu_{0}\sigma_{y}+\xi_{x}\otimes\Delta,$
(130)

with the symmetry operators given in Eq. 74. Due to the valley-U(1), $\Delta$
must be proportional to $\tau_{0,z}$. For simplicity, we only consider
$\Delta$’s that are diagonal in the Moiré valley index, i.e., pairings that
preserve the Moiré translation symmetry. Since the triplet pairing is even
under the spinful time-reversal ($i\hat{s}_{y}T$) and odd under the
$i\hat{s}_{y}$ spin-rotation, the triplet pairing must be odd under the
spinless time-reversal symmetry $T=\tau_{x}\mu_{x}K$. Thus the options for
$\Delta$ include: $\tau_{z}\mu_{0}\sigma_{0,x,z}$,
$\tau_{z}\mu_{z}\sigma_{y}$, $\tau_{0}\mu_{z}\sigma_{0,x,z}$,
$\tau_{0}\mu_{0}\sigma_{y}$.

We now determine the forms of the $B_{1}$ and $B_{2}$ pairings. For the
pairing to be invariant under $C_{3z}=e^{i\frac{2\pi}{3}\tau_{z}\sigma_{z}}$,
$\Delta$ must be proportional to $\sigma_{0}$ or $\sigma_{z}$, hence only
$\tau_{z}\mu_{0}\sigma_{0,z}$ and $\tau_{0}\mu_{z}\sigma_{0,z}$ are valid.
Because the $B_{1}$ and $B_{2}$ by definition are odd under $C_{2z}$ and they
are also odd under the spinless time-reversal ($T$) as discussed in the last
paragraph, they must be even under $C_{2z}T$. For it to be invariant under
$C_{2z}T=\sigma_{x}K$, only $\tau_{z}\mu_{0}\sigma_{0}$ and
$\tau_{0}\mu_{z}\sigma_{0}$ are valid. $\tau_{z}\mu_{0}\sigma_{0}$ and
$\tau_{0}\mu_{z}\sigma_{0}$ are even and odd under $C_{2x}=\mu_{x}\sigma_{x}$,
thus they form the $B_{1}$ and $B_{2}$ representations, respectively. These
forms of $B_{1}$ and $B_{2}$ pairings are also given in Table VI of Ref.
[164]. The BdG Hamiltonian of the $B_{1}$ pairing and the symmetries are given
by

$\mathcal{H}=k_{x}\xi_{z}\tau_{z}\mu_{0}\sigma_{x}+k_{y}\xi_{z}\tau_{0}\mu_{0}\sigma_{y}+\Delta_{0}\xi_{x}\tau_{z}\mu_{0}\sigma_{0},$
(131)

and Eq. 74, respectively. Since the Hamiltonian and symmetries in the
$\eta=+1$ block are the same as those of the $A_{1}$ pairing (Eqs. 72 and 74),
the $B_{1}$ pairing must have the same topology and the same boundary states
as the $A_{1}$ pairing. An identical argument holds in the opposite $\eta=-1$
valley.

The BdG Hamiltonian of the $B_{2}$ pairing is

$\mathcal{H}=k_{x}\xi_{z}\tau_{z}\mu_{0}\sigma_{x}+k_{y}\xi_{z}\tau_{0}\mu_{0}\sigma_{y}+\Delta_{0}\xi_{x}\tau_{0}\mu_{z}\sigma_{0}.$
(132)

It breaks the single valley symmetries $C_{2x}=\mu_{x}\sigma_{x}$ and
$\mathcal{P}=i\xi_{z}\mu_{y}\sigma_{x}K$. Nevertheless, there exists a
projective $C_{2x}$ and a projective $\mathcal{P}$

$C_{2x}^{\prime}=\xi_{z}\tau_{0}\mu_{x}\sigma_{x},\qquad\mathcal{P}^{\prime}=i\xi_{0}\tau_{0}\mu_{y}\sigma_{x}K\
,$ (133)

which correspond to $C_{2x}$ and $\mathcal{P}_{c}$ (Section A.3) followed by a
charge-U(1) rotation $c\rightarrow ic$, respectively. The Hamiltonian and
symmetries in the $\eta=+1$ block are the same those of the $A_{2}$ pairing
(Section D.1), thus the $B_{2}$ pairing must have the same topology as the
$A_{2}$ pairing, which is trivial.

In the above we have shown that the Dirac BdG Hamiltonians of the $B_{1}$ and
$B_{2}$ pairings are the same as the Dirac BdG Hamiltonians of $A_{1}$ and
$A_{2}$ pairings respectively, up to a sign depending on the valley (which
will not affect the physics of TBG because the valleys are independent). In
fact, this correspondence between singlet pairing and triplet pairing is in
general true if the TBG Hamiltonian has two independent spin-SU(2) symmetries
in the two valleys, which is usually assumed in theoretical works on TBG.
Suppose we have a singlet ($\uparrow\downarrow-\downarrow\uparrow$) pairing

$\sum_{\mathbf{k}}\sum_{\eta}\sum_{ll^{\prime}\alpha\beta}\Delta_{l,\alpha;l^{\prime},\beta}^{(\eta)}(\mathbf{k})c^{\dagger}_{\mathbf{k},\eta,l,\alpha,\uparrow}c^{\dagger}_{-\mathbf{k},-\eta,-l^{\prime},\beta,\downarrow},$
(134)

where $\Delta$ is Hermitian due to spinful time-reversal symmetry
($i\tau_{x}\mu_{x}s_{y}K$) and satisfies
$\Delta^{(\eta)}_{l,\alpha;l^{\prime},\beta}(\mathbf{k})=\Delta^{(-\eta)}_{-l^{\prime},\beta;-l,\alpha}(-\mathbf{k})$
due to the $i\hat{s}_{y}$ spin-rotation. The pairing must be degenerate under
any charge and spin rotations in the two valleys separately. We apply charge-
spin rotation $e^{i\frac{\pi}{4}(s_{0}-s_{z})}$ in the valley $\eta=+$ and
$e^{i\frac{\pi}{2}(s_{0}+s_{z})}$ in the valley $\eta=-$ and obtain

$\sum_{\mathbf{k}}\sum_{\eta}\sum_{ll^{\prime}\alpha\beta}\eta\Delta_{l,\alpha;l^{\prime},\beta}^{(\eta)}(\mathbf{k})c^{\dagger}_{\mathbf{k},\eta,l,\alpha,\uparrow}c^{\dagger}_{-\mathbf{k},-\eta,-l^{\prime},\beta,\downarrow}.$
(135)

The rotated pairing becomes triplet $(\uparrow\downarrow+\downarrow\uparrow)$
since it is odd under $i\hat{s}_{y}$. Therefore, the triplet pairings and
singlet pairings are related by a unitary transformation. The BdG Hamiltonian
in the valley $\eta=+$ is invariant under this transformation. Thus the
single-valley symmetries and topologies of the singlet pairing and the
corresponding rotated triplet pairing must be the same.

##  Appendix E Valley-U(1) breaking on a hard edge

In this section we discuss the valley-U(1) symmetry, and how in a system of
TBG with no pairing either the valley symmetry or the anti-unitary particle-
hole $\mathcal{P}$ _must_ be strongly broken on the edge.

Consider a circular sample of TBG from radius $0<r<a$. The Hamiltonian per
spin sector without pairing reads

$\displaystyle
H=k_{x}\tau_{z}\mu_{0}\sigma_{x}+k_{y}\tau_{0}\mu_{0}\sigma_{y},$ (136)

with the symmetries

$\displaystyle C_{2z}T=\sigma_{x}K,~{}~{}{\cal P}=i\mu_{y}\sigma_{x}K.$ (137)

Outside the sample is a vacuum and is fully gapped. The vacuum must respect
$C_{2z}T$ and we assume that the vacuum also respects $\mathcal{P}$. Then
there must be terms consistent with the symmetries in Eq. 137 to gap out the
Dirac terms _outside_ the sample to properly model the symmetry-preserving
vacuum.

Working in graphene valley $\tau_{z}=+1$, we wish to gap the spectrum by
adding terms that anti-commute with the Hamiltonian: $\mu_{i}\sigma_{z}$.
Terms that commute with $H$ and affect Moiré valley only, like $\mu_{x,z}$, do
not open up gaps; they simply shift the Dirac crossing points away from each
other. Forcing consistency with $\cal{P}$ further restricts the gap terms to
the form

$\displaystyle M_{1}(\mathbf{r})\mu_{0}\sigma_{z},$ (138)

which to be symmetric under $C_{2z}T$ forces
$M_{1}(\mathbf{r})=-M_{1}(-\mathbf{r})$, with origin at $\mathbf{r}=0$. Adding
this term will gap out the Dirac cones, except along lines where
$M_{1}(\mathbf{r})$ vanishes. By $M_{1}(\mathbf{r})=-M_{1}(-\mathbf{r})$ and
the intermediate value theorem, there must be at least two lines extending
from the sample to infinity that are gapless, leading to a contradiction. For
a simple way to see the gapless lines, take $\mathbf{r}=re^{i\theta}$, and
consider a circle of constant radius $r_{0}$ so that $\mathbf{r}$ lies outside
the sample of TBG. $C_{2z}T$ symmetry forces
$M_{1}(r_{0},\theta)=-M_{1}(r_{0},\theta+\pi)$. If
$M_{1}(r_{0},\theta_{0})>0$, then $M_{1}(r_{0},\theta_{0}+\pi)<0$, and there
must be at least two diametric points along the circle of constant radius
$r_{0}$ where $M_{1}$ vanishes. This argument holds for all $r_{0}$ outside
the sample of TBG.

What this analysis implies is that outside the TBG sample, in order to
reproduce the fully gapped vacuum that respects $\mathcal{P}$ we must involve
terms that mix the valleys: e.g.
$M_{2}(\mathbf{r})\tau_{x}\mu_{x,z}\sigma_{x}$. Alternatively, breaking
$\mathcal{P}$ allows the term $M_{3}(\mathbf{r})\mu_{y}\sigma_{z}$. Both
$M_{2}(\mathbf{r})$ and $M_{3}(\mathbf{r})$ will gap out the Dirac cone
everywhere, because they are symmetric under $C_{2z}T$ and is not required to
vanish anywhere in the vacuum. This means that along the edge of TBG, at those
points where $M_{1}(\mathbf{r})$ vanishes we must gap out the system instead
with $M_{2}$ or $M_{3}$, and will strongly break either valley-$U(1)$ or
$\mathcal{P}$.

##  Appendix F Shiba States

We study the response of TBG-TSC to impurities. Similar techniques involving
Shiba states [166] have been proposed to detect topological superconductivity
in 1D and topological insulators in higher dimensions [167, 102], and even
build topological phases [168]. These works study the Green’s function of the
Hamiltonian and how it reacts to delta-function impurities. However, the Dirac
version of our Hamiltonian has a problem: it does not appropriately
distinguish between the trivial and topological phases on its own. The Dirac
equation possess a symmetry $\Delta\rightarrow-\Delta$, with $\Delta$ the
mass, while the two phases are meant to be distinct. The standard technique is
to modify the Dirac Hamiltonian mass (i.e. superconducting gap) to $\Delta-
mk^{2}$ for some mass $m$ (Section A.4) [169]. This term is employed to
properly calculate, for example, the Chern number for the Dirac theory of
topological insulators and superconductors: with $m<0$ for $\Delta>0$, the
system is trivial, and with $m>0$, the system is topological [169]. Notice
that the topological phase has a band inversion while the trivial phase does
not. In this Appendix, we calculate the behavior of Shiba states for the Dirac
Hamiltonian with quadratic correction, and demonstrate a qualitative
difference between the $m<0$ and $m>0$ phases.

For TBG-TSC, however, we argue that $m<0$ and $m>0$ _do not_ correspond to the
topological and trivial phases. The distinction between topological and
trivial phases works for non-anomalous systems, like a Chern insulator in two
dimensions. In that situation $m<0$ and $m>0$ indeed do correspond to
topological and trivial phases. In TBG-TSC however, the single-valley
Hamiltonian is anomalous, analogous to the surface of a 3DTI. The two distinct
phases correspond to $\Delta>0$ and $\Delta<0$ (not $m$), and even then are
not topological/trivial themselves but only have topological distinction when
compared to one another. Properly capturing the topological-trivial
distinction requires $|m\Lambda^{2}|\gg\Delta$, with $\Lambda$ being the UV
cutoff for $k$. In this limit the band inversion is present and the
Hamiltonian is well-defined even as $k\rightarrow\infty$.

The opposite limit, $|m\Lambda^{2}|\ll|\Delta|$, does not have a well-defined
large momentum limit. Only the difference between $\Delta<0$ and $\Delta>0$
has a well-defined topological invariant. In this work we have assumed the
second limit because the nearly uniform spin-singlet pairing (Section A.4) has
no $\mathbf{k}$-dependence or very weak $\mathbf{k}$-dependence (Section A.4):
$\Delta-m\mathbf{k}^{2}$ is close to being uniform, so
$m\mathbf{k}^{2}\ll|\Delta|$ for all $\mathbf{k}$ in the BZ. This means that
while Shiba states do distinguish between $m<0$ and $m>0$ phases (the former
corresponding to singlet-pairing TBG-TSC), it does not offer a distinguishing
signature of TBG-TSC. Nevertheless, we examine the behavior of impurities in
the TBG-TSC Hamiltonian for different values of $m$.

###  F.1 Dirac Hamiltonian with quadratic correction

Given a superconducting Hamiltonian ${\cal H}_{0}$ and an impurity potential
$V\delta(r)X$, for some matrix $X$ with eigenvalues $\pm 1$, zero modes can be
found by solving the equation

$\displaystyle({\cal H}_{0}({\bf r})$ $\displaystyle+V\delta({\bf
r})X)\psi=E\psi,$ (139) $\displaystyle\text{Det}[$ $\displaystyle
V^{-1}-G(E,{\bf r}=0)X]=0.$ (140)

where the Green’s function $G$ has been integrated over momentum to yield
$G(E,{\bf r}=0)$. A solution to Eq. (140) is indicative of Shiba states. If
there is a solution for $E$ smaller than the gap, then there are in-gap
states. Take Hamiltonian

$\displaystyle{\cal H}(\mathbf{k})$
$\displaystyle=\xi_{z}(k_{x}\tau_{z}\mu_{0}\sigma_{x}+k_{y}\tau_{0}\mu_{0}\sigma_{y})$
$\displaystyle+(\Delta-mk^{2})\xi_{x}\tau_{0}\mu_{0}\sigma_{0},$ (141)

which describes TBG-TSC.

For convenience we set $\tau_{z}=\mu_{z}=1$. The Green’s function reads

$\displaystyle G(E,\mathbf{k})=(E-{\cal H}(\mathbf{k}))^{-1}$ (142)
$\displaystyle=\dfrac{E+\xi_{z}(k_{x}\sigma_{x}+k_{y}\sigma_{y})+(\Delta-
mk^{2})\xi_{x}}{E^{2}-k_{x}^{2}-k_{y}^{2}-(\Delta-mk^{2})^{2}},$ (143)

and the terms odd in $k$ will vanish over the momentum integral required to
obtain $G(E,\mathbf{r}=0)=\int_{k}G(E,\mathbf{k})$. To prove that there are
Shiba states, we follow the procedure performed in Ref. 102.

In Ref. 102, the authors examine the Green’s function at the edges of the
positive and negative energy bands and show they may diverge to $\pm\infty$.
In certain cases one may prove the existence of Shiba states via the
intermediate value theorem: if the Green’s function diverges to $\infty$ at
the positive energy band edge and $-\infty$ at the negative energy band edge,
then for _any_ nonzero impurity strength $V$ Eq. (140) has at least one
solution in the gap. We apply their method to the Dirac Hamiltonian with
quadratic correction. At the band edges, we expect van Hove singularities,
though the exact location and nature depends upon the relative strength of
$\Delta,m$. As $k^{2}+(\Delta-
mk^{2})^{2}=\Delta^{2}+(1-2m\Delta)k^{2}+m^{2}k^{4}$, the nature of the gap
depends on the sign of $1-2m\Delta$, where we have set $v_{F}=0$. If
$m<\frac{1}{2\Delta}$, then the dispersion about $k=0$ is still parabolic. If
$m>\frac{1}{2\Delta}$, the dispersion at the band edge is actually about
$k=k_{c}$ and resembles the Mexican hat potential, see Fig. 9 333We have set
Fermi velocity $v_{F}=1$. Restoring $v_{F}$ gives the two limits for the mass
as $m>v_{F}^{2}/(2\Delta)$ and $m<v_{F}^{2}/(2\Delta)$. In this situation, we
actually expand around

$\displaystyle\Delta^{2}+(1-2m\Delta)k^{2}$
$\displaystyle+m^{2}k^{4}=m^{2}(k^{2}-k_{c}^{2})^{2}+\Delta^{\prime 2},$ (144)
$\displaystyle k_{c}^{2}$ $\displaystyle=\frac{1}{2m^{2}}(2m\Delta-1)$ (145)
$\displaystyle\Delta^{\prime}$ $\displaystyle=\frac{1}{2m}\sqrt{4m\Delta-1}.$
(146)

Figure 9: Locations where we evaluate the Green’s function. In (a), for
$m<\frac{1}{2\Delta}$, the dispersion is parabolic. In (b), for
$m<\frac{1}{2\Delta}$, the dispersion is a Mexican hat.

Notice that the only Pauli matrix that enters in final equation for $G$ is
$\xi_{x}$. This means while the impurity potential $X$ may be composed of
Pauli matrices $X=\xi_{i}\tau_{j}\mu_{k}\sigma_{l}$, the term
$\tau_{j}\mu_{k}\sigma_{l}$ will not affect the integrals and can be chosen to
be $\pm 1$. We thus only need to consider the case with $X=\xi_{0}$ and
$\xi_{x,y,z}$. The Green’s function has eigenvalues

$\displaystyle\int_{k}G(E,k)=g_{0}(E)+g_{1}(E)\xi_{x}$ (147) Eig
$\displaystyle[\int_{k}G(E,k)]=g_{0}(E)\pm g_{1}(E),$ (148) $\displaystyle
g_{0}$ $\displaystyle=\int\dfrac{kdk}{2\pi}\dfrac{E}{E^{2}-k^{2}-(\Delta-
mk^{2})^{2}}$ (149) $\displaystyle g_{1}$
$\displaystyle=\int\dfrac{kdk}{2\pi}\dfrac{\Delta-mk^{2}}{E^{2}-k^{2}-(\Delta-
mk^{2})^{2}}.$ (150)

###  F.2 $m<\frac{1}{2\Delta}$, Parabolic Dispersion

As in Ref. 102, we calculate $g_{0},g_{1}$ at the band edges. If
$m<\frac{1}{2\Delta}$ (we have chosen $\Delta>0$), the bands are parabolic and
the band edges occur at $k=0$. We evaluate the Green’s function at energy
$E=\Delta-\delta\epsilon$ for the positive energy band and
$E=-\Delta+\delta\epsilon$ for the negative energy band, $\delta\epsilon$
infinitesimal and positive. The integrals are

$\displaystyle g_{0}$
$\displaystyle=\int_{0}^{\Lambda}\dfrac{kdk}{2\pi}\dfrac{\pm\Delta\mp\delta\epsilon}{-2\Delta\delta\epsilon-(1-2m\Delta)k^{2}-m^{2}k^{4}}$
(151) $\displaystyle g_{1}$
$\displaystyle=\int_{0}^{\Lambda}\dfrac{kdk}{2\pi}\dfrac{\Delta-
mk^{2}}{-2\Delta\delta\epsilon-(1-2m\Delta)k^{2}-m^{2}k^{4}}.$ (152)

The upper bound $\Lambda$ is the UV cutoff, chosen to be much smaller than
$|\mathbf{K}|$ with $\mathbf{K}$ the $C_{3z}$-invariant K point in the
graphene Brillouin zone. We have thrown away $\delta\epsilon^{2}$ terms, as
they are small.

Separating the integrals into $R_{1},R_{2}$, and capturing leading
divergences:

$\displaystyle
R_{1}=\int_{0}^{\Lambda}\dfrac{kdk}{2\pi}\dfrac{1}{-2\Delta\delta\epsilon-(1-2m\Delta)k^{2}-m^{2}k^{4}}$
$\displaystyle\sim-\dfrac{1}{2\pi(1-2m\Delta)}\int_{0}^{0+}\dfrac{dxx}{2\Delta\delta\epsilon+x^{2}}$
$\displaystyle\sim\dfrac{1}{4\pi(1-2m\Delta)}\log{\delta\epsilon}\qquad{\rightarrow-\infty}.$
(153) $\displaystyle
R_{2}=\int_{0}^{\Lambda}\dfrac{kdk}{2\pi}\dfrac{k^{2}}{-2\Delta\delta\epsilon-(1-2m\Delta)k^{2}-m^{2}k^{4}}$
$\displaystyle=\text{finite},<0$ (154) $g_{0}=\pm\Delta R_{1}\mp\delta\epsilon
R_{1},~{}~{}g_{1}=\Delta R_{1}-mR_{2}.$ (155)

$R_{2}$ is negative as $\Delta>0,(1-2m\Delta)>0,m^{2}>0$. We finally calculate
the value of the eigenvalues of the Green’s function at the band edges:

conductance band edge: (156) $\displaystyle\qquad{g_{0}+g_{1}\sim 2\Delta
R_{1}\rightarrow-\infty}$ $\displaystyle\qquad{g_{0}-g_{1}=mR_{2}}$ valence
band edge: (157) $\displaystyle\qquad{g_{0}+g_{1}=-mR_{2}}$
$\displaystyle\qquad{g_{0}-g_{1}\sim-2\Delta R_{1}\rightarrow\infty}$

If $m>0$, the intermediate value theorem guarantees solutions to the Shiba
state equation Eq. (140) for arbitrary impurity potential $V$. If $m<0$, there
is no guarantee, consistent with topological/trivial phases of the Dirac
equation. Again, in order for $m>0$ and $m<0$ to correspond to topological and
trivial phases, we require that $|m|\Lambda^{2}\gg|\Delta|$. Without this
condition, the analysis of the Shiba states still holds, but the two cases
$m>0,m<0$ do not correspond to topological and trivial phases.

Quantity | Positive energy band edge | Negative energy band edge | Band center
---|---|---|---
Parabolic dispersion, $m<0$ |  |  |
$g_{0}+g_{1}$ | $-\infty$ | $-mR_{2}<0$ | —
$g_{0}-g_{1}$ | $mR_{2}>0$ | $\infty$ | —
$g_{0}^{2}-g_{1}^{2}$ | $-\infty$ | $-\infty$ | $-\infty$
Parabolic dispersion, $0<m<\frac{1}{2\Delta}$ |  |  |
$g_{0}+g_{1}$ | $-\infty$ | $-mR_{2}>0$ | —
$g_{0}-g_{1}$ | $mR_{2}<0$ | $\infty$ | —
$g_{0}^{2}-g_{1}^{2}$ | $\infty$ | $\infty$ | $-\infty$
Mexican hat, $m>\frac{1}{2\Delta}$ |  |  |
$g_{0}+g_{1}$ | $-\infty$ | $\infty$ | —
$g_{0}-g_{1}$ | $-\infty$ | $\infty$ | —
$g_{0}^{2}-g_{1}^{2}$ | $\infty$ | $\infty$ | $-\infty$
Table 2: Eigenvalues/determinant of the Green’s function evaluated at the band
edges and band center $E=0$. Horizontal lines indicate that the value of said
quantity is unknown, but unnecessary for general arguments.

###  F.3 $m>\dfrac{1}{2\Delta}$, Mexican Hat Dispersion

We now consider the case $m>\frac{1}{2\Delta}$. In this situation the gap
resembles the Mexican hat potential illustrated in Fig. 9. We again expand
about the conductance and valence band edges, using the formula for
$g_{0},g_{1}$ in Eq. (149), (150), now located at energy $\pm\Delta^{\prime}$
and momentum $k_{c}$:

$\displaystyle g_{0}$
$\displaystyle=\int_{0}^{\Lambda}\dfrac{kdk}{2\pi}\dfrac{\pm\Delta^{\prime}\mp\delta\epsilon}{-2\Delta^{\prime}\delta\epsilon-m^{2}(k^{2}-k_{c}^{2})^{2}}$
(158) $\displaystyle g_{1}$
$\displaystyle=\int_{0}^{\Lambda}\dfrac{kdk}{2\pi}\dfrac{\Delta-
mk^{2}}{-2\Delta^{\prime}\delta\epsilon-m^{2}(k^{2}-k_{c}^{2})^{2}}.$ (159)

with
$mk_{c}^{2}=\Delta-\frac{1}{2m},\Delta^{\prime}=\sqrt{\frac{\Delta}{m}-\frac{1}{4m^{2}}}$.
Define

$\displaystyle
R_{3}=\int_{0}^{\Lambda}\dfrac{kdk}{2\pi}\dfrac{1}{-2\Delta^{\prime}\delta\epsilon-m^{2}(k^{2}-k_{c}^{2})^{2}}$
(160)
$\displaystyle\sim\int_{k_{c}-\epsilon}^{k_{c}+\epsilon}\dfrac{kdk}{2\pi}\dfrac{1}{-2\Delta^{\prime}\delta\epsilon-m^{2}(k^{2}-k_{c}^{2})^{2}}$
(161) $\displaystyle\sim$
$\displaystyle-\dfrac{1}{2\pi}\int_{0}^{0^{+}}du\dfrac{1}{2\Delta^{\prime}\delta\epsilon+m^{2}u^{2}}$
(162) $\displaystyle\propto$ $\displaystyle-\frac{1}{2\pi
m}(2\Delta^{\prime}\delta\epsilon)^{-1/2}\qquad\rightarrow-\infty.$ (163)

The integrals become

$\displaystyle g_{0}$
$\displaystyle=(\pm\Delta^{\prime}\mp\delta\epsilon)R_{3},~{}~{}$ (164)
$\displaystyle\sim\pm\frac{1}{2m}\sqrt{4m\Delta-1}R_{3}.$ (165) $\displaystyle
g_{1}$
$\displaystyle\sim\int_{k_{c}-\epsilon}^{k_{c}+\epsilon}\dfrac{kdk}{2\pi}\dfrac{\Delta-
mk^{2}}{-2\Delta^{\prime}\delta\epsilon-m^{2}(k^{2}-k_{c}^{2})^{2}}$
$\displaystyle\approx\int_{k_{c}-\epsilon}^{k_{c}+\epsilon}\dfrac{kdk}{2\pi}\dfrac{\Delta-
mk_{c}^{2}}{-2\Delta^{\prime}\delta\epsilon-m^{2}(k^{2}-k_{c}^{2})^{2}}$ (166)
$\displaystyle=\frac{1}{2m}R_{3}.$ (167)

As $\Delta-mk_{c}^{2}=\frac{1}{2m}$,
$\Delta^{\prime}=\frac{1}{2m}\sqrt{4m\Delta-1}$. We thus have

conductance band edge: (168)
$\displaystyle{g_{0}+g_{1}\sim(\frac{1}{2m}[\sqrt{4m\Delta-1}+1])R_{3}\rightarrow-\infty}$
(169)
$\displaystyle{g_{0}-g_{1}\sim(\frac{1}{2m}[\sqrt{4m\Delta-1}-1])R_{3}\rightarrow-\infty}$
(170) valence band edge: (171)
$\displaystyle{g_{0}+g_{1}\sim(\frac{1}{2m}[-\sqrt{4m\Delta-1}+1])R_{3}\rightarrow\infty}$
(172)
$\displaystyle{g_{0}-g_{1}\sim(\frac{1}{2m}[-\sqrt{4m\Delta-1}-1])R_{3}\rightarrow\infty},$
(173)

as $m>\frac{1}{2\Delta}$, $\sqrt{4m\Delta-1}>1$.

###  F.4 Shiba States

We are now in a position to argue the existence/nonexistence of zero modes.
For impurity potential $VX=V\xi_{0}\tau_{i}\mu_{j}\sigma_{k}$, for arbitrary
$i,j,k$, the impurity violates chiral symmetry $S$ and corresponds to a time-
reversal breaking perturbation. In both scenarios of $m\Delta>0$ (either in
the parabolic or Mexican hat dispersion) the eigenvalues of the Green’s
function both go through $0$ somewhere inside the gap. Contrast this to
$m\Delta<0$, which never reaches $0$ inside the gap. Eq. (140) may still have
solutions. However, similar to Ref. 102, increasing $V$ arbitarily high will
drive the Shiba states into the bulk modes. For the case with $m\Delta>0$,
however, regardless of how strongly we increase $V$ since $G$ has a zero there
is always a solution. This is a qualitative difference between the two cases.

For impurities $X=\xi_{z}\tau_{i}\mu_{j}\sigma_{k}$, however, the situation
greatly changes. We are required to solve the equation (taking
$\tau_{i}\mu_{j}\sigma_{k}=+1)$

$\displaystyle\text{Det}[V^{-1}-G_{0}(E,{\bf r}=0)\xi_{z}]=0,$ (174)

and $G_{0}(E,{\bf r}=0)\xi_{z}$ is not Hermitian! Its eigenvalues are no
longer real, and we are not guaranteed solutions even if $G$ crosses $0$.

Solving Eq. (174) yields

$\displaystyle V^{-2}=g_{0}^{2}-g_{1}^{2},$ (175)

with $g_{0},g_{1}$ defined as before in Eq. (149), (150). To find zero modes,
observe that at $E=0$ $g_{0}=0$ (trivially by Eq. (149)), so the quantity
$g_{0}^{2}-g_{1}^{2}$ is negative at the band center $E=0$. In the case
$m\Delta<0$, $g_{0}^{2}-g_{1}^{2}=(g_{0}+g_{1})(g_{0}-g_{1})$ diverges to
$-\infty$ in both conductance and valence edges. In fact it is negative
throughout the gap, and thus there are no solutions. In the case $m\Delta>0$,
$g_{0}^{2}-g_{1}^{2}=(g_{0}+g_{1})(g_{0}-g_{1})$ diverges to $\infty$ in both
conductance and valence edges. Thus two in-gap modes are found, at $\pm
E^{\prime}$ for some $E^{\prime}$, again regardless of impurity strength. The
quantity $g_{0}^{2}-g_{1}^{2}$ is easily calculated by multiplying together
the eigenvalues $(g_{0}+g_{1}),(g_{0}-g_{1})$, calculated for the parabolic
and Mexican hat dispersions. We have collected these results in Table 2.

###  F.5 Multiple Valleys

We have demonstrated that zero modes exist for arbitrary impurity strengths
for $m\Delta>0$, while for $m\Delta<0$ Shiba states will either not exist or
can be pushed into the bulk with decreasing impurity strength. As calculated
in App. A.4, in TBG-TSC with nearly uniform pairing we have $m\Delta<0$, where
the on-site pairing strength is much greater than the NNN pairing strength,
$\bar{\Delta}_{0}\gg\bar{\Delta}_{2}$. However, in TBG-TSC the observation of
Shiba states is complicated by the presence of multiple valleys and zero
modes. A delta-function impurity itself will break valley-U(1) symmetry,
allowing for the zero modes of different valleys to hybridize.

Take the real-space perturbation

$\displaystyle{~{}~{}~{}~{}~{}~{}m\sum_{\alpha,l,s}\sum_{\mathbf{R}\in
l}\delta_{\mathbf{R},0}c^{\dagger}_{\mathbf{R},l,\alpha,s}c_{\mathbf{R},l,\alpha,s}}$
(176)
$\displaystyle={\frac{m}{N_{M}N_{0}}\sum_{\alpha,l,s,\eta,\eta^{\prime},\mathbf{k},\mathbf{k}^{\prime},\mathbf{Q},\mathbf{Q}^{\prime}}e^{-i(\eta\mathbf{K}_{l}+\mathbf{k}-\mathbf{Q})\cdot\mathbf{t}_{\alpha}}}$
$\displaystyle{\times
e^{i(\eta^{\prime}\mathbf{K}_{l}+\mathbf{k}^{\prime}-\mathbf{Q}^{\prime})\cdot\mathbf{t}_{\alpha}}c^{\dagger}_{\mathbf{k},\mathbf{Q},\eta,\alpha,s}c_{\mathbf{k}^{\prime},\mathbf{Q}^{\prime},\eta^{\prime},\alpha,s}}.$
(177)

We see that this perturbation mixes the valleys, and thus require that the
perturbation not be a delta-function impurity, but ramp up gradually over the
length scale of the graphene lattice (but quickly on the scale of the Moiré
lattice).

##  Appendix G Bosonization

We examine the fate of the gapless edge modes along the $x$-direction (Section
B.1 in the presence of strong interactions. We begin with an edge
perpendicular to the x-axis. As the dispersion is linear in free fermions and
the edge is one-dimensional, interactions can be captured via bosonization.
Under the symmetries of TBG-TSC, the fermionic movers listed in Eqs. (80)-(83)
transform as

$\displaystyle U(1)_{\text{valley}}$ $\displaystyle:\psi_{\alpha
ns}\rightarrow e^{i\gamma}\psi_{\alpha ns}$ (178) $\displaystyle
U(1)_{\text{spin}}$ $\displaystyle:\psi_{\alpha ns}\rightarrow
e^{is\gamma}\psi_{\alpha ns}$ (179) $\displaystyle T$
$\displaystyle:\psi_{\alpha ns}\rightarrow is\psi_{-\alpha-n-s}^{\dagger}$
(180) $\displaystyle i\hat{s}_{y}$ $\displaystyle:\psi_{\alpha ns}\rightarrow
s\psi_{\alpha n-s}$ (181) $\displaystyle C_{2x}$ $\displaystyle:\psi_{\alpha
ns}\rightarrow\psi_{\alpha-ns}$ (182) $\displaystyle{\cal P}_{c}$
$\displaystyle:\psi_{\alpha ns}\rightarrow i\alpha n\psi_{-\alpha
ns}^{\dagger},$ (183)

with $\alpha=1,2,n=R,L,s=\uparrow,\downarrow$, where again $\mathcal{P}_{c}$
has been defined in Eq. (A.3). There are $8$ total modes: two pairs of helical
modes per valley. Note that using the BdG basis has turned $T$ into a BdG
symmetry (Eq. (180)). For example:

$\displaystyle\psi_{1R\uparrow}(y)$
$\displaystyle=\frac{\sqrt{\lambda}}{2}\int
dxe^{-\lambda|x|}(c_{+,+,A,\uparrow}(x,y)-ic_{+,+,B,\uparrow}(x,y)-c^{\dagger}_{-,-,A,\downarrow}(x,y)-ic^{\dagger}_{-,-,B,\downarrow}(x,y))$
(184) $\displaystyle T\psi_{1R\uparrow}(y)T^{-1}$
$\displaystyle=\frac{\sqrt{\lambda}}{2}\int
dxe^{-\lambda|x|}(c_{-,-,A,\uparrow}(x,y)+ic_{-,-,B,\uparrow}(x,y)-c^{\dagger}_{+,+,A,\downarrow}(x,y)+ic^{\dagger}_{+,+,B,\downarrow}(x,y))$
$\displaystyle=i\psi^{\dagger}_{2L\downarrow}(y).$ (185)

The kinetic energy has a linear dispersion, and is simply

$\displaystyle H_{\text{kinetic}}=\sum_{\alpha ns}\int_{y}-inv_{F}\psi_{\alpha
ns}^{\dagger}\partial_{y}\psi_{\alpha ns}.$ (186)

$C_{2x}$ is preserved as it switches $y\rightarrow-y$. We first consider the
possible bilinears allowed under the above symmetries, and determine under
what conditions the edge can be fully gapped. The kinetic term has spin
rotation invariance, but when we add interactions we will not enforce full
spin-$SU(2)$ invariance, and work with the simpler $i\hat{s}_{y}$ rotation for
ease in bosonization. By spin and valley symmetry, we can add

$\displaystyle H_{\text{gap}}=\sum_{\alpha\beta mn}\int_{y}a_{\alpha\beta
mn}(\psi^{\dagger}_{\alpha m\uparrow}\psi_{\beta n\uparrow}$
$\displaystyle+\psi^{\dagger}_{\alpha m\downarrow}\psi_{\beta n\downarrow})$
$\displaystyle+H.c.,$ (187)

and the remaining symmetries further constrain the form of $a$. We further
restrict ourselves to bilinears that scatter right-movers to left-movers, as
terms mixing two right movers and two left movers do not open up gaps in the
spectrum. The effect of such scattering terms that hop between modes of the
same chirality, e.g. $\psi_{1L\uparrow}^{\dagger}\psi_{2L\uparrow}$, will only
shift the linearized spectrum to the left or right in momentum space. They
will not create gaps, so we ignore them.

$\displaystyle H_{\text{gap}}$
$\displaystyle=\sum_{\alpha\beta}\int_{y}a_{\alpha\beta}(\psi^{\dagger}_{\alpha
R\uparrow}\psi_{\beta L\uparrow}+\psi^{\dagger}_{\alpha
R\downarrow}\psi_{\beta L\downarrow})$
$\displaystyle+a_{\alpha\beta}^{*}(\psi^{\dagger}_{\beta
L\uparrow}\psi_{\alpha R\uparrow}+\psi^{\dagger}_{\beta
L\downarrow}\psi_{\alpha R\downarrow}),$ (188)

and under $T$ and ${\cal P}_{c}$:

$\displaystyle TH_{\text{gap}}T^{-1}$
$\displaystyle=\sum_{\alpha\beta}\int_{y}-a_{\alpha\beta}^{*}(\psi^{\dagger}_{-\beta
R\downarrow}\psi_{-\alpha L\downarrow}$ $\displaystyle+\psi^{\dagger}_{-\beta
R\uparrow}\psi_{-\alpha L\uparrow})+H.c.$ (189) $\displaystyle{\cal
P}_{c}H_{\text{gap}}{\cal P}_{c}^{-1}$
$\displaystyle=\sum_{\alpha\beta}\int_{y}\alpha\beta
a_{\alpha\beta}(\psi^{\dagger}_{-\beta L\uparrow}\psi_{-\alpha R\uparrow}$
$\displaystyle+\psi^{\dagger}_{-\beta L\downarrow}\psi_{-\alpha
R\downarrow})+H.c.,$ (190)

which constrains

$\displaystyle a_{\alpha\beta}=-a_{-\beta-\alpha}^{*}=\alpha\beta
a_{-\alpha-\beta}^{*}.$ (191)

The diagonal entries with $\alpha\beta=1$ vanish, and the off-diagonal entries
$a_{1,-1}=-a_{1,-1}^{*}=-a_{-1,1}^{*}$. Thus

$\displaystyle
H_{\text{gap}}=\int_{y}ia(\psi^{\dagger}_{1R\uparrow}\psi_{2L\uparrow}+\psi^{\dagger}_{2R\uparrow}\psi_{1L\uparrow})+H.c.+(\uparrow\leftrightarrow\downarrow).$
(192)

This term is odd under $C_{2x}$. Thus, if the edge preserves $C_{2x}$
symmetry, it is gapless on the level of free fermions, as we know from our
symmetry analysis in App. B.1.

To study the gapless edge in the presence of interactions, we bosonize,
following Ref. 171. Take

$\displaystyle\psi_{\alpha ns}\sim e^{i(\varphi_{\alpha s}+n\theta_{\alpha
s})},$ (193)

with bosonized fields $\varphi,\theta$ obeying

$\displaystyle[\varphi_{\alpha s}(y),\theta_{\beta
s^{\prime}}(y^{\prime})]=\frac{i\pi}{2}\delta_{\alpha\beta}\delta_{ss^{\prime}}\text{sgn}(y-y^{\prime}),$
(194)

to ensure anti-commutation between fermions at different positions. We do not
enforce translation symmetry, meaning that we will not concern ourselves with
momentum conservation when studying gap terms that may appear in the
interacting Hamiltonian.

This yields

$\displaystyle U(1)_{\text{valley}}:\varphi_{\alpha
s}\rightarrow\varphi_{\alpha s}+\gamma,~{}\theta_{\alpha
s}\rightarrow\theta_{\alpha s}$ (195) $\displaystyle
U(1)_{\text{spin}}:\varphi_{\alpha s}\rightarrow\varphi_{\alpha
s}+s\gamma,~{}\theta_{\alpha s}\rightarrow\theta_{\alpha s}$ (196)
$\displaystyle T:\varphi_{\alpha
s}\rightarrow\varphi_{-\alpha-s}-\frac{\pi}{2}s,~{}\theta_{\alpha
s}\rightarrow-\theta_{-\alpha-s}$ (197) $\displaystyle
i\hat{s}_{y}:\varphi_{\alpha
s}\rightarrow\varphi_{\alpha-s}-\frac{\pi}{2}(1-s),~{}\theta_{\alpha
s}\rightarrow\theta_{\alpha-s}$ (198) $\displaystyle C_{2x}:\varphi_{\alpha
s}(y)\rightarrow\varphi_{\alpha s}(-y),~{}\theta_{\alpha
s}(y)\rightarrow-\theta_{\alpha s}(-y)$ (199) $\displaystyle{\cal
P}_{c}:\varphi_{\alpha s}\rightarrow-\varphi_{-\alpha s},~{}\theta_{\alpha
s}\rightarrow-\theta_{-\alpha s}+\frac{\pi}{2}\alpha\text{~{}(incorrect)}.$
(200)

However, there is a subtlety regarding the transformation law for
$\mathcal{P}_{c}$. It is not consistent with the bosonization scheme: consider
a bilinear which bosonizes as

$\displaystyle e^{i\gamma}$ $\displaystyle\psi^{\dagger}_{\alpha
ms}\psi_{\beta ns}+H.c.$ $\displaystyle\sim$
$\displaystyle\sin(-\varphi_{\alpha s}-m\theta_{\alpha s}+\varphi_{\beta
s}+n\theta_{\beta s}+\gamma).$ (201)

The bilinears transform under $\mathcal{P}_{c}$ as:

$\displaystyle{\cal P}_{c}:-e^{i\gamma}$ $\displaystyle\alpha\beta
mn\psi^{\dagger}_{-\beta ns}\psi_{-\alpha ms}+H.c.\sim-\sin(-\varphi_{-\beta
s}-n\theta_{-\beta s}+\varphi_{-\alpha s}+m\theta_{-\alpha
s}+\gamma+\frac{\pi}{2}(\alpha m-\beta n)).$ (202)

However, directly transforming the bosonized cosines gives

$\displaystyle\mathcal{P}_{c}:$ $\displaystyle\sin(\varphi_{-\alpha
s}+m\theta_{-\alpha s}-\varphi_{-\beta s}-n\theta_{-\beta
s}+\gamma-\frac{\pi}{2}(m\alpha-n\beta)).$ (203)

These two results differ by an overall a minus sign. Because we concern
ourselves only with four-fermion cosines, there will always be two pairs of
fermions, so the overall sign of the quartic form is restored. However, it is
interesting to see how to resolve the apparent inconsistency in the
$\mathcal{P}_{c}$ symmetry. We resolve this issue below, by introducing
additional transformation properties that resemble Klein factors, as performed
in Ref. [172]. Define an ordering on the 8 flavors of fermion operators
$\psi_{\alpha ns}$ by

$\displaystyle f(\alpha,n,s)=j,~{}j\in\\{0,1,\cdots 7\\}.$ (204)
$\displaystyle f(1,R,\uparrow)=0,f(1,L,\uparrow)=1,$ (205) $\displaystyle
f(2,R,\uparrow)=2,f(2,L,\uparrow)=3,$ (206) $\displaystyle
f(1,R,\downarrow)=4,f(1,L,\downarrow)=5,$ (207) $\displaystyle
f(2,R,\downarrow)=6,f(2,L,\downarrow)=7.$ (208)

Define operator $\xi_{j}$ as

$\displaystyle\xi_{j}$ $\displaystyle=2\pi\sum_{j^{\prime}<j}N_{j^{\prime}},$
(209) $\displaystyle[\xi_{j},\varphi_{\alpha s}+n\theta_{\alpha s}]$
$\displaystyle=2\pi i\Theta(j-f(\alpha,n,s)-0.5)$ (210) $\displaystyle
N_{j^{\prime}}=\int_{y}\psi^{\dagger}_{\alpha ns}\psi_{\alpha ns}$
$\displaystyle,~{}j^{\prime}=f(\alpha,n,s).$ (211)

Here $\xi_{j}$ is a Jordan-Wigner type string that plays a role similar to a
Klein factor, giving anti-commutation of different flavors of fermions. This
$\xi_{j}$ will introduce an additional negative sign that fixes the
transformation rules for $\mathcal{P}_{c}$.

$\displaystyle\mathcal{P}_{c}:\varphi_{\alpha s}+n\theta_{\alpha
s}\rightarrow$ $\displaystyle-(\varphi_{-\alpha s}+n\theta_{-\alpha s})$
$\displaystyle+\frac{\pi}{2}\alpha n+\xi_{f(-\alpha,n,s)}.$ (212)

Under this new transformation, the $\mathcal{P}_{c}$ transformed cosine
becomes

$\displaystyle\mathcal{P}_{c}:$ $\displaystyle\sin(\varphi_{-\alpha
s}+m\theta_{-\alpha s}-\varphi_{-\beta s}-n\theta_{-\beta
s}+\gamma-\frac{\pi}{2}(m\alpha-n\beta)-\xi_{f(-\alpha,m,s)}+\xi_{f(-\beta,n,s)}),$
(213)

which using the BCH formula,

$\displaystyle=\frac{1}{2i}[e^{i(\varphi_{-\alpha s}+m\theta_{-\alpha
s}-\varphi_{-\beta s}-n\theta_{-\beta
s}+\gamma-\frac{\pi}{2}(m\alpha-n\beta)-\xi_{f(-\alpha,m,s)}+\xi_{f(-\beta,n,s)})}-H.c.],$
(214) $\displaystyle=\frac{1}{2i}[e^{i(\varphi_{-\alpha s}+m\theta_{-\alpha
s}-\varphi_{-\beta s}-n\theta_{-\beta
s}+\gamma-\frac{\pi}{2}(m\alpha-n\beta))}e^{i(-\xi_{f(-\alpha,m,s)}+\xi_{f(-\beta,n,s)})}$
$\displaystyle e^{\frac{1}{2}[(\varphi_{-\alpha s}+m\theta_{-\alpha
s}-\varphi_{-\beta s}-n\theta_{-\beta
s}+\gamma-\frac{\pi}{2}(m\alpha-n\beta)),(-\xi_{f(-\alpha,m,s)}+\xi_{f(-\beta,n,s)})]}-H.c.]$
(215) $\displaystyle=\frac{1}{2i}[e^{i(\varphi_{-\alpha s}+m\theta_{-\alpha
s}-\varphi_{-\beta s}-n\theta_{-\beta
s}+\gamma-\frac{\pi}{2}(m\alpha-n\beta))}$ $\displaystyle
e^{\frac{1}{2}([\varphi_{-\alpha s}+m\theta_{-\alpha
s},\xi_{f(-\beta,n,s)})]+[\varphi_{-\beta s}+n\theta_{-\beta
s},\xi_{f(-\alpha,m,s)}])}-H.c.]$ (216) $\displaystyle=-\sin(\varphi_{-\alpha
s}+m\theta_{-\alpha s}-\varphi_{-\beta s}-n\theta_{-\beta
s}+\gamma-\frac{\pi}{2}(m\alpha-n\beta))$ (217)

where in the third line we have used $e^{i\xi_{j}}=1$ (as $N_{j}$ are integer-
valued), and in the last line we have used that
$[\xi_{f(\alpha,n,s)},\varphi_{\alpha^{\prime}s^{\prime}}+n^{\prime}\theta_{\alpha^{\prime}s^{\prime}}]+[\xi_{f(\alpha^{\prime},n^{\prime},s^{\prime})},\varphi_{\alpha
s}+n\theta_{\alpha s}]=2\pi
i(\Theta(f(\alpha,n,s)-f(\alpha^{\prime},n^{\prime},s^{\prime})-0.5)+\Theta(f(\alpha^{\prime},n^{\prime},s^{\prime})-f(\alpha,n,s)-0.5))=2\pi
i$ (when $f(\alpha,n,s)\neq f(\alpha^{\prime},n^{\prime},s^{\prime})$), so the
commutator yields the required sign $-1$. Adding the $\xi$ factor fixes the
sign in the transformation law for bilinears. We will, however, only consider
4-fermion interactions (there are no cosine bilinears in our problem), so the
problematic sign fixes itself (as two negative signs cancel).

Now we add interactions. According to the rules of bosonization,

$\displaystyle n_{\alpha ms}$ $\displaystyle=\psi_{\alpha
ms}^{\dagger}\psi_{\alpha ms}-\langle\psi_{\alpha ms}^{\dagger}\psi_{\alpha
ms}\rangle$ (218) $\displaystyle=\dfrac{\partial_{y}(m\varphi_{\alpha
s}+\theta_{\alpha s})}{2\pi}.$ (219)

Beginning with the bosonized Hamiltonian for a free wire:

$\displaystyle H_{\text{kin}}=\sum_{\alpha
s}\frac{v_{F}}{2\pi}\int_{y}(\partial_{y}\varphi_{\alpha
s})^{2}+(\partial_{y}\theta_{\alpha s})^{2},$ (220)

we add symmetric density-density interactions (obeying Eq. 178-183) of the
form

$\displaystyle
u_{1}[n_{1L\uparrow}n_{1R\uparrow}+n_{1L\downarrow}n_{1R\downarrow}+n_{2L\uparrow}n_{2R\uparrow}+n_{2L\downarrow}n_{2R\downarrow}],$
$\displaystyle
u_{2}[n_{1L\uparrow}n_{1R\downarrow}+n_{1L\downarrow}n_{1R\uparrow}+n_{2L\uparrow}n_{2R\downarrow}+n_{2L\downarrow}n_{2R\uparrow}],$
$\displaystyle
u_{3}[n_{1L\uparrow}n_{1L\downarrow}+n_{1R\uparrow}n_{1R\downarrow}+n_{2L\uparrow}n_{2L\downarrow}+n_{2R\uparrow}n_{2R\downarrow}],$
$\displaystyle
u_{4}[n_{1L\uparrow}n_{2L\uparrow}+n_{1R\uparrow}n_{2R\uparrow}+n_{1L\downarrow}n_{2L\downarrow}+n_{1R\downarrow}n_{2R\downarrow}],$
$\displaystyle
u_{5}[n_{1L\uparrow}n_{2R\uparrow}+n_{1L\downarrow}n_{2R\downarrow}+n_{1R\uparrow}n_{2L\uparrow}+n_{1R\downarrow}n_{2L\downarrow}],$
$\displaystyle
u_{6}[n_{1L\uparrow}n_{2L\downarrow}+n_{1L\downarrow}n_{2L\uparrow}+n_{1R\uparrow}n_{2R\downarrow}+n_{1R\downarrow}n_{2R\uparrow}],$
$\displaystyle
u_{7}[n_{1L\uparrow}n_{2R\downarrow}+n_{1L\downarrow}n_{2R\uparrow}+n_{1R\uparrow}n_{2L\downarrow}+n_{1R\downarrow}n_{2L\uparrow}].$
(221)

These interactions renormalize the quadratic terms in bosonized language. For
example,

$\displaystyle
u_{1}[n_{1L\uparrow}n_{1R\uparrow}+n_{1L\downarrow}n_{1R\downarrow}+n_{2L\uparrow}n_{2R\uparrow}+n_{2L\downarrow}n_{2R\downarrow}]$
$\displaystyle=\frac{u_{1}}{4\pi^{2}}[(\partial_{y}(\varphi_{1\uparrow}+\theta_{1\uparrow}))(\partial_{y}(-\varphi_{1\uparrow}+\theta_{1\uparrow}))$
$\displaystyle+(1\leftrightarrow 2)+(\uparrow\leftrightarrow\downarrow)]$
$\displaystyle=\frac{u_{1}}{4\pi^{2}}[-(\partial_{y}\varphi_{1\uparrow})^{2}+(\partial_{y}\theta_{1\uparrow})^{2}+(1\leftrightarrow
2)+(\uparrow\leftrightarrow\downarrow)].$ (222)

In a similar calculation to Ref. 173, we gather all quadratic fields into
matrices $M_{\varphi,\theta}$.

$\displaystyle
H_{\text{int}}=\frac{1}{4\pi^{2}}\int_{x}\partial_{x}{\vec{\varphi}}^{T}M_{\varphi}\partial_{x}{\vec{\varphi}}+\partial_{x}{\vec{\theta}}^{T}M_{\theta}\partial_{x}{\vec{\theta}},$
(223)

with

$\displaystyle\partial_{x}{\vec{\varphi}}^{T}=\begin{bmatrix}\partial_{x}\varphi_{1\uparrow}&\partial_{x}\varphi_{2\uparrow}&\partial_{x}\varphi_{1\downarrow}&\partial_{x}\varphi_{2\downarrow}\end{bmatrix},$
(224)
$\displaystyle\partial_{x}{\vec{\theta}}^{T}=\begin{bmatrix}\partial_{x}\theta_{1\uparrow}&\partial_{x}\theta_{2\uparrow}&\partial_{x}\theta_{1\downarrow}&\partial_{x}\theta_{2\downarrow}\end{bmatrix}.$
(225)

The matrices are given by

$\displaystyle M_{\varphi/\theta}=\begin{bmatrix}\mp u_{1}&u_{4}\mp u_{5}&\mp
u_{2}+u_{3}&u_{6}\mp u_{7}\\\ u_{4}\mp u_{5}&\mp u_{1}&u_{6}\mp u_{7}&\mp
u_{2}+u_{3}\\\ \mp u_{2}+u_{3}&u_{6}\mp u_{7}&\mp u_{1}&u_{4}\mp u_{5}\\\
u_{6}\mp u_{7}&\mp u_{2}+u_{3}&u_{4}\mp u_{5}&\mp u_{1}\end{bmatrix}.$ (226)

It will be convenient for us to switch into the spin-charge basis by defining

$\displaystyle\varphi_{c}$
$\displaystyle=(\varphi_{1\uparrow}+\varphi_{2\uparrow}+\varphi_{1\downarrow}+\varphi_{2\downarrow})/2$
(227) $\displaystyle\varphi_{s}$
$\displaystyle=(\varphi_{1\uparrow}+\varphi_{2\uparrow}-\varphi_{1\downarrow}-\varphi_{2\downarrow})/2$
(228) $\displaystyle\varphi_{a}$
$\displaystyle=(\varphi_{1\uparrow}-\varphi_{2\uparrow}+\varphi_{1\downarrow}-\varphi_{2\downarrow})/2$
(229) $\displaystyle\varphi_{b}$
$\displaystyle=(\varphi_{1\uparrow}-\varphi_{2\uparrow}-\varphi_{1\downarrow}+\varphi_{2\downarrow})/2$
(230) $\displaystyle\theta_{c}$
$\displaystyle=(\theta_{1\uparrow}+\theta_{2\uparrow}+\theta_{1\downarrow}+\theta_{2\downarrow})/2$
(231) $\displaystyle\theta_{s}$
$\displaystyle=(\theta_{1\uparrow}+\theta_{2\uparrow}-\theta_{1\downarrow}-\theta_{2\downarrow})/2$
(232) $\displaystyle\theta_{a}$
$\displaystyle=(\theta_{1\uparrow}-\theta_{2\uparrow}+\theta_{1\downarrow}-\theta_{2\downarrow})/2$
(233) $\displaystyle\theta_{b}$
$\displaystyle=(\theta_{1\uparrow}-\theta_{2\uparrow}-\theta_{1\downarrow}+\theta_{2\downarrow})/2,$
(234)

which obeys the same commutation relations as the original $\varphi,\theta$
variables. This basis was obtained by diagonalizing the $M_{\varphi,\theta}$
matrix. In this new basis the Hamiltonian reads

$\displaystyle
H_{\text{kin}}+H_{\text{int}}=\sum_{i}\frac{v_{i}}{2\pi}\int_{y}K_{i}(\partial_{y}\varphi_{i})^{2}+\frac{1}{K_{i}}(\partial_{y}\theta_{i})^{2},$
(235)

with$v_{i}K_{i}$ and $v_{i}/K_{i}$ being the eigenvalues of

$\displaystyle v_{F}+2\pi M_{\varphi},v_{F}+2\pi M_{\theta}.$ (236)
$\displaystyle v_{c}$ $\displaystyle=\frac{1}{2\pi}[(2\pi
v_{F}-u_{1}-u_{2}+u_{3}+u_{4}-u_{5}+u_{6}-u_{7})$ $\displaystyle(2\pi
v_{F}+u_{1}+u_{2}+u_{3}+u_{4}+u_{5}+u_{6}+u_{7})]^{1/2},$ $\displaystyle
K_{c}$ $\displaystyle=\sqrt{\dfrac{(2\pi
v_{F}-u_{1}-u_{2}+u_{3}+u_{4}-u_{5}+u_{6}-u_{7})}{(2\pi
v_{F}+u_{1}+u_{2}+u_{3}+u_{4}+u_{5}+u_{6}+u_{7})}}$ $\displaystyle v_{s}$
$\displaystyle=\frac{1}{2\pi}[(2\pi
v_{F}-u_{1}+u_{2}-u_{3}+u_{4}-u_{5}-u_{6}+u_{7})$ $\displaystyle(2\pi
v_{F}+u_{1}-u_{2}-u_{3}+u_{4}+u_{5}-u_{6}-u_{7})]^{1/2},$ $\displaystyle
K_{s}$ $\displaystyle=\sqrt{\dfrac{(2\pi
v_{F}-u_{1}+u_{2}-u_{3}+u_{4}-u_{5}-u_{6}+u_{7})}{(2\pi
v_{F}+u_{1}-u_{2}-u_{3}+u_{4}+u_{5}-u_{6}-u_{7})}}$ $\displaystyle v_{a}$
$\displaystyle=\frac{1}{2\pi}[(2\pi
v_{F}-u_{1}-u_{2}+u_{3}-u_{4}+u_{5}-u_{6}+u_{7})$ $\displaystyle(2\pi
v_{F}+u_{1}+u_{2}+u_{3}-u_{4}-u_{5}-u_{6}-u_{7})]^{1/2},$ $\displaystyle
K_{a}$ $\displaystyle=\sqrt{\dfrac{(2\pi
v_{F}-u_{1}-u_{2}+u_{3}-u_{4}+u_{5}-u_{6}+u_{7})}{(2\pi
v_{F}+u_{1}+u_{2}+u_{3}-u_{4}-u_{5}-u_{6}-u_{7})}}$ $\displaystyle v_{b}$
$\displaystyle=\frac{1}{2\pi}[(2\pi
v_{F}-u_{1}+u_{2}-u_{3}-u_{4}+u_{5}+u_{6}-u_{7})$ $\displaystyle(2\pi
v_{F}+u_{1}-u_{2}-u_{3}-u_{4}-u_{5}+u_{6}+u_{7})]^{1/2},$ $\displaystyle
K_{b}$ $\displaystyle=\sqrt{\dfrac{(2\pi
v_{F}-u_{1}+u_{2}-u_{3}-u_{4}+u_{5}+u_{6}-u_{7})}{(2\pi
v_{F}+u_{1}-u_{2}-u_{3}-u_{4}-u_{5}+u_{6}+u_{7})}}.$ (237)

The previous symmetry-preserving density-density interactions simply
renormalized the Luttinger parameters. To fully gap the Luttinger liquid, one
must add symmetric commuting cosines involving $\varphi_{i},\theta_{i}$.
Notice that the effect of valley-U(1) symmetry is to shift
$\varphi_{c}\rightarrow\varphi_{c}+2\gamma$, with $\gamma$ the phase factor.
No other bosonized variables in Eqs. (227)-(234) are altered! Hence
$\varphi_{c}$ cannot appear in a cosine – including $\varphi_{c}$ would
violate valley charge conservation. The same holds true for $\varphi_{s}$ and
spin-U(1) conservation. Because $\varphi_{c},\varphi_{s}$ cannot appear in any
cosines, in the limit of strong pairing their conjugate variables
$\theta_{c},\theta_{s}$ may be pinned, for example via the terms:

$\displaystyle~{}w_{1}[\cos(2\theta_{s}+2\theta_{c})+\cos(2\theta_{s}-2\theta_{c})]$
$\displaystyle\propto$
$\displaystyle~{}w_{1}[\psi_{1R\uparrow}^{\dagger}\psi_{2R\uparrow}^{\dagger}\psi_{2L\uparrow}\psi_{1L\uparrow}+(\uparrow\leftrightarrow\downarrow)+H.c.]$
(238)

Note that lower harmonics of $\theta_{c},\theta_{s}$, such as
$\cos(2\theta_{s})$, cannot be expressed in terms of fermions. For example

$\displaystyle\cos(2\theta_{s})=\cos(\theta_{1\uparrow}+\theta_{2\uparrow}-\theta_{1\downarrow}-\theta_{2\downarrow}),$
(239)

while a generic hopping term is constructed out of $\psi$ operators that
bosonize as $e^{i(\varphi_{\alpha s}\pm\theta_{\alpha s})}$. This implies, for
a gap term reading

$\displaystyle\cos(\sum_{\alpha s}(m_{\alpha s}\varphi_{\alpha s}+n_{\alpha
s}\theta_{\alpha s})),$ (240)

that the parity of the $\varphi$ and $\theta$ coefficients must be equal, that
is:

$\displaystyle m_{\alpha s}=n_{\alpha s}~{}\text{mod}~{}2.$ (241)

$\cos(2\theta_{s})$ clearly violates this constraint. With the charge and spin
sectors gapped, we gap out the $a,b$ modes. There are many choices; for
example, gapping with the $\theta$ terms

$\displaystyle~{}w_{2}[\cos{(2\theta_{a}+2\theta_{b})}+\cos{(2\theta_{a}-2\theta_{b})}],$
$\displaystyle\propto$
$\displaystyle~{}w_{2}[\psi_{1R\uparrow}^{\dagger}\psi_{2L\uparrow}^{\dagger}\psi_{2R\uparrow}\psi_{1L\uparrow}+(\uparrow\leftrightarrow\downarrow)+H.c.],$
(242)

or the $\varphi$ terms,

$\displaystyle~{}w_{3}[\cos{(2\varphi_{a}+2\varphi_{b})}+\cos{(2\varphi_{a}-2\varphi_{b})}],$
$\displaystyle\propto$
$\displaystyle~{}w_{3}[\psi_{1R\uparrow}^{\dagger}\psi_{1L\uparrow}^{\dagger}\psi_{2L\uparrow}\psi_{2R\uparrow}+(\uparrow\leftrightarrow\downarrow)+H.c.],$
(243)

or a mix of both,

$\displaystyle~{}w_{4}[\cos{(2\theta_{a}+2\varphi_{b})}+\cos{(2\theta_{a}-2\varphi_{b})}],$
$\displaystyle\propto$
$\displaystyle~{}w_{4}[\psi_{1L\downarrow}^{\dagger}\psi_{2R\uparrow}^{\dagger}\psi_{2L\downarrow}\psi_{1R\uparrow}+(\uparrow\leftrightarrow\downarrow)+H.c.],$
(244)
$\displaystyle~{}w_{5}[\cos{(2\theta_{b}+2\varphi_{a})}+\cos{(2\theta_{b}-2\varphi_{a})}],$
$\displaystyle\propto$
$\displaystyle~{}w_{5}[\psi_{2L\downarrow}^{\dagger}\psi_{2R\uparrow}^{\dagger}\psi_{1R\uparrow}\psi_{1L\downarrow}+(\uparrow\leftrightarrow\downarrow)+H.c.].$
(245)

The charge and spin sectors can be gapped with the $w_{1}$ cosine, while we
have a choice of different gap terms to choose to gap out the $a,b$ sectors.
All five terms respect the symmetries of the problem, and for any given choice
of $w_{2,3,4,5}$, we can pin both the $a$ and $b$ sectors, whether by
$\varphi$ or $\theta$. The cosines commute with the $w_{1}$ cosine, and so in
principle it is possible to fully gap out the $C_{2x}$-symmetric wire.

As the Luttinger Hamiltonian is quadratic in these new variables Eq. 234, the
scaling of the cosine terms follows from standard analysis of the sine-Gordon
model, see Ref. 174.

Our RG flow equations thus read

$\displaystyle\frac{dw_{1}}{dl}$ $\displaystyle=(2-K_{s}-K_{c})w_{1}$ (246)
$\displaystyle\frac{dw_{2}}{dl}$
$\displaystyle=(2-K_{a}-K_{b})w_{2},~{}~{}\frac{dw_{3}}{dl}=(2-\frac{1}{K_{a}}-\frac{1}{K_{b}})w_{3}$
(247) $\displaystyle\frac{dw_{4}}{dl}$
$\displaystyle=(2-K_{a}-\frac{1}{K_{b}})w_{4},~{}~{}\frac{dw_{5}}{dl}=(2-\frac{1}{K_{a}}-K_{b})w_{5}.$
(248)

Regardless of the Luttinger parameters $K_{a},K_{b}$, there will always be a
choice (except for the point where $K_{a},K_{b}=1$ that makes at least one of
the terms $w_{2,3,4,5}$ relevant; thus these sectors will always be gapped.
The charge and spin sectors will be gapped so long as $K_{c},K_{s}<1$, or
repulsive interactions.

##  Appendix H Interacting zero modes

| $U(1)_{v}$ | $U(1)_{s}$ | $i\hat{s}_{y}$ | $C_{2x}$ | $T$ | ${\cal P}_{c}$
---|---|---|---|---|---|---
$\Phi_{1\uparrow}$ | $e^{i\gamma}\Phi_{1\uparrow}$ | $e^{i\gamma}\Phi_{1\uparrow}$ | $\Phi_{1\downarrow}$ | $\Phi_{2\uparrow}$ | $i\Phi_{1\downarrow}^{\dagger}$ | $-\Phi_{2\uparrow}^{\dagger}$
$\Phi_{2\uparrow}$ | $e^{i\gamma}\Phi_{2\uparrow}$ | $e^{i\gamma}\Phi_{2\uparrow}$ | $\Phi_{2\downarrow}$ | $\Phi_{1\uparrow}$ | $i\Phi_{2\downarrow}^{\dagger}$ | $\Phi_{1\uparrow}^{\dagger}$
$\Phi_{1\downarrow}$ | $e^{i\gamma}\Phi_{1\downarrow}$ | $e^{-i\gamma}\Phi_{1\downarrow}$ | $-\Phi_{1\uparrow}$ | $\Phi_{2\downarrow}$ | $-i\Phi_{1\uparrow}^{\dagger}$ | $-\Phi_{2\downarrow}^{\dagger}$
$\Phi_{2\downarrow}$ | $e^{i\gamma}\Phi_{2\downarrow}$ | $e^{-i\gamma}\Phi_{2\downarrow}$ | $-\Phi_{2\uparrow}$ | $\Phi_{1\downarrow}$ | $-i\Phi_{2\uparrow}^{\dagger}$ | $\Phi_{1\downarrow}^{\dagger}$
Table 3: Transformation laws of the zero modes $\Phi$, to be employed in the
interacting formalism.

We now consider the corner states and whether interactions can gap them out
completely in the presence of all symmetries. This is to be expected, as the
edge modes can be completely gapped, as shown in Appendix G. The zero modes
$\Phi$ transform as Table 3. For example, under $T$,

$\displaystyle\Phi_{1\uparrow}$
$\displaystyle=\frac{1}{\sqrt{\mathcal{N}}}\int dxdy\
e^{-\Delta_{0}|x|-\int^{y}_{0}dy^{\prime}M_{1}(y^{\prime})}(c_{+,+,A,\uparrow}(x,y)+c_{+,+,B,\uparrow}(x,y)-ic_{-,-,A,\downarrow}^{\dagger}(x,y)-ic_{-,-,B,\downarrow}^{\dagger}(x,y))$
(249) $\displaystyle\rightarrow\frac{1}{\sqrt{\mathcal{N}}}\int dxdy\
e^{-\Delta_{0}|x|-\int^{y}_{0}dy^{\prime}M_{1}(y^{\prime})}(c_{-,-,A,\uparrow}(x,y)+c_{-,-,B,\uparrow}(x,y)+ic_{+,+,A,\downarrow}^{\dagger}(x,y)+ic_{+,+,B,\downarrow}^{\dagger}(x,y))$
(250) $\displaystyle=i\Phi_{1\downarrow}^{\dagger}$ (251)

Under a ${\cal P}_{c}$ breaking perturbation the zero modes
$\Phi_{1\uparrow},\Phi_{1\downarrow}$ will shift away from
$\Phi_{2\uparrow},\Phi_{2\downarrow}$. Let us first work in this limit, and
ask if interactions can gap out the modes
$\Phi_{1\uparrow},\Phi_{1\downarrow}$. The density terms
$\Phi_{1\uparrow}^{\dagger}\Phi_{1\uparrow}$ and
$\Phi_{1\downarrow}^{\dagger}\Phi_{1\downarrow}$ are forbidden by
$i\hat{s}_{y}T$. We can however add the interaction

$\displaystyle(\Phi_{1\uparrow}^{\dagger}\Phi_{1\uparrow}-\frac{1}{2})(\Phi_{1\downarrow}^{\dagger}\Phi_{1\downarrow}-\frac{1}{2})$
(252)

to partially gap out the system. A twofold degeneracy still remains, and this
cannot be lifted. Note that this interaction is invariant under the full spin-
SU(2) rotation symmetry. The subtraction of $-1/2$ in the densities is
required to be consistent with time-reversal, which exchanges creation and
annihilation operators as listed in Table 3:

$\displaystyle(\Phi_{1\uparrow}^{\dagger}\Phi_{1\uparrow}-\frac{1}{2})\rightarrow(\Phi_{1\downarrow}\Phi_{1\downarrow}^{\dagger}-\frac{1}{2})=-(\Phi_{1\downarrow}^{\dagger}\Phi_{1\downarrow}-\frac{1}{2}).$
(253)

Restoring ${\cal P}_{c}$ symmetry brings all $4$ corner modes at a single
point. In class BDI the system with eight Majorana modes has been analyzed in
the presence of interactions and has been shown to be fully gapped [175].
However, our system is in class CII, with different symmetries. Density-
density interactions partially gap the system:

$\displaystyle-
r_{1}(\Phi_{1\uparrow}^{\dagger}\Phi_{1\uparrow}-\frac{1}{2})(\Phi_{1\downarrow}^{\dagger}\Phi_{1\downarrow}-\frac{1}{2})$
$\displaystyle-
r_{1}(\Phi_{2\uparrow}^{\dagger}\Phi_{2\uparrow}-\frac{1}{2})(\Phi_{2\downarrow}^{\dagger}\Phi_{2\downarrow}-\frac{1}{2}),$
(254)

which for $r_{1}>0$ leave behind four degenerate ground states

$\displaystyle\ket{1,0,1,0},~{}~{}\ket{0,1,0,1},~{}~{}\ket{1,1,1,1},~{}~{}\ket{0,0,0,0},$
(255)

where $0,1$ label the occupation of the
$\Phi_{1\uparrow},\Phi_{2\uparrow},\Phi_{1\downarrow},\Phi_{2\downarrow}$
fermions respectively. To lift this degeneracy completely we add

$\displaystyle
r_{2}\Phi_{1\uparrow}^{\dagger}\Phi_{2\uparrow}\Phi_{1\downarrow}^{\dagger}\Phi_{2\downarrow}+H.c.,$
(256)

which, with $r_{1},r_{2}>0$, leaves behind the unique ground state

$\displaystyle\ket{G.S.}=\frac{1}{\sqrt{2}}(\ket{1,0,1,0}-\ket{0,1,0,1}).$
(257)

All perturbations are symmetric, even under spin-SU(2) rotation, and as both
terms commute we have completely gapped out the zero modes. In fact, we do not
even need the density-density terms to find a unique ground state: the $r_{2}$
term is sufficient. Our results mirror the situation in class BDI: 8 Majoranas
can be fully gapped but 4 cannot [175].

##  Appendix I Josephson Junction

###  I.1 Setup

Consider a Josephson junction built out of two sheets of TBG-TSC; on the left
the superconductor is held at phase $0$, while the right the superconductor is
held at phase $\phi$. The Hamiltonian reads

$\displaystyle H[\phi]$
$\displaystyle=\xi_{z}\mu_{0}(-i\partial_{x}\tau_{z}\sigma_{x}-i\partial_{y}\sigma_{y})+\Delta\Theta(-x)\xi_{x}\mu_{0}\tau_{0}\sigma_{0}$
$\displaystyle+\Delta\Theta(x)[\cos{\phi}\xi_{x}\mu_{0}\tau_{0}\sigma_{0}+\sin{\phi}\xi_{y}\mu_{0}\tau_{0}\sigma_{0}],$
(258)

with $\Theta(x)$ the Heaviside step function. This Hamiltonian describes a
$(\phi_{L},\phi_{R})=(0,\phi)$ domain wall, where the first number of the pair
is the phase on the left and the second the phase on the right.

For $\phi=0$ the spectrum of the TBG-TSC Josephson junction is simple – the
system is one continuous sheet of TBG-TSC. There is a gap everywhere on order
of the superconducting pairing strength $\Delta$. This argument assumes that
the domain wall is perfectly sharp; eventually it will be advantageous to
relax this condition, allowing the phase to vary smoothly across the
transition. As we will show in Section I.3, this allows for the formation of
in-gap states even at $\phi=0,2\pi$.

For $\phi=\pi$ the Josephson junction reduces to the pairing domain wall; we
understand there are four _total_ zero modes for a single corner, two per
graphene valley. We wish to understand the system away from $\phi=\pi$; what
is the fate of the zero modes?

###  I.2 Symmetries of the Josephson junction

We examine which symmetries of TBG-TSC persist away from $\phi=\pi$. Chiral
symmetry $S$, time-reversal $T$, and particle-hole symmetry $\mathcal{P}$ are
broken by the inclusion of the $\Delta\xi_{y}\sin{\phi}$ term. However, they
can be restored if we take the additional operation $\phi\rightarrow-\phi$:
thus, these three symmetries conspire to relate the Hamiltonian $H[\phi]$ to
the Hamiltonian at $H[-\phi]$.

$C_{2x}$, being a mirror symmetry that flips $y\rightarrow-y$, is unaffected
by the Josephson junction (as the phase difference is along the
$x$-direction). Interestingly, $C_{2z}T$ is _also_ preserved, even when
$\phi=\pi$, when the origin of rotation is located at the junction. Even
though $T$ does not commute with the additional $\xi_{y}$ term (effectively
sending $\phi\rightarrow-\phi$), the inversion $C_{2z}$ (when centered at the
domain wall) switches the left and right sections of the domain wall. So
beginning with a domain wall with phase $0$ on the left and $\phi$ on the
right (a $(0,\phi)$ domain wall for short), $T$ will change this to a
$(0,-\phi)$ domain wall, and $C_{2z}$ flips this to $(-\phi,0)$. This is
identical to a $(0,\phi)$ domain wall via a gauge rotation
$U(\phi)=e^{-i\phi(\xi_{z}-1)/2}$. The symmetries are listed in Table 4.

Symmetry | Action on $H[\phi]$ | $\mathbf{k}\rightarrow$ | $\phi\rightarrow$
---|---|---|---
$C_{2x}$ | $\mu_{x}\sigma_{x}$ | $C_{2x}\mathbf{k}$ | $\phi$
$U(\phi)C_{2z}T$ | $e^{-i\phi(\xi_{z}-1)/2}\sigma_{x}K$ | $\mathbf{k}$ | $\phi$
$T$ | $\tau_{x}\mu_{x}K$ | $-\mathbf{k}$ | $-\phi$
${\cal P}$ | $i\xi_{z}\mu_{y}\sigma_{x}K$ | $-\mathbf{k}$ | $-\phi$
$S$ | $\xi_{y}$ | $\mathbf{k}$ | $-\phi$
Table 4: Table of symmetries in TBG-TSC Josephson junction, along with their
action on $\phi$. The latter $3$ symmetries $T,\mathcal{P},S$ are not
symmetries of $H[\phi]$; they relate $H[\phi]$ to $H[-\phi]$.

Though $C_{2z}T$ is present at $\phi\neq\pi$, chiral symmetry $S$ is not, and
both were required to enforce the four zero modes at $0$ energy. Thus, the
modes will move away from $0$. The exact behavior of the four in-gap states,
which we call $\Phi_{as}[\phi]$, can be deduced from general arguments:

  * •

1\. Because of spin rotation $i\hat{s}_{y}$ symmetry, the four in-gap states
can be thought of as two pairs $\Phi_{as}[\phi]$, with
$\Phi_{1\uparrow}[\phi]$ and $\Phi_{1\downarrow}[\phi]$ having identical
spectra, and the same for the
$\Phi_{2\uparrow}[\phi],\Phi_{2\downarrow}[\phi]$ pair. When $\mathcal{P}_{c}$
is broken, these two pairs of in-gap states separate in position space, with
one pair moving up and one moving down.

  * •

2\. Even if $\mathcal{P}_{c}$ is broken, the two pairs of corner in-gap states
are related by $C_{2x}$ and _also_ have identical spectra. If $C_{2x}$ is
broken then this is no longer true, but for weak $C_{2x}$-breaking the
qualitative nature of the spectrum will not change dramatically.

  * •

3\. At $\phi=\pi$, the in-gap states are zero modes
$\Phi_{as}[\phi=\pi]=\Phi_{as}$, and at $\phi=0,2\pi$, they are fully merged
with the bulk states. This will no longer be true when we consider the soft
domain wall; in this situation at $\phi=0,2\pi$ there are in-gap corner
states.

  * •

4\. Because of $S$ symmetry relating the Hamiltonians at $\phi$ and $-\phi$,
the spectrum at $\phi,E[\phi]$, is related to the opposite spectrum by
$E[\phi]=-E[-\phi]$.

We thus deduce the form of the in-gap state spectra to be as depicted in Fig.
10: there are four in-gap modes corresponding to four zero modes at
$\phi=\pi$: two $(\mu_{z}=\pm 1)$ per graphene valley. When $\mathcal{P}_{c}$
is broken, the two pairs of in-gap states will separate in space; one pair
moves up and the other moves down in real space along the domain wall.
Nevertheless, the energy spectra of the in-gap states will not change: they
are forced to all be identical to one another by either $C_{2x}$ or spin-U(1).
For the rest of this Appendix, we assume that the two pairs of in-gap states
are separated and only consider a single pair.

###  I.3 Soft Domain Wall

Figure 10: A soft potential at the domain wall creates in-gap states. For a
hard edge, when the superconducting phase difference is $0$, if the sheet of
TBG-TSC is homogeneous (as depicted in (a)), the system is merely a continuous
sheet of TBG-TSC and the in-gap states (depicted in (c)) flow into the bulk.
Spin-rotation symmetry keeps the in-gap states corresponding to
$\uparrow,\downarrow$ degenerate. We model a soft potential in (b), where the
left and right superconductors are separated by a pairing region at
superconducting phase $\pi$. Corner modes are bound to both domain walls, and
can hybridize as the boundary is shrunk, as depicted in (d). The in-gap energy
states are now well-separated from the bulk. The gap arising between the
dashed line stems from hybridization in the dashed line in (b), and can be in
principle smaller than the gap $\Delta$.

Our previous calculations have assumed that the transition between the two
superconductors forming the Josephson junction is a sharp step, while in
reality the transition will be smooth. In fact, in order to preserve
valley-U(1), we actually require that the phase jump be smooth on the order of
the graphene lattice (but sharp on the order of the Moiré lattice).

When the domain wall separating the superconductors is spread out instead of a
sharp step, additional in-gap states are pulled from the continuum. An easy
way to see this is to pick a judicious form of the domain wall: consider a
$(0,\pi,0)$ junction, where there are _three_ superconducting regions at phase
$\phi=0,\pi,0$ respectively, for the left, middle, and right regions, as
depicted in Fig. 10(b). The middle region is new and is short on the order of
the Moiré length scale. There are thus _two_ domain walls, separated by a
small region. Because the left and the right regions are both at phase
$\phi=0$, this is a Josephson junction with phase difference $0$; a $(0,0)$
domain wall. The purpose of the middle $\phi=\pi$ region is to mimic a soft
domain wall, with the great advantage that we already know the low-lying
states that exist at the $(0,\pi)$ and $(\pi,0)$ domain walls: call them
$\Phi_{as},\Phi_{as}^{\prime}$ respectively. These are allowed to hybridize
because the center region is narrow; if this hybridization is small then there
are in-gap states even at $\phi=0,2\pi$.

To calculate the zero modes and their properties under symmetries at the
$(\pi,0)$ domain wall, we perform the gauge transformation
$U(\pi)=e^{-i\pi(\xi_{z}-1)/2}=\xi_{z}$ that takes $c\rightarrow ic$. This has
no effect on the hopping terms $c^{\dagger}c$, but it reverses the sign of the
pairing terms $\Delta cc\rightarrow-\Delta cc$. In the Dirac Hamiltonian
notation, $\xi_{z}$ commutes with the kinetic term but anti-commutes with the
pairing. The gauge transformation thus maps a $(0,\pi)$ junction to a
$(\pi,0)$ junction. The primed zero mode operators $\Phi^{\prime}$ read

$\Phi_{1\uparrow}^{\prime}=\frac{1}{\sqrt{\mathcal{N}}}\int dxdy\
e^{-\Delta_{0}|x|-\int^{y}_{0}dy^{\prime}M_{1}(y^{\prime})}(ic_{+,+,A,\uparrow}(x,y)+ic_{+,+,B,\uparrow}(x,y)-c_{-,-,A,\downarrow}^{\dagger}(x,y)-c_{-,-,B,\downarrow}^{\dagger}(x,y)),$
(259) $\Phi_{2\uparrow}^{\prime}=\frac{1}{\sqrt{\mathcal{N}}}\int dxdy\
e^{-\Delta_{0}|x|-\int^{y}_{0}dy^{\prime}M_{1}(y^{\prime})}(ic_{+,-,A,\uparrow}(x,y)+ic_{+,-,B,\uparrow}(x,y)-c_{-,+,A,\downarrow}^{\dagger}(x,y)-c_{-,+,B,\downarrow}^{\dagger}(x,y)),$
(260)

and under $T$ transform as
$\Phi_{1\uparrow}^{\prime}\rightarrow-i\Phi_{1\downarrow}^{\prime\dagger},\Phi_{2\uparrow}^{\prime}\rightarrow-i\Phi_{2\downarrow}^{\prime\dagger}.$
We are allowed to hybridize corner states across the soft domain wall as

$\displaystyle
H_{\text{domain}}=v_{0}\Phi_{2\uparrow}^{\dagger}\Phi_{2\uparrow}^{\prime}+v_{0}\Phi_{2\downarrow}^{\dagger}\Phi_{2\downarrow}^{\prime}+H.c.,$
(261)

with four eigenstates, at energies $E=\pm|v_{0}|$. Without loss of generality,
we take $v_{0}$ real and the low-energy eigenstates as

$\displaystyle\Psi_{a\uparrow}=\frac{1}{\sqrt{2}}(\Phi_{2\uparrow}+\Phi_{2\uparrow}^{\prime}),$
(262)
$\displaystyle\Psi_{b\uparrow}=\frac{1}{\sqrt{2}}(\Phi_{2\uparrow}-\Phi_{2\uparrow}^{\prime}).$
(263)

Under $T$,

$\displaystyle\Psi_{a\uparrow}\rightarrow-i\Psi_{b\downarrow}^{\dagger}$ (264)
$\displaystyle\Psi_{b\uparrow}\rightarrow-i\Psi_{a\downarrow}^{\dagger}.$
(265)

$\displaystyle
H_{\text{domain}}=v_{0}(\Psi_{a\uparrow}^{\dagger}\Psi_{a\uparrow}-\Psi_{b\uparrow}^{\dagger}\Psi_{b\uparrow})+(\uparrow\leftrightarrow\downarrow).$
(266)

Thus, if the hybridization is small $v_{0}<|\Delta|$, then there are in-gap
states at energy $\pm v_{0}$.

###  I.4 Interactions

Figure 11: Single-particle and multi-particle spectra of the fractional
Josephson effect in TBG-TSC. (a) Single-particle spectra. Multiple in-gap
states have been created by modifying the potential. The spectrum is doubled
for both spins. If $\mathcal{P}$ is preserved, the spectrum is quadrupled. (b)
Many-body spectrum, $\mathcal{P}$ broken. The many-body states at $\phi=0$ are
labeled by occupation of $\Psi_{-2s},\Psi_{-1s},\Psi_{1s},\Psi_{2s}$
respectively. The subscript $\uparrow\downarrow$ indicates the exact same
occupation in both spin sectors. At $\phi=\pi$, the occupation is labeled by
$\Phi_{-1s},\Phi_{0s},\Phi_{1s}$, respectively. The crossing in the shaded
circle can be avoided by the term
$\Phi_{-1\uparrow}\Phi_{0\uparrow}\Phi_{1\uparrow}\Phi_{-1\downarrow}\Phi_{0\downarrow}\Phi_{1\downarrow}+H.c.$.
There are many other states in the multi-particle spectra, but their presence
does not affect the pumping cycle qualitatively.

Consider now a pumping cycle that winds $\phi$ by $2\pi$. Beginning in the
ground state at $\phi=0$, what is the fate of the system after $\phi$ is
adiabatically changed by $2\pi$? First, let us pull down additional in-gap
states so that at $\phi=0$, there are $4$ in-gap fermion states per spin
sector, or $8$ total; there are an additional $8$ modes at the
$C_{2x}$-symmetric position that we do not consider. The eight fermions at
$\phi=0$ are $\Psi_{b\uparrow},\Psi_{b\downarrow}$, where $b$ labels the
flavor of $\Psi$ (_not_ the corners separated in space!). At $\phi=\pi$ the
six in-gap states are $\Phi_{b\uparrow},\Phi_{b\downarrow}$; again, $b$ refers
to the flavor of in-gap state, not the different $C_{2x}$-related corners.
This is depicted in Fig. 11. Without interactions, we take the non-interacting
Hamiltonian

$\displaystyle H_{\text{in-
gap}}[\phi=0]=\sum_{s}E_{1}(\Psi_{1s}^{\dagger}\Psi_{1s}-\frac{1}{2})+E_{2}(\Psi_{2s}^{\dagger}\Psi_{2s}-\frac{1}{2})$
$\displaystyle\qquad-
E_{1}(\Psi_{-1s}^{\dagger}\Psi_{-1s}-\frac{1}{2})-E_{2}(\Psi_{-2s}^{\dagger}\Psi_{-2s}-\frac{1}{2}),$
(267) $\displaystyle H_{\text{in-
gap}}[\phi=\pi]=\sum_{s}E_{1}^{\prime}(\Phi_{1s}^{\dagger}\Phi_{1s}-\frac{1}{2})-E_{1}^{\prime}(\Phi_{-1s}^{\dagger}\Phi_{-1s}-\frac{1}{2}).$
(268)

The hybridization energies satisfy $E_{1},E_{2},E_{1}^{\prime}<|\Delta|$. This
spectrum is depicted in Fig. 11. Following the single-particle spectrum, the
ground state initially consists of all negative energy states unoccupied.
After evolution, the states $\Psi_{1s}$ go from unoccupied to occupied, and
the ground state is not mapped back to its original form. An additional
evolution of $\phi$ by $2\pi$ will then turn $\Psi_{2s}$ from unoccupied to
occupied. Each evolution of $2\pi$, the valley number of the state changes by
$+2$, as states are dragged from the negative energy continuum and pumped into
the positive energy continuum. So long as valley-U(1) is preserved, there is
no way for the ground state to return to its original form. The Josephson
cycle is aperiodic. To remedy this issue, we must consider interactions, which
will serve to reduce the valley-U(1) symmetry down to $Z_{3}$.

The valley-U(1) symmetry originates from the fact that the twist angle between
the bilayers of graphene is small. This means that the Moiré potential is
smooth on the scale of the graphene lattice vectors, and the graphene
translation operators $T_{a_{i}}$ are approximately good operators. In
momentum space, the low-energy excitations are near the Dirac points
$K,K^{\prime}$, meaning that one requires a large number of Moiré BZ vectors
to scatter a state in one valley to the other. In other words,

$\displaystyle T_{a_{i}}\psi_{+}\approx e^{-iK\cdot a_{i}}\psi_{+},$ (269)
$\displaystyle T_{a_{i}}\psi_{-}\approx e^{iK\cdot a_{i}}\psi_{-}.$ (270)

The states around the positive valley $\psi_{+}$ are located at $K$, while the
opposite valley at $K^{\prime}=-K$. Thus hoppings between the two valleys are
forbidden. At the free fermion level, we thus have independent number
conservation of both the $+$ valley fermions and $-$ valley fermions. Note
that pairing between opposite valleys is allowed, and is the form of the
pairing Hamiltonian we take.

Once we turn to the multi-particle language, more terms are allowed. Notice
that $3K=0$, that is, three times the $K$ point yields a reciprocal lattice
vector in the graphene BZ. This means that, despite not conserving
valley-U(1), terms of the form
$\psi_{1+}\psi_{2+}\psi_{3+}\psi_{4+}\psi_{5+}\psi_{6+}$ preserve translation
symmetry in the graphene lattice. We thus claim that in the many-body
formalism, valley-U(1) symmetry reduces to $Z_{3}$, with an extra $Z_{2}$ from
fermion parity yielding $Z_{6}$.

With valley-U(1) broken to $Z_{3}$, the Josephson effect is now $6\pi$
periodic. The multi-particle ground state evolution is depicted in Fig. 11(b).
In the multi-particle spectrum, all crossings are protected by $Z_{3}$
symmetry, except for the shaded circle, which can be gapped by

$\displaystyle
u\Phi_{-1\uparrow}\Phi_{0\uparrow}\Phi_{1\uparrow}\Phi_{-1\downarrow}\Phi_{0\downarrow}\Phi_{1\downarrow}+H.c.,$
(271)

which preserves $Z_{3}$ symmetry. The $\Phi_{0s}$ operator is the zero mode at
$\phi=\pi$ and the $\Phi_{1s},\Phi_{-1s}$ are two in-gap states, as depicted
in Fig. 11(a). Tracing out the evolution of the ground state reveals a $6\pi$
periodicity, or $Z_{3}$ fractional Josephson effect.

We explain the evolution of the ground state explicitly. Consider the two spin
sectors to be identical copies that behave identically under evolution of
$\phi$: in the spin up sector, $\ket{1,1,0,0}_{\uparrow}$ is the ground state,
the numbers denoting the occupation of
$\Psi_{-2\uparrow},\Psi_{-1\uparrow},\Psi_{1\uparrow},\Psi_{2\uparrow}$. After
$2\pi$ evolution, this state flows to $\ket{1,1,1,0}_{\uparrow}$, and then to
$\ket{1,1,1,1}_{\uparrow}$, and then into the continuum. In order to prevent
this, we must open a gap between two states of different valley number, which
we have argued can only occur when the valley number differs by a multiple of
$3$.

To see the cycle explicitly, define the ground state at $\phi=0$ as
$\ket{\Omega(\phi=0)}$, which satisfies the following:

$\displaystyle\Psi_{-1s}^{\dagger}\ket{\Omega(\phi=0)}=0$ (272)
$\displaystyle\Psi_{-2s}^{\dagger}\ket{\Omega(\phi=0)}=0$ (273)
$\displaystyle\Psi_{1s}\ket{\Omega(\phi=0)}=0$ (274)
$\displaystyle\Psi_{2s}\ket{\Omega(\phi=0)}=0,$ (275)

and at $\phi=\pi$, the two ground states
$\ket{\Omega(\phi=\pi)},\Phi_{0s}^{\dagger}\ket{\Omega(\phi=\pi)}$ satisfy

$\displaystyle\Phi_{0s}\ket{\Omega(\phi=\pi)}=0$ (276)
$\displaystyle\Phi_{1s}\ket{\Omega(\phi=\pi)}=0$ (277)
$\displaystyle\Phi_{-1s}^{\dagger}\ket{\Omega(\phi=\pi)}=0.$ (278)

We now study the flow of the ground state under winding $\phi$. The ground
state evolves as

$\displaystyle\ket{\Omega(\phi=0)}$
$\displaystyle\rightarrow\Phi_{0\uparrow}^{\dagger}\Phi_{0\downarrow}^{\dagger}\ket{\Omega(\phi=\pi)}\rightarrow\Psi_{1\uparrow}^{\dagger}\Psi_{1\downarrow}^{\dagger}\ket{\Omega(\phi=0)}$
$\displaystyle\rightarrow\Phi_{1\uparrow}^{\dagger}\Phi_{0\uparrow}^{\dagger}\Phi_{1\downarrow}^{\dagger}\Phi_{0\downarrow}^{\dagger}\ket{\Omega(\phi=\pi)},$
(279)

or in the occupation number notation

$\displaystyle\ket{1,1,0,0}_{\uparrow,\downarrow}\rightarrow\ket{1,1,0}_{\uparrow,\downarrow}\rightarrow\ket{1,1,1,0}_{\uparrow,\downarrow}\rightarrow\ket{1,1,1}_{\uparrow,\downarrow}.$
(280)

But now with interactions, the state
$\Phi_{1\uparrow}^{\dagger}\Phi_{0\uparrow}^{\dagger}\Phi_{1\downarrow}^{\dagger}\Phi_{0\downarrow}^{\dagger}\ket{\Omega(\phi=\pi)}$
can tunnel to $\Phi_{-1\uparrow}\Phi_{-1\downarrow}\ket{\Omega(\phi=\pi)}$, as
they differ by $6$ operators: the perturbation
$u\Phi_{-1\uparrow}\Phi_{0\uparrow}\Phi_{1\uparrow}\Phi_{-1\downarrow}\Phi_{0\downarrow}\Phi_{1\downarrow}+H.c.$
serves this purpose. Continuing the evolution gives

$\displaystyle\Phi_{-1\uparrow}\Phi_{-1\downarrow}\ket{\Omega(\phi=\pi)}$
$\displaystyle\rightarrow\Psi_{-1\uparrow}\Psi_{-1\downarrow}\ket{\Omega(\phi=0)}$
$\displaystyle\rightarrow\ket{\Omega(\phi=\pi)}\rightarrow\ket{\Omega(\phi=0)},$
(281)

or in the occupation number notation,

$\displaystyle\ket{0,0,0}_{\uparrow,\downarrow}\rightarrow\ket{1,0,0,0}_{\uparrow,\downarrow}\rightarrow\ket{1,0,0}_{\uparrow,\downarrow}\rightarrow\ket{1,1,0,0}_{\uparrow,\downarrow}.$
(282)

This is the $Z_{3}$ fractional Josephson effect.

Two important considerations needed to protect the $Z_{3}$ fractional
Josephson effect are that the states
$\ket{1,1,1,0}_{\uparrow\downarrow},\ket{1,0,0,0}_{\uparrow\downarrow}$ have a
protected crossing, as well as the states
$\ket{1,1,0}_{\uparrow\downarrow},\ket{1,0,0}_{\uparrow\downarrow}$. Both
crossings are protected by valley symmetry – their total valley numbers do not
differ by a multiple of $3$.