††thanks: These two authors contributed equally††thanks: These two authors
contributed equally

# van Hove, Rashba, and Hubbard meet to form
first-order and higher-order topological superconductors

Pietro M. Bonetti  Max Planck Institute for Solid State Research,
Heisenbergstrasse 1, D-70569 Stuttgart, Germany  Debmalya Chakraborty
Department of Physics and Astronomy, Uppsala University, Box 516, S-751 20
Uppsala, Sweden  Xianxin Wu  Institute for Theoretical Physics, Chinese
Academy of Sciences, Beijing, China  Andreas P. Schnyder  Max Planck Institute
for Solid State Research, Heisenbergstrasse 1, D-70569 Stuttgart, Germany
Yukawa Institute for Theoretical Physics (YITP), Kyoto University, Kyoto
606-8502, Japan

(August 27, 2024)

###### Abstract

We investigate topological superconductivity in the Rashba-Hubbard model,
describing heavy-atom superlattice and van der Waals materials with broken
inversion. We focus in particular on fillings close to the van Hove
singularities, where a large density of states enhances the superconducting
transition temperature. To determine the topology of the superconducting gaps
and to analyze the stability of their surface states in the presence of
disorder and residual interactions, we develop an fRG+MFT approach, which
combines the unbiased functional renormalization group (fRG) with a real-space
mean-field theory (MFT). Our approach uncovers a cascade of topological
superconducting states, including $A_{1}$ and $B_{1}$ pairings, whose wave
functions are of dominant $p$\- and $d$-wave character, respectively, as well
as a time-reversal breaking $A_{1}+iB_{1}$ pairing. While the $A_{1}$ and
$B_{1}$ states have first order topology with helical and flat-band Majorana
states, respectively, the $A_{1}+iB_{1}$ pairing exhibits second-order
topology with Majorana corner modes. Implications for experimental systems and
potential applications for quantum technologies are being discussed.

Topological superconductors (TSCs) are of high current interest due to their
exceptional properties and potential for applications in quantum information
technologies [1, 2, 3, 4, 5, 6, 7]. There are two different types of TSCs: (i)
Engineered heterostructures in which topological superconductivity is
proximity induced by a conventional $s$-wave superconductor and (ii) intrinsic
TSCs which arise as an intrinsic property of the material. While a great
variety of different heterostructures realizing TSCs has been proposed [7],
progress in the field of intrinsic TSCs has been less swift. One of the most
promising recent examples of intrinsic TSCs are perhaps the iron-based
superconductors [8, 9], in which evidence for zero-energy vortex bound states
of topological origin has been reported [10, 11]. Intrinsic TSCs have several
advantages compared to their engineered counterparts: For one, they exhibit
higher transition temperatures and are less affected by impurities. Moreover,
their topological superconductivity exists in the entire volume of the
material, and not just at a heterostructure interface.

With the recent progress in the MBE synthesis of superlattices [12] and in the
fabrication of 2D van der Waals materials [13, 14, 15], the design of
intrinsic TSCs is now becoming a reality. Particularly promising are van der
Waals heavy-atom materials [16, 17, 18] and heavy-atom superlattices, such as
LAO/STO [19, 20], EuO/KTO [21], and CeCoIn5/YbCoIn5 [22, 23]. The latter shows
signs of unconventional superconductivity below the relatively high transition
temperature $T_{c}\simeq 2$ K. An important feature of these designer
materials is their high tunability. For example, the Rashba spin-orbit
coupling in the CeCoIn5/YbCoIn5 superlattices can be adjusted by the width of
the YbCoIn5 blocks. In light of these recent developments, it is now becoming
paramount to investigate which types of Fermi surface shapes and spin-orbit
couplings are conducive for topological superconductivity at large $T_{c}$.

Motivated by this, we theoretically analyze the conditions under which
topological superconductivity emerges in the Rashba-Hubbard model [24, 25, 26,
27, 28, 29, 30, 31], which captures the electronic correlation physics of the
above mentioned 2D heavy-atom materials. Focusing on the $t$-$t^{\prime}$
square lattice with nearest-neighbor Rashba spin-orbit coupling, we
investigate the superconducting pairing propensity as a function of
$t^{\prime}$ and hole doping. We are in particular interested in the pairing
states near the van Hove fillings, since there a large density of states
enhances $T_{c}$. Due to the mutual interference between particle-particle and
particle-hole instabilities at the van Hove singularities, an unbiased
numerical method is required to reliably determine the superconducting pairing
states. We therefore employ the truncated-unity functional renormalization
group (fRG) [32, 33], which treats all instability channels on an equal
footing. In order to investigate the topological characteristics as well as
the edge and disorder properties of the dominant pairing states, we introduce
here, in addition, a novel fRG+MFT approach, which combines the fRG analysis
with a real-space mean-field theory (MFT) calculation [34, 35, 36, 37]. Within
this framework, we are able to deduce in an impartial manner not only the
pairing symmetries, but also the topological and edge properties of the
superconducting states, as well as the (in-)stabilities of the edge states
against residual interactions. This allows us to quantitatively predict the
topological superconducting features of the Rashba-Hubbard model, which in
turn informs the design of superconducting 2D van der Waals and superlattice
materials with nontrivial topology.

Specifically, we find that ferro- and antiferromagnetic fluctuations, which
are strongly enhanced near the van Hove singularities, lead to $A_{1}$ and
$B_{1}$ superconducting states, with dominant $p$\- and $d$-wave like pairing
character, respectively, and a $T_{c}$ estimated to be as high as $\sim
10^{-2}t$ for $U=3t$. Both of these pairing states exhibit nontrivial first
order topology [3] and host helical and flat-band Majorana edge states,
respectively. Remarkably, we find near the phase boundary between these two
states in addition a time-reversal breaking $A_{1}+iB_{1}$ pairing state with
higher-order topology and Majorana corner modes [38, 39, 40, 41, 42, 43, 44,
45, 46]. Our combined fRG+MFT approach reveals that while both the helical
Majorana and corner edge states are robust to residual interactions, the flat-
band Majorana states are unstable towards the formation of a 1D phase crystal
[47, 48, 49]. Our findings provide a route towards the design of intrinsic
topological superconductors with potential applications in quantum information
technologies [50, 51].

Results
Before analyzing the different pairing instabilities in detail, we first
briefly describe the model Hamiltonian and the employed fRG + MFT approach.

Model and methods. We start from the Rashba-Hubbard Hamiltonian on the square
lattice (with lattice spacing $a=1$) given by $H=H_{0}+H_{U}$, with

$\begin{split}&H_{0}=\sum_{j,j^{\prime},\sigma}t_{jj^{\prime}}c^{\dagger}_{j,\sigma}c_{j^{\prime},\sigma}-\mu\sum_{j,\sigma}n_{j,\sigma}\\\
&+i\sum_{\begin{subarray}{c}j,j^{\prime}\\\
\sigma,\sigma^{\prime}\end{subarray}}\lambda_{jj^{\prime}}\left[\left(\mathbf{r}_{j}-\mathbf{r}_{j^{\prime}}\right)\times
c^{\dagger}_{j,\sigma}\vec{\tau}_{\sigma\sigma^{\prime}}c_{j^{\prime},\sigma^{\prime}}\right]_{z}+\textrm{H.c.},\\\
&H_{U}=U\sum_{j}n_{j,\uparrow}n_{j,\downarrow},\end{split}$ (1)

where $c_{j,\sigma}$ ($c_{j,\sigma}^{\dagger}$) is the annihilation (creation)
operator of an electron at site $j$ with spin projection $\sigma$,
$\vec{\tau}$ are the Pauli matrices, $\mathbf{r}_{j}$ is the lattice
coordinate of site $j$, $n_{j,\sigma}=c^{\dagger}_{j,\sigma}c_{j,\sigma}$ is
the spin-resolved particle number operator, $\mu$ is the chemical potential,
$H_{0}$ is the non-interacting part and $H_{U}$ is the Hubbard interaction
term. Hamiltonian (1) is described by few parameters: the hopping amplitudes
$t_{jj^{\prime}}$, the Rashba couplings $\lambda_{jj^{\prime}}$, and the
onsite Hubbard repulsive interaction $U>0$. In the following, we choose the
hopping amplitudes such that $t_{\langle j,j^{\prime}\rangle}=-t$ when $j$ and
$j^{\prime}$ are nearest neighbors, $t_{\langle\langle
j,j^{\prime}\rangle\rangle}=-t^{\prime}$ when $j$ and $j^{\prime}$ are second
neighbors, and zero otherwise. We also consider the Rashba spin-orbit coupling
to be nonzero only for nearest neighbors $\lambda_{\langle
j,j^{\prime}\rangle}=\lambda t$. Fourier transforming the quadratic part of
the Hamiltonian, one obtains

$H_{0}=\int_{\mathbf{k}}c^{\dagger}_{\mathbf{k}}\hat{\mathcal{H}}_{\mathbf{k}}c_{\mathbf{k}},$
(2)

where $\int_{\mathbf{k}}=\int_{\mathrm{BZ}}\frac{d^{2}k}{(2\pi)^{2}}$ is a
shorthand for the integral over the first Brillouin zone.
$\hat{\mathcal{H}}_{\mathbf{k}}$ is a $2\times 2$ matrix that reads as

$\hat{\mathcal{H}}_{\mathbf{k}}=(\epsilon_{\mathbf{k}}-\mu)\mathbb{1}+\lambda\vec{g}_{\mathbf{k}}\cdot\vec{\tau},$
(3)

with $\epsilon_{\mathbf{k}}=-2t(\cos k_{x}+\cos k_{y})-4t^{\prime}\cos
k_{x}\cos k_{y}$ the Fourier transform of $t_{jj^{\prime}}$, and
$\vec{g}_{\mathbf{k}}=2t(-\sin k_{y},\sin k_{x},0)$. Hamiltonian (3) can be
easily diagonalized introducing the so-called _helicity basis_ ,
$c_{\mathbf{k}}\to U_{\mathbf{k}}c_{\mathbf{k}}$, with

$\displaystyle
U_{\mathbf{k}}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{cc}1&e^{-i\phi_{\mathbf{k}}}\\\
-e^{i\phi_{\mathbf{k}}}&1\end{array}\right),$ (6)

with
$e^{i\phi_{\mathbf{k}}}=(g_{1,\mathbf{k}}+ig_{2,\mathbf{k}})/|\vec{g}_{\mathbf{k}}|$,
leading to the quasiparticle bands
$E^{\pm}_{\mathbf{k}}=\epsilon_{\mathbf{k}}-\mu\pm\lambda|\vec{g}_{\mathbf{k}}|$.
Model (1) explicitly breaks parity (inversion) symmetry, hence its
superconducting states cannot be classified as even-parity (spin-singlet) or
odd-parity (spin-triplet). Rather, they are a superposition of the two.

Figure 1:  Superconducting phase diagram. Phase diagram of the Rashba-Hubbard
model on the square lattice as a function of chemical potential $\mu$ and
second-neighbor hopping $t^{\prime}$ for $U=3t$. $t^{\prime}$ and $\mu$ are
measured in units of the nearest-neighbor hopping amplitude $t$. The dashed
lines indicate the values of the chemical potential
$\mu_{\mathrm{vH_{1}}}(t^{\prime})$, $\mu_{\mathrm{vH_{2}}}(t^{\prime})$ where
the van Hove singularities VHS1 and VHS2 occur.  Figure 2:  Triplet-to-singlet
ratio and ferromagnetic fluctuations. Panel (a): triplet strength $\eta$ of
the superconducting state as a function of chemical potential $\mu$ and
second-neighbor hopping $t^{\prime}$. Panel (b): leading superconducting
eigenvalue $\lambda_{\mathrm{SC}}$ computed at $T=10^{-6}t$. Panels (c) and
(d): ferromagnetic fluctuations in (c) the transverse and (d) the longitudinal
spin susceptibility. The insets are colormaps of the momentum dependence of
the relative susceptibility.

We investigate the superconducting (SC) phases of the model in Eq. (1) by
means of the functional renormalization (fRG) group, combined with mean-field
theory (fRG+MFT) [34, 35, 36, 37]. This results in a renormalized
superconducting gap equation of the form

$\Delta^{\sigma\sigma^{\prime}}_{\mathbf{k}}=\sum_{s,s^{\prime}}\int_{\mathbf{k}^{\prime}}\widetilde{V}^{\sigma\sigma^{\prime}ss^{\prime}}_{\mathbf{k}\mathbf{k}^{\prime}}\left[T\sum_{n}F^{ss^{\prime}}(\mathbf{k}^{\prime},\nu_{n})\right],$
(7)

where $\nu_{n}$ are fermionic Matsubara frequencies, $T$ is the temperature,
$\Delta^{\sigma\sigma^{\prime}}_{\mathbf{k}}$ is the SC order parameter, and
$F^{ss^{\prime}}(\mathbf{k}^{\prime},\nu_{n})$ is the spin-resolved anomalous
propagator, describing the propagation of a hole that gets reflected into a
particle or vice-versa. The function
$\widetilde{V}^{\sigma\sigma^{\prime}ss^{\prime}}_{\mathbf{k}\mathbf{k}^{\prime}}$
describes an effective interaction, computed by means of the fRG (see Methods
Section and Supplementary Information SI for a deeper explanation). To
determine the symmetry (and not the size) of the SC gap, we linearize Eq. (7)
in $\Delta^{\sigma\sigma^{\prime}}_{\mathbf{k}}$. The equation becomes
therefore an eigenvalue problem, where the eigenvector corresponding to the
largest positive eigenvalue gives the information on the symmetry of the
leading superconducting state.

Phase diagram and magnetic fluctuations. In Fig. 1, we show the computed
superconducting phase diagram as a function of chemical potential $\mu\leq 0$
and second-neighbor hopping $t^{\prime}\leq 0$. We also show the different
topologies of the two Fermi surfaces, defined as the zeros of
$E^{-}_{\mathbf{k}}$ (blue) and $E^{+}_{\mathbf{k}}$ (red), as well as the
lines along which the van Hove singularities of the two quasiparticle bands
reach the Fermi level. We note that the density of states and position in $k$
space of the two van Hove singularities are markedly different: VHS1 (left
line) has a larger density of states than VHS2 (right line) and is also
located further away from the $(\pi,0)$ point, see Supplementary Information
SIV. A large portion of the phase diagram displays a leading superconducting
state living in the $B_{1}$ representation of the combined spin-lattice
symmetry group of Hamiltonian (1) (see Method Section and Supplementary
Information SII for more details). The resulting gap function can be expressed
as

$\displaystyle\Delta^{B_{1}}_{\mathbf{k}}=$ $\displaystyle\alpha_{s}(\cos
k_{x}-\cos k_{y})t^{0}+\alpha_{t}\left(\sin k_{y}t^{1}+\sin k_{x}t^{2}\right)$
(8) $\displaystyle+\,\text{higher harmonics},$

where $t^{\mu}=i\tau^{\mu}\tau^{2}$ ($\tau^{\mu}$, $\mu=0,..,3$ are the Pauli
matrices), and $\alpha_{s}$ and $\alpha_{t}$ are free parameters. Note that
the singlet component of the gap (the one proportional to $t^{0}$) is the
$d_{x^{2}-y^{2}}$-wave SC order parameter one would obtain in the Hubbard
model without introducing the spin-orbit coupling. At larger absolute values
of the chemical potential, we obtain a phase belonging to the $A_{1}$
representation of the discrete symmetry group of the Hamiltonian. The SC gap
in this state is a superposition of an extended $s$-wave in the singlet
component and a helical $p$-wave in the triplet:

$\begin{split}\Delta^{A_{1}}_{\mathbf{k}}=&\alpha_{s}(\cos k_{x}+\cos
k_{y})t^{0}+\alpha_{t}\left(-\sin k_{y}t^{1}+\sin k_{x}t^{2}\right)\\\
&+\,\text{higher harmonics}.\end{split}$ (9)

For larger $-\mu$ and small $-t^{\prime}$, we also find a $B_{2}$ phase,
characterized, however, by rather small values of the leading eigenvalue,
resulting in low transition temperatures. The gap function in this phase is
given by:

$\begin{split}\Delta^{B_{2}}_{\mathbf{k}}=&\alpha_{s}\sin k_{x}\sin
k_{y}t^{0}+\alpha_{t}\left(\sin k_{x}t^{1}-\sin k_{y}t^{2}\right)\\\
&+\,\text{higher harmonics}.\end{split}$ (10)

We remark that the transition lines in Fig. 1 only mark the crossing between
two eigenvalues, and a superconducting state resulting from the superposition
of two neighboring phases can be expected near the phase boundaries. At the
bottom left of the phase diagram in Fig. 1, that is, for small $-\mu$ and
$-t^{\prime}$, we labeled as $iAF$ (incommensurate antiferromagnetism) the
region where the fRG flow finds magnetism as the leading instability.

In the upper left panel of Fig. 2, we plot the parameter $\eta$,

$\eta=\frac{\sum_{\mu=1}^{3}\int_{\mathbf{k}}\mathrm{Tr}\left[\frac{1}{2}\Delta_{\mathbf{k}}(t^{\mu})^{\dagger}\right]^{2}}{\sum_{\mu=0}^{3}\int_{\mathbf{k}}\mathrm{Tr}\left[\frac{1}{2}\Delta_{\mathbf{k}}(t^{\mu})^{\dagger}\right]^{2}},$
(11)

quantifying the strength of spin-triplet in the leading superconducting state.
Values of $\eta$ that are close to zero imply a singlet-dominated SC state,
while $\eta\sim 1$ implies a triplet-dominated state. We note that the $B_{1}$
phase is mostly singlet-dominated, reminiscent of the $d$-wave SC state
emerging in the Hubbard model without spin-orbit coupling. Conversely, the
$A_{1}$ phase is triplet-dominated in a broad region in between the two van-
Hove singularities for large values of $-t^{\prime}$. This property can be
understood by analyzing the spin fluctuations in the model, which are
responsible for pairing [24, 25]. In the insets of panels (c) and (d) of Fig.
2, we plot the momentum dependence of the transverse and longitudinal static
spin susceptibility, probing the spin fluctuations in and out of the plane
defined by the lattice, respectively. In the main panels, we show the ratio
between the susceptibility evaluated at $\mathbf{q}=\mathbf{0}$ and its
maximum in the Brillouin zone, which gives a good indication of the strength
of ferromagnetic fluctuations. When this ratio gets close to one, the largest
peak of the relative susceptibility is to be found close to
$\mathbf{q}=\mathbf{0}$, viceversa, when it is close to zero the peak is far
away from $\mathbf{q}=\mathbf{0}$. This is the case of the region where a
magnetic instability is found, due to a divergence of the transverse
susceptibility for $\mathbf{q}$ near the antiferromagnetic wave-vector
$(\pi,\pi)$. The comparison between panels (a), (c) and (d) of Fig. 2 is very
instructive. Indeed, we note that spin fluctuations peaked around $(\pi,\pi)$
tend to favor a $B_{1}$ singlet-dominated pairing, while ferromagnetic
fluctuations enhance a triplet-dominated $A_{1}$ state.

In panel (b) of Fig. 2, we also report the value of the leading
superconducting eigenvalue $\lambda_{\mathrm{SC}}$, which provides information
on the relative size of the superconducting transition temperature. Note that
$\lambda_{\mathrm{SC}}>1$ implies that the critical temperature is higher than
the temperature at which the calculations have been performed, that is,
$T=10^{-6}t$. The superconducting eigenvalue is maximum along the van Hove
lines and close to the magnetically unstable region, where spin fluctuations
are stronger in absolute value, and it decays by moving away from it. Along
the leftmost van Hove line, it decays mildly, remaining large even down to
$t^{\prime}=-0.5t$, producing robust spin-triplet-dominated superconductivity
in this region. Within the fRG+MFT scheme, the transition temperature $T_{c}$
can be estimated either as the temperature at which the vertex function
diverges or as the lowest temperature at which the mean-field gap equation
returns a zero order parameter, corresponding to the point where
$\lambda_{\mathrm{SC}}=1$. By construction, both criteria yield the same
transition temperature. The behavior of the critical temperature as a function
of the chemical potential and second neighbor hopping resembles the one of the
superconducting eigenvalue (panel (b) of Fig. 2) and it reaches its maximum
values, of the order of $10^{-2}t$, along the van Hove line, next to the
antiferromagnetic instability. For large values of $|t^{\prime}|$, $T_{c}$ is
of the order of $10^{-3}t$-$10^{-4}t$ in the vicinity of the van Hove
singularity, and it decays away from it.

Figure 3:  B1 pairing state. (a) Momentum dependence of the B1 gap function
for intraband order parameters at $T=0$. $k_{x}>0$ and $k_{x}<0$ shows the
contributions from the two different helicity bands. Dashed lines indicate the
Fermi surface of the corresponding helicity band. (b) Phase crystal forming at
the (11) edge. The sine of the phase $\theta$ of the $d$-wave superconducting
order parameter is plotted in real space at $T=0$ with the bulk value of
$\sin\theta$ being subtracted. At the (11) edge there are clear modulations
visible in $\sin\theta$. (c) Eigenvalues $E$ around zero energy for different
temperature $T/T_{c}$. $n$ is the eigen index with twice the total number of
sites $2N$ being subtracted. All energies are in units of t. (d) Density of
states $N(E)$ around zero energy for different temperatures $T/T_{c}$. The
inset shows the density of states at $T/T_{c}=0.23$ in a wider range.

Having understood the origin of different superconducting phases in Fig. 1, we
then focus on each of them individually by investigating the characteristic
points marked by black crosses in Fig. 1. In particular, we characterize the
topological properties of the superconducting phases by looking at the
corresponding edge states. We use a hybrid method of fRG+MFT where we obtain
the effective interaction strengths in different superconducting channels from
fRG. Using these interactions we perform a renormalized mean-field theory in
real space, as detailed in the Method Section and Supplementary Information
SIII. Furthermore, we scale the interaction strengths by 10 in order to
investigate the edge properties.

B1 pairing state. In Fig. 3, we investigate the properties of the $B_{1}$
phase for the parameters marked by cross I ($\mu=-0.5$ and $t^{\prime}=-0.5$)
in Fig. 1. We first plot in Fig. 3(a) the absolute value of the intraband
superconducting gap functions $\Delta^{++}_{k}$ (shown for $k_{x}<0$) and
$\Delta^{--}_{k}$ (shown for $k_{x}>0$) with a periodic boundary condition.
The interband components are small and hence we do not show them here. Since
the leading instability is $B_{1}$ with mostly singlet components (see Fig.
2(a)), the order parameter has primarily $d$-wave symmetry and consequently
the gap functions show nodes along the diagonal $k_{x}=\pm k_{y}$. This nodal
structure suggests that a (11) edge is pair-breaking [52, 53] and will form
topological zero-energy states [54, 55, 56, 57] protected by time-reversal
symmetry and translational symmetry along (11). To investigate the edge
properties, we consider an open boundary geometry (see Fig. 3(b)) with one
(11) edge. At intermediate temperatures ($1>T/T_{c}>0.17$), we find a flat
band of zero-energy states as shown by red dots in Fig. 3(c). The large
degeneracy of these zero-energy states make them thermodynamically unstable
and prone to symmetry breaking. Thus at lower temperatures, $T/T_{c}<0.17$,
the (11) edge hosts oscillations in the phase $\theta$ of the $d$-wave
superconducting order parameter, see (b). Such oscillations are often referred
to as phase crystals breaking both time-reversal symmetry and translational
symmetry along the (11) edge [47, 48, 49]. In comparison to earlier literature
on phase crystals, one remarkable feature of the phase crystals obtained here
is that they are robust to the additional presence of a subdominant $A_{1}$
superconducting order parameter originating in the bulk due to the spin-orbit
coupling. The formation of phase crystals Doppler shifts the zero-energy
states to finite energies, as shown in Fig. 3(c). With lowering temperature,
the shift increases. We also calculate the spatially averaged density of
states
$N(E)=1/N\sum_{i,n,\sigma}|u_{i\sigma}^{n}|^{2}\delta(E-E_{n})+|v_{i\sigma}^{n}|^{2}\delta(E+E_{n})$,
where $N$ is the total number of lattice sites, and $u_{i\sigma}^{n}$ and
$v_{i\sigma}^{n}$ are the eigenfunctions with eigenvalues $E$. To numerically
evaluate $N(E)$, we use a Lorentzian with fixed width 0.01 to calculate the
delta-function. As shown in Fig. 3(d), $N(E)$ shows a large zero-bias peak for
$T/T_{c}>0.17$ which does not change with increasing temperature. Due to the
formation of the phase crystal at $T/T_{c}\approx 0.17$, the zero-bias peak
gets suppressed for lower temperatures. We do not find a full gap in $N(E)$
even at the lowest temperature, suggesting that there is no formation of
extended $s$-wave phase at the edge.

Figure 4:  $\mathbf{A}_{1}$ pairing state. (a) Momentum dependence of the
$A_{1}$ gap function for intraband order paramters. $k_{x}>0$ and $k_{x}<0$
shows the contributions from the two different helicity bands. Dashed lines
indicate the Fermi surface of the corresponding helicity band. (b) Position
dependence of the wave-function amplitude of the lowest positive eigenstate
near zero-energy for a nanoribbon with open boundary conditions along $x$ and
periodic boundary conditions along $y$. The inset shows the eigenvalues. (c)
Momentum resolved density of states averaged over $x$, $\langle
A(k_{y},E,i_{x})\rangle_{i_{x}}$, showing linearly dispersing edge states. (d)
Momentum resolved spin density of states averaged over a width near the left
(10) edge of the nanoribbon, $\langle A_{s}(k_{y},E,i_{x})\rangle_{i_{x}}$,
showing the helical nature of the edge states. All the plots are for $T=0$.
All energies are in units of t. The maximum values of the colorbars in (c) and
(d) are cut for better visibility.

A1 pairing state. We now investigate the topological edge states of the
$A_{1}$ phase for the parameters marked by cross III ($\mu=-1.6$ and
$t^{\prime}=-0.5$) in Fig. 1. In Fig. 4(a), we first show the absolute values
of intraband superconducting gap functions obtained with periodic boundary
condition. As seen, the gap functions feature point nodes, but away from the
Fermi surface. These nodes are accidental as they can be moved by changing
parameters of the model. Due to the lack of nodes on the Fermi surface, we
find a full gap in the spectrum (not shown). With this $A_{1}$ pairing, the
system is a time-reversal symmetric superconductor and belongs to the DIII
class, characterized by a $\mathbb{Z}_{2}$ invariant in two dimension [3]. In
our case, there is one pocket around the M point with negative pairing and the
system is topologically nontrivial according to
$N_{2D}=\Pi_{s}[\text{sgn}(\Delta_{s})]^{m_{s}}$, where $m_{s}$ is the number
of time-reversal invariant points enclosed by the $s$th Fermi surface [58]. To
study the topological edge states, we then consider a nanoribbon geometry with
open boundary conditions along $x$ and periodic boundary conditions along $y$.
Here, we show the results for the nanorribbon geometry with no self-
consistency on the order parameters, in order to reach a large system size. We
obtain the self-consistent order parameters with periodic boundary conditions
along both $x$ and $y$, and using this self-consistent solution we employ the
nanoribbon geometry with no further self-consistency. We have however verified
for a smaller system size that all the edge features are insensitive to self-
consistency. In the inset of Fig. 4(b), we show the eigenvalues for the
nanoribbon geometry. As seen, there are in-gap edge states formed for
$|E|<0.2$. We also show the wavefunction
$|\psi_{n}(i)|^{2}=\sum_{\sigma}(|u_{i\sigma}^{n}|^{2}+|v_{i\sigma}^{n}|^{2})$
for the lowest positive eigenstate with $E\approx 0$ in the main panel of Fig.
4(b). The wavefunction is localized to the open (10) edges with oscillating
weights going away from the edge. The localization length is less than the
half of the length along $x$-direction, thus prohibiting any interference of
the two edges. Using translational symmetry along the $y$-direction, we look
at the edge spectrum by calculating the momentum resolved density of states
$A(k_{y},E,i_{x})=\sum_{n,\sigma}|u_{i_{x},k_{y},\sigma}^{n}|^{2}\delta(E-E_{n})$
and spin density of states
$A_{s}(k_{y},E,i_{x})=\sum_{n,\sigma}\mathrm{sgn}(\sigma)|u_{i_{x},k_{y},\sigma}^{n}|^{2}\delta(E-E_{n})$,
where $u_{i_{x},k_{y},\sigma}^{n}$ is the Fourier transform of
$u_{i,\sigma}^{n}$ along the $y$-direction with $i_{x}$ being the
$x$-coordinate of the lattice site $i$ and $E_{n}$ is the eigenenergy of the
state $n$. In Fig. 4(c), we show $A(k_{y},E,i_{x})$ averaged over all $i_{x}$.
We find a linearly dispersing Dirac edge spectrum inside the gap $E=0.2$ at
$k_{y}=\pi$. To identify the spin character of the edge states, we show in
Fig. 4(d) $A_{s}(k_{y},E,i_{x})$ averaged over $i_{x}=0,10$ from the left (10)
edge. The spectrum clearly shows a helical nature of the edge states with
positive values for positive velocities and negative for negative velocities.
The corresponding $A_{s}(k_{y},E,i_{x})$ averaged over $i_{x}=0,10$ from the
right (10) edge shows an opposite behavior, not shown. The helical Majorana
edge modes clearly demonstrate the nontrivial topology of the $A_{1}$
superconducting phase.

Figure 5:  $\mathbf{A}_{1}\mathbf{+}\mathbf{i}\mathbf{B}_{1}$ pairing state.
(a) Momentum dependence of the $A_{1}+iB_{1}$ gap function for intraband order
parameters. $k_{x}>0$ and $k_{x}<0$ shows the contributions from the two
different helicity bands. Dashed lines indicate the Fermi surface of the
corresponding helicity band. (b) Condenstation energy as a function of the
relative phase $\phi$ between $A_{1}$ and $B_{1}$ order parameters, showing
$\pi/2$ as the lowest in energy. (c) Eigenvalues with open boundary conditions
along both $x$ and $y$, showing four Majorana states. Inset shows a zoomed
view of the four Majorana states. (d) Spin projected wave function of the four
Majarana states as a function of position. All the plots are for $T=0$ and the
energies are in units of t.

A${}_{1}+\mathbf{i}$B1 pairing state. The transition from the $B_{1}$
superconducting state to the $A_{1}$ state with increasing $|\mu|$ in Fig. 1
gives the possibility of a coexisting phase where both order parameters are
comparable. To explore this possibility, we fix the parameters to the values
marked by cross II ($\mu=-1.2$ and $t^{\prime}=-0.5$) in Fig. 1, which is near
the phase boundary. Interestingly, we find a time-reversal symmetry breaking
$A_{1}+iB_{1}$ superconducting phase as the lowest energy state. The
corresponding gap functions are shown in Fig. 5(a). Again, the nodes do not
lie on the Fermi surface making the state fully gapped with an energy gap of
$E=0.25$. To verify that the relative phase of $\pi/2$ between the $A_{1}$ and
the $B_{1}$ order parameters is indeed the global minimum, we compare the
condensation energies of the superconducting states with different relative
phase $\phi$. As shown in Fig. 5(b), $\phi=\pm\pi/2$ gives the largest
condensation energy and consequently $A_{1}\pm iB_{1}$ is lowest in energy.
Due to the imaginary $B_{1}$ pairing component, time-reversal symmetry and
four-fold rotational symmetry $C_{4}$ are broken in this pairing state but
their combination is preserved, leading to a gap opening in the helical
Majorana edge states. However, owing to the sign change of the $B_{1}$ pairing
under $C_{4}$, the mass terms for adjacent edges [(10) and (01) edges] have
opposite sign. Therefore, when these two edges meet at the corners, the mass
term vanishes and Majorana corner modes are generated, realizing a second-
order topological superconductor [38] (a detailed analysis is given in the
Supplementary Information SV). To demonstrate this nontrivial topology, we
study the corner states in a geometry with open boundary conditions along both
$x$\- and $y$-directions. Similar to the $A_{1}$ phase, we only consider non
self-consistent eigenstates with open boundaries taking the fully self-
consistent bulk solutions. Remarkably, we find four Majorana zero-energy
states located in the superconducting gap, as shown in Fig. 5(c). The
wavefunctions corresponding to these Majorana states are localized at the four
corners, not shown. Spin characteristics of these Majorana states can be
visualized by looking at the spin projected wavefunctions
$\varrho_{s}(i)=\sum_{\sigma
n^{\prime}}\mathrm{sgn}(\sigma)\left(|u_{i\sigma}^{n^{\prime}}|^{2}+|v_{i\sigma}^{n^{\prime}}|^{2}\right)$
where $n^{\prime}$ are the four zero-energy states. In Fig. 5(d), we see that
$\varrho_{s}(i)$ is localized at the four corners with alternating signs for
alternating corners. We have also verified the presence of Majorana corner
states in a self-consistent calculation with edges for smaller system sizes.
Therefore, the Majorana corner states in our numerical calculations confirm
that the $A_{1}+iB_{1}$ state is a higher-order topological superconductor.

Discussion
We have shown that the combined effects of van Hove singularities, Rashba
spin-orbit coupling (SOC), and repulsive Hubbard interactions give rise to a
cascade of topological superconducting states, including a nodal $B_{1}$ state
($d$-wave like) with flat-band Majorana edge modes, a fully gapped $A_{1}$
state ($p$-wave like) with helical Majorana edge modes, and a time-reversal
breaking $A_{1}+iB_{1}$ state with Majorana corner modes. While the first two
pairing states have first order topology, characterized by a winding number
and a Pfaffian ${Z}_{2}$ invariant, respectively, the third one exhibits
higher-order topology. Since Rashba SOC breaks inversion symmetry, all three
phases consist of a mixture of spin-singlet and spin-triplet Cooper pairs,
whose relative ratio, however, is vastly different: The $B_{1}$ and $A_{1}$
states have dominant singlet and triplet pairing, respectively, whereas the
$A_{1}+iB_{1}$ is composed of an almost equal amount of singlet and triplet
Cooper pairs.The pairing glue for these superconductors is provided by ferro-
and antiferromagnetic spin fluctuations, which are strongly enhanced near the
van Hove fillings. Interestingly, the Rashba SOC splits the van Hove
singularity at the $X$ point of the square lattice into two _different types_
, namely, VHS1 which shifts along the $X-\Gamma$ direction and VHS2 which is
moved in the $X-M$ direction, see Supplementary Information SIV. For fillings
near VHS2 the spin fluctuations are largely ferromagnetic leading to $A_{1}$
superconductivity with predominantly triplet pairing, whereas for fillings
close to VHS1 incommensurate spin fluctuations give rise to $A_{1}+iB_{1}$
superconductivity. Remarkably, the divergence at VHS1 is considerably stronger
than the one at VHS2, implying that the higher-order topological
superconductor has the largest $T_{c}$. As a consequence, correlated Rashba
materials whose Fermi level can be tuned to VHS1 are good candidates for
intrinsic higher-order topological superconductors.

It is remarkable that we find the above three topological superconductors in a
single phase diagram with only onsite Hubbard interactions. We emphasize that
they emerge robustly in the Rashba-Hubbard model, independent of the size of
the Rashba SOC and survive the inclusion of further neighbor hoppings.
Although we have presented here only results for the square lattice, we expect
these topological superconductors to arise also in other 2D lattices, e.g.,
the triangular or the honeycomb lattice, although with modified irreps, due to
the different spin-lattice symmetry groups. For example, in the triangular
lattice the $B_{1}$ phase would become an $E_{2}+iE^{\prime}_{2}$ pairing
state (chiral $d$-wave like). Finally, we remark that due to their topological
nature these superconducting phases survive moderately strong disorder [59,
60, 61, 62, 63, 64], which we exemplified in the Supplementary Information SVI
for nonmagnetic onsite disorder with box distributions.

Our results provide a guide to understand and design topological
superconductivity in heavy-atom superlattices and van der Waals materials. The
high variability of these materials may allow to tune the Fermi level to the
van Hove fillings, such that an intrinsic topological superconductor with
large $T_{c}$ can be realized. For example, in LAO/STO it is possible to tune
the carrier density, and therefore the Fermi surface, by electric gating [20].
This may open up a route to use intrinsic topological superconductors for the
design of quantum information devices. Moreover, the $A_{1}+iB_{1}$
superconductor may show interesting diode [65, 66] and piezoelectric effects
[67], since it breaks both inversion and time-reversal symmetry.

Methods

fRG+MFT approach. To determine the leading instabilities of the Rashba-Hubbard
model, we employ a truncated unity renormalization procedure using the
temperature flow. For simplicity we truncate the effective action at the two-
particle level and ingore the flow of the self-energy. The particle-particle
and particle-hole channels are expanded in a complete basis of form factors
(truncated unity approach). To continue the flow below the stopping
temperature and to detect the superconducting instabilities, the functional
renormalization group is combined with a momentum-space mean field theory. In
the presence of disorder or edges the renormalized mean-field equations are
solved in real space. For more details see Supplementary Information SI and
SIII.

Symmetry classification of superconducting gap. The spin-lattice symmetries of
the square lattice Rashba-Hubbard model form a group with eight elements that
is equivalent to the point group $C_{4v}$. Its irreducible represnetations
consist of four one-dimensional irreps ($A_{1}$, $A_{2}$, $B_{1}$, $B_{2}$)
and one two-dimensional irrep ($E$). The symmetries of the superconducting
gaps can therefore be classified according to these irreps, see SII for more
details.

Acknowledgements. We are thankful to A. Greco and M. Hirschmann for
enlightening discussions. D.C. and A.P.S. thank the Max-Planck-Institute FKF
Stuttgart and the YITP Kyoto, respectively, for hospitality during the course
of this project.

Data availability. The numerical data of the plots within this paper are
available from the authors upon reasonable request.

## References

  * [1] Beenakker, C.  Search for majorana fermions in superconductors.  _Annual Review of Condensed Matter Physics_ 4, 113–136 (2013).  URL https://doi.org/10.1146/annurev-conmatphys-030212-184337.
  * [2] Schnyder, A. P. & Brydon, P. M. R.  Topological surface states in nodal superconductors.  _Journal of Physics: Condensed Matter_ 27, 243201 (2015).  URL https://dx.doi.org/10.1088/0953-8984/27/24/243201.
  * [3] Chiu, C.-K., Teo, J. C. Y., Schnyder, A. P. & Ryu, S.  Classification of topological quantum matter with symmetries.  _Rev. Mod. Phys._ 88, 035005 (2016).  URL https://link.aps.org/doi/10.1103/RevModPhys.88.035005.
  * [4] Sato, M. & Ando, Y.  Topological superconductors: a review.  _Reports on Progress in Physics_ 80, 076501 (2017).  URL https://dx.doi.org/10.1088/1361-6633/aa6ac7.
  * [5] Sharma, M. M., Sharma, P., Karn, N. K. & Awana, V. P. S.  Comprehensive review on topological superconducting materials and interfaces.  _Superconductor Science and Technology_ 35, 083003 (2022).  URL https://dx.doi.org/10.1088/1361-6668/ac6987.
  * [6] von Rohr, F. O.  Chemical principles of intrinsic topological superconductors (2023).  eprint 2304.00090.
  * [7] Flensberg, K., von Oppen, F. & Stern, A.  Engineered platforms for topological superconductivity and majorana zero modes.  _Nature Reviews Materials_ 6, 944–958 (2021).  URL https://doi.org/10.1038/s41578-021-00336-6.
  * [8] Hao, N. & Hu, J.  Topological quantum states of matter in iron-based superconductors: from concept to material realization.  _National Science Review_ 6, 213–226 (2018).  URL https://doi.org/10.1093/nsr/nwy142.
  * [9] Wu, X., Zhang, R.-X., Xu, G., Hu, J. & Liu, C.-X.  In the pursuit of majorana modes in iron-based high-tc superconductors.  _Memorial Volume for Shoucheng Zhang_ 35–60 (2021).  URL https://www.worldscientific.com/doi/abs/10.1142/9789811231711_0005.
  * [10] Yin, J.-X. _et al._ Observation of a robust zero-energy bound state in iron-based superconductor fe(te,se).  _Nature Physics_ 11, 543–546 (2015).  URL https://doi.org/10.1038/nphys3371.
  * [11] Wang, D. _et al._ Evidence for majorana bound states in an iron-based superconductor.  _Science_ 362, 333–335 (2018).  URL https://www.science.org/doi/abs/10.1126/science.aao1797.
  * [12] Shimozawa, M., Goh, S. K., Shibauchi, T. & Matsuda, Y.  From kondo lattices to kondo superlattices.  _Reports on Progress in Physics_ 79, 074503 (2016).  URL https://dx.doi.org/10.1088/0034-4885/79/7/074503.
  * [13] Novoselov, K. S., Mishchenko, A., Carvalho, A. & Neto, A. H. C.  2d materials and van der waals heterostructures.  _Science_ 353, aac9439 (2016).  URL https://www.science.org/doi/abs/10.1126/science.aac9439.
  * [14] Liu, Y. _et al._ Van der waals heterostructures and devices.  _Nature Reviews Materials_ 1, 16042 (2016).  URL https://doi.org/10.1038/natrevmats.2016.42.
  * [15] Bian, R. _et al._ Recent progress in the synthesis of novel two-dimensional van der Waals materials.  _National Science Review_ 9 (2021).  URL https://doi.org/10.1093/nsr/nwab164.  Nwab164.
  * [16] Jang, B. G., Lee, C., Zhu, J.-X. & Shim, J. H.  Exploring two-dimensional van der waals heavy-fermion material: Data mining theoretical approach.  _npj 2D Materials and Applications_ 6, 80 (2022).  URL https://doi.org/10.1038/s41699-022-00357-x.
  * [17] Hegger, H. _et al._ Pressure-induced superconductivity in quasi-2d ${\mathrm{cerhin}}_{5}$.  _Phys. Rev. Lett._ 84, 4986–4989 (2000).  URL https://link.aps.org/doi/10.1103/PhysRevLett.84.4986.
  * [18] Posey, V. & et al.  Existence of heavy fermion behavior in a van der waals layered material (2023).  URL https://meetings.aps.org/Meeting/MAR23/Session/G29.5.  APS March Meeting 2023.
  * [19] Reyren, N. _et al._ Superconducting interfaces between insulating oxides.  _Science_ 317, 1196–1199 (2007).  URL https://www.science.org/doi/abs/10.1126/science.1146006.
  * [20] Caviglia, A. D. _et al._ Electric field control of the laalo3/srtio3 interface ground state.  _Nature_ 456, 624–627 (2008).  URL https://doi.org/10.1038/nature07576.
  * [21] Liu, C. _et al._ Two-dimensional superconductivity and anisotropic transport at KTa$\mathrm{O}_{3}$ (111) interfaces.  _Science_ 371, 716–721 (2021).  URL https://www.science.org/doi/abs/10.1126/science.aba5511.
  * [22] Mizukami, Y. _et al._ Extremely strong-coupling superconductivity in artificial two-dimensional kondo lattices.  _Nature Physics_ 7, 849–853 (2011).  URL https://doi.org/10.1038/nphys2112.
  * [23] Shimozawa, M. _et al._ Controllable rashba spin-orbit interaction in artificially engineered superlattices involving the heavy-fermion superconductor ${\mathrm{cecoin}}_{5}$.  _Phys. Rev. Lett._ 112, 156404 (2014).  URL https://link.aps.org/doi/10.1103/PhysRevLett.112.156404.
  * [24] Greco, A. & Schnyder, A. P.  Mechanism for Unconventional Superconductivity in the Hole-Doped Rashba-Hubbard Model.  _Phys. Rev. Lett._ 120, 177002 (2018).  URL https://link.aps.org/doi/10.1103/PhysRevLett.120.177002.
  * [25] Greco, A., Bejas, M. & Schnyder, A. P.  Ferromagnetic fluctuations in the rashba-hubbard model.  _Phys. Rev. B_ 101, 174420 (2020).  URL https://link.aps.org/doi/10.1103/PhysRevB.101.174420.
  * [26] Ghadimi, R., Kargarian, M. & Jafari, S. A.  Competing superconducting phases in the interacting two-dimensional electron gas with strong rashba spin-orbit coupling.  _Phys. Rev. B_ 99, 115122 (2019).  URL https://link.aps.org/doi/10.1103/PhysRevB.99.115122.
  * [27] Wolf, S. & Rachel, S.  Spin-orbit coupled superconductivity: Rashba-Hubbard model on the square lattice.  _Phys. Rev. B_ 102, 174512 (2020).  URL https://link.aps.org/doi/10.1103/PhysRevB.102.174512.
  * [28] Nogaki, K. & Yanase, Y.  Strongly parity-mixed superconductivity in the rashba-hubbard model.  _Phys. Rev. B_ 102, 165114 (2020).  URL https://link.aps.org/doi/10.1103/PhysRevB.102.165114.
  * [29] Biderang, M., Zare, M.-H. & Sirker, J.  Proximity-driven ferromagnetism and superconductivity in the triangular rashba-hubbard model.  _Phys. Rev. B_ 105, 064504 (2022).  URL https://link.aps.org/doi/10.1103/PhysRevB.105.064504.
  * [30] Kawano, M. & Hotta, C.  Phase diagram of the square-lattice hubbard model with rashba-type antisymmetric spin-orbit coupling.  _Phys. Rev. B_ 107, 045123 (2023).  URL https://link.aps.org/doi/10.1103/PhysRevB.107.045123.
  * [31] Beyer, J. _et al._ Rashba spin-orbit coupling in the square-lattice hubbard model: A truncated-unity functional renormalization group study.  _Phys. Rev. B_ 107, 125115 (2023).  URL https://link.aps.org/doi/10.1103/PhysRevB.107.125115.
  * [32] Kopietz, P., Bartosch, L. & Schütz, F.  _Introduction to the Functional Renormalization Group_ (Springer Berlin Heidelberg, 2010).  URL https://doi.org/10.1007/978-3-642-05094-7.
  * [33] Metzner, W., Salmhofer, M., Honerkamp, C., Meden, V. & Schönhammer, K.  Functional renormalization group approach to correlated fermion systems.  _Rev. Mod. Phys._ 84, 299–352 (2012).  URL https://link.aps.org/doi/10.1103/RevModPhys.84.299.
  * [34] Wang, J., Eberlein, A. & Metzner, W.  Competing order in correlated electron systems made simple: Consistent fusion of functional renormalization and mean-field theory.  _Phys. Rev. B_ 89, 121116(R) (2014).  URL https://link.aps.org/doi/10.1103/PhysRevB.89.121116.
  * [35] Yamase, H., Eberlein, A. & Metzner, W.  Coexistence of Incommensurate Magnetism and Superconductivity in the Two-Dimensional Hubbard Model.  _Phys. Rev. Lett._ 116, 096402 (2016).  URL https://link.aps.org/doi/10.1103/PhysRevLett.116.096402.
  * [36] Vilardi, D., Bonetti, P. M. & Metzner, W.  Dynamical functional renormalization group computation of order parameters and critical temperatures in the two-dimensional Hubbard model.  _Phys. Rev. B_ 102, 245128 (2020).  URL https://link.aps.org/doi/10.1103/PhysRevB.102.245128.
  * [37] Bonetti, P. M.  Accessing the ordered phase of correlated Fermi systems: Vertex bosonization and mean-field theory within the functional renormalization group.  _Phys. Rev. B_ 102, 235160 (2020).  URL https://link.aps.org/doi/10.1103/PhysRevB.102.235160.
  * [38] Wang, Y., Lin, M. & Hughes, T. L.  Weak-pairing higher order topological superconductors.  _Phys. Rev. B_ 98, 165144 (2018).  URL https://link.aps.org/doi/10.1103/PhysRevB.98.165144.
  * [39] Yan, Z., Song, F. & Wang, Z.  Majorana corner modes in a high-temperature platform.  _Phys. Rev. Lett._ 121, 096803 (2018).  URL https://link.aps.org/doi/10.1103/PhysRevLett.121.096803.
  * [40] Wang, Q., Liu, C.-C., Lu, Y.-M. & Zhang, F.  High-temperature majorana corner states.  _Phys. Rev. Lett._ 121, 186801 (2018).  URL https://link.aps.org/doi/10.1103/PhysRevLett.121.186801.
  * [41] Wu, X., Liu, X., Thomale, R. & Liu, C.-X.  High-Tc superconductor Fe(Se,Te) monolayer: an intrinsic, scalable and electrically tunable Majorana platform.  _National Science Review_ 9 (2021).  URL https://doi.org/10.1093/nsr/nwab087.  Nwab087.
  * [42] Hsu, Y.-T., Cole, W. S., Zhang, R.-X. & Sau, J. D.  Inversion-protected higher-order topological superconductivity in monolayer ${\mathrm{wte}}_{2}$.  _Phys. Rev. Lett._ 125, 097001 (2020).  URL https://link.aps.org/doi/10.1103/PhysRevLett.125.097001.
  * [43] Kheirkhah, M., Yan, Z., Nagai, Y. & Marsiglio, F.  First- and second-order topological superconductivity and temperature-driven topological phase transitions in the extended hubbard model with spin-orbit coupling.  _Phys. Rev. Lett._ 125, 017001 (2020).  URL https://link.aps.org/doi/10.1103/PhysRevLett.125.017001.
  * [44] Wu, X. _et al._ Boundary-obstructed topological high-${\mathit{t}}_{c}$ superconductivity in iron pnictides.  _Phys. Rev. X_ 10, 041014 (2020).  URL https://link.aps.org/doi/10.1103/PhysRevX.10.041014.
  * [45] Ikegaya, S., Rui, W. B., Manske, D. & Schnyder, A. P.  Tunable majorana corner modes in noncentrosymmetric superconductors: Tunneling spectroscopy and edge imperfections.  _Phys. Rev. Res._ 3, 023007 (2021).  URL https://link.aps.org/doi/10.1103/PhysRevResearch.3.023007.
  * [46] Li, T., Geier, M., Ingham, J. & Scammell, H. D.  Higher-order topological superconductivity from repulsive interactions in kagome and honeycomb systems.  _2D Materials_ 9, 015031 (2021).  URL https://dx.doi.org/10.1088/2053-1583/ac4060.
  * [47] Håkansson, M., Löfwander, T. & Fogelström, M.  Spontaneously broken time-reversal symmetry in high-temperature superconductors.  _Nat. Phys._ 11, 755–760 (2015).  URL https://doi.org/10.1038/nphys3383.
  * [48] Holmvall, P., Fogelström, M., Löfwander, T. & Vorontsov, A. B.  Phase crystals.  _Phys. Rev. Research_ 2, 013104 (2020).  URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.013104.
  * [49] Chakraborty, D., Löfwander, T., Fogelström, M. & Black-Schaffer, A. M.  Disorder-robust phase crystal in high-temperature superconductors stabilized by strong correlations.  _npj Quantum Materials_ 7, 44 (2022).  URL https://doi.org/10.1038/s41535-022-00450-w.
  * [50] Nayak, C., Simon, S. H., Stern, A., Freedman, M. & Das Sarma, S.  Non-abelian anyons and topological quantum computation.  _Rev. Mod. Phys._ 80, 1083–1159 (2008).  URL https://link.aps.org/doi/10.1103/RevModPhys.80.1083.
  * [51] Zhang, S.-B. _et al._ Topological and holonomic quantum computation based on second-order topological superconductors.  _Phys. Rev. Res._ 2, 043025 (2020).  URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.043025.
  * [52] Kashiwaya, S. & Tanaka, Y.  Tunnelling effects on surface bound states in unconventional superconductors.  _Rep. Prog. Phys._ 63, 1641–1724 (2000).  URL https://doi.org/10.1088%2F0034-4885%2F63%2F10%2F202.
  * [53] Löfwander, T., Shumeiko, V. S. & Wendin, G.  Andreev bound states in high-tcsuperconducting junctions.  _Supercond. Sci. Technol._ 14, R53–R77 (2001).  URL https://doi.org/10.1088%2F0953-2048%2F14%2F5%2F201.
  * [54] Ryu, S. & Hatsugai, Y.  Topological origin of zero-energy edge states in particle-hole symmetric systems.  _Phys. Rev. Lett._ 89, 077002 (2002).  URL https://link.aps.org/doi/10.1103/PhysRevLett.89.077002.
  * [55] Sato, M., Tanaka, Y., Yada, K. & Yokoyama, T.  Topology of andreev bound states with flat dispersion.  _Phys. Rev. B_ 83, 224511 (2011).  URL https://link.aps.org/doi/10.1103/PhysRevB.83.224511.
  * [56] Schnyder, A. P. & Brydon, P. M. R.  Topological surface states in nodal superconductors.  _J. Condens. Matter Phys._ 27, 243201 (2015).  URL https://doi.org/10.1088%2F0953-8984%2F27%2F24%2F243201.
  * [57] Schnyder, A. P., Timm, C. & Brydon, P. M. R.  Edge currents as a signature of flatbands in topological superconductors.  _Phys. Rev. Lett._ 111, 077001 (2013).  URL https://link.aps.org/doi/10.1103/PhysRevLett.111.077001.
  * [58] Qi, X.-L., Hughes, T. L. & Zhang, S.-C.  Topological invariants for the fermi surface of a time-reversal-invariant superconductor.  _Phys. Rev. B_ 81, 134508 (2010).  URL https://link.aps.org/doi/10.1103/PhysRevB.81.134508.
  * [59] Queiroz, R. & Schnyder, A. P.  Stability of flat-band edge states in topological superconductors without inversion center.  _Phys. Rev. B_ 89, 054501 (2014).  URL https://link.aps.org/doi/10.1103/PhysRevB.89.054501.
  * [60] Queiroz, R. & Schnyder, A. P.  Helical majorana surface states of strongly disordered topological superconductors with time-reversal symmetry.  _Phys. Rev. B_ 91, 014202 (2015).  URL https://link.aps.org/doi/10.1103/PhysRevB.91.014202.
  * [61] Lu, Y., Virtanen, P. & Heikkilä, T. T.  Effect of disorder on majorana localization in topological superconductors: A quasiclassical approach.  _Phys. Rev. B_ 102, 224510 (2020).  URL https://link.aps.org/doi/10.1103/PhysRevB.102.224510.
  * [62] Garg, A., Randeria, M. & Trivedi, N.  Strong correlations make high-temperature superconductors robust against disorder.  _Nat. Phys._ 4, 762–765 (2008).  URL https://doi.org/10.1038/nphys1026.
  * [63] Chakraborty, D. & Ghosal, A.  Fate of disorder-induced inhomogeneities in strongly correlated d-wave superconductors.  _New J. Phys._ 16, 103018 (2014).  URL https://doi.org/10.1088%2F1367-2630%2F16%2F10%2F103018.
  * [64] Ghosal, A., Chakraborty, D. & Kaushal, N.  Prospects of anderson’s theorem for disordered cuprate superconductors.  _Physica B Condens. Matter_ 536, 867 – 876 (2018).  URL http://www.sciencedirect.com/science/article/pii/S0921452617305252.
  * [65] Ando, F. _et al._ Observation of superconducting diode effect.  _Nature_ 584, 373–376 (2020).  URL https://doi.org/10.1038/s41586-020-2590-4.
  * [66] Daido, A., Ikeda, Y. & Yanase, Y.  Intrinsic superconducting diode effect.  _Phys. Rev. Lett._ 128, 037001 (2022).  URL https://link.aps.org/doi/10.1103/PhysRevLett.128.037001.
  * [67] Chazono, M., Watanabe, H. & Yanase, Y.  Superconducting piezoelectric effect.  _Phys. Rev. B_ 105, 024509 (2022).  URL https://link.aps.org/doi/10.1103/PhysRevB.105.024509.
  * [68] Honerkamp, C. & Salmhofer, M.  Temperature-flow renormalization group and the competition between superconductivity and ferromagnetism.  _Phys. Rev. B_ 64, 184516 (2001).  URL https://link.aps.org/doi/10.1103/PhysRevB.64.184516.
  * [69] Wetterich, C.  Exact evolution equation for the effective potential.  _Phys. Lett. B_ 301, 90 – 94 (1993).  URL http://www.sciencedirect.com/science/article/pii/037026939390726X.
  * [70] Husemann, C. & Salmhofer, M.  Efficient parametrization of the vertex function, $\Omega$ scheme, and the $t,{t}^{{}^{\prime}}$ Hubbard model at van Hove filling.  _Phys. Rev. B_ 79, 195125 (2009).  URL https://link.aps.org/doi/10.1103/PhysRevB.79.195125.
  * [71] Husemann, C., Giering, K.-U. & Salmhofer, M.  Frequency-dependent vertex functions of the ($t,{t}^{{}^{\prime}}$) Hubbard model at weak coupling.  _Phys. Rev. B_ 85, 075121 (2012).  URL https://link.aps.org/doi/10.1103/PhysRevB.85.075121.
  * [72] Lichtenstein, J. _et al._ High-performance functional Renormalization Group calculations for interacting fermions.  _Comput. Phys. Commun._ 213, 100 – 110 (2017).  URL http://www.sciencedirect.com/science/article/pii/S0010465516303927.
  * [73] Schober, G. A. H., Ehrlich, J., Reckling, T. & Honerkamp, C.  Truncated-Unity Functional Renormalization Group for Multiband Systems With Spin-Orbit Coupling.  _Front. Phys._ 6 (2018).  URL https://www.frontiersin.org/articles/10.3389/fphy.2018.00032.
  * [74] Platt, C., Hanke, W. & Thomale, R.  Functional renormalization group for multi-orbital fermi surface instabilities.  _Adv. Phys._ 62, 453–562 (2013).  URL https://doi.org/10.1080/00018732.2013.862020.
  * [75] Zhu, J.  _Bogoliubov-de Gennes Method and Its Applications_.  Lecture Notes in Physics (Springer International Publishing, 2016).  URL https://books.google.se/books?id=Tep6DAAAQBAJ.
  * [76] Schnyder, A. P., Brydon, P. M. R. & Timm, C.  Types of topological surface states in nodal noncentrosymmetric superconductors.  _Phys. Rev. B_ 85, 024522 (2012).  URL https://link.aps.org/doi/10.1103/PhysRevB.85.024522.

Supplementary Information - van Hove, Rashba, and Hubbard meet to form first-
order and higher-order topological superconductors

This Supplementary Information contains six sections. In SI we discuss the
details of the truncated unity renormalization group procedure and explain how
it is combined with a momentum-space mean-field theory. In SII we present the
symmetry classification of the superconducting states in terms of the irreps
of the spin-lattice group. In SIII we provide the details of the real space
renormalized mean-field theory and show how it is matched to the fRG. In SIV
effective models for the van Hove singularities are derived and analyzed. In
SV the corner modes of the $A_{1}+iB_{1}$ pairing state are discussed. The
effects of disorder are studied in SVI.

##  SI fRG+MFT method

In this section we provide some details on the fRG+MFT method, employed to
determine the leading superconducting symmetry in the Rashba-Hubbard model.

###  SI.1 Functional renormalization group

In order to determine the system’s instabilities and the leading
superconducting symmetry, we run a truncated functional renormalization group
flow [32, 33]. The flow is implemented by a progressive integration of
fermionic modes by introducing a flowing cutoff in the bare propagator. Here,
we choose the so-called _temperature flow_ [68]. By rescaling the Grassmann
fields according to

$\displaystyle\eta_{\mathbf{k},\sigma}=T^{\frac{3}{4}}\psi_{\mathbf{k},\sigma},$
(S1a)
$\displaystyle\overline{\eta}_{\mathbf{k},\sigma}=T^{\frac{3}{4}}\overline{\psi}_{\mathbf{k},\sigma},$
(S1b)

where $T$ is the temperature, we obtain the action

$\mathcal{S}^{T}[\eta,\overline{\eta}]=\sum_{n}\int_{\mathbf{k}}\overline{\eta}_{\mathbf{k},n}\left[T^{-\frac{1}{2}}\left(i\nu_{n}-\hat{\mathcal{H}}_{\mathbf{k}}\right)\right]\eta_{\mathbf{k},n}+U\sum_{n_{1},n_{2},n_{3}}\int_{\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3}}\overline{\eta}_{\mathbf{k}_{1},n_{1},\uparrow}\overline{\eta}_{\mathbf{k}_{2},n_{2},\downarrow}\eta_{\mathbf{k}_{4},n_{4},\downarrow}\eta_{\mathbf{k}_{3},n_{3},\uparrow},$

where $n$ labels fermionic Matsubara frequencies $\nu_{n}=(2n+1)\pi T$
($n\in\mathbb{Z}$),
$\mathbf{k}_{4}=\mathbf{k}_{1}+\mathbf{k}_{2}-\mathbf{k}_{3}$ and
$n_{4}=n_{1}+n_{2}-n_{3}$, and $\hat{\mathcal{H}}_{\mathbf{k}}$ is defined in
Eq. (3). We note that upon rescaling (S1), only the quadratic part of the
action remains temperature-dependent, with bare propagator

$\hat{G}_{0}^{T}(\mathbf{k},\nu_{n})=T^{\frac{1}{2}}\left[i\nu_{n}-\hat{\mathcal{H}}_{\mathbf{k}}\right]^{-1}.$
(S2)

The temperature therefore takes the role of a renormalization group scale and
one can derive an exact evolution equation for the effective action functional
[69, 33]

$\frac{d}{dT}\Gamma^{T}[\eta,\overline{\eta}]=-\sum_{n}\int_{\mathbf{k}}\overline{\eta}_{\mathbf{k},n}\frac{d}{dT}[\hat{G}_{0}^{T}(\mathbf{k},\nu_{n})]^{-1}\eta_{\mathbf{k},n}-\frac{1}{2}\Tr\left[\frac{d}{dT}[\mathbf{G}_{0}^{T}]^{-1}\left(\mathbf{\Gamma}^{(2)T}[\eta,\overline{\eta}]\right)^{-1}\right].$
(S3)

Here, we have defined

$\mathbf{G}_{0}^{T}(\mathbf{k},\nu_{n})=\left(\begin{array}[]{cc}\hat{G}_{0}^{T}(\mathbf{k},\nu_{n})&0\\\
0&-\left[\hat{G}_{0}^{T}(-\mathbf{k},-\nu_{n})\right]^{\mathrm{T}}\end{array}\right),$
(S4)

and

$\mathbf{\Gamma}^{(2)T}[\eta,\overline{\eta}](x,x^{\prime})=\left(\begin{array}[]{cc}\frac{\delta\Gamma^{T}[\eta,\overline{\eta}]}{\delta\overline{\eta}(x^{\prime})\delta\eta(x)}&\frac{\delta\Gamma^{T}[\eta,\overline{\eta}]}{\delta\overline{\eta}(x^{\prime})\overline{\delta}\eta(x)}\\\
\frac{\delta\Gamma^{T}[\eta,\overline{\eta}]}{\eta(x^{\prime})\delta\eta(x)}&\frac{\delta\Gamma^{T}[\eta,\overline{\eta}]}{\delta\eta(x^{\prime})\overline{\delta}\eta(x)}\end{array}\right),$
(S5)

where $x=(\mathbf{k},n,\sigma)$ and the symbol $[\bullet]^{\mathrm{T}}$
denotes matrix transposition and does not have to be confused with the
temperature $T$. Expanding the effective action functional
$\Gamma^{T}[\eta,\overline{\eta}]$ on both sides of (S3), one can derive an
infinite hierarchy of flow equations for the $n$ particle vertices

$\Gamma^{(2n)T}(x_{1},\dots,x_{n},x_{1}^{\prime},\dots,x_{n}^{\prime})=\frac{\delta\Gamma^{T}[\eta,\overline{\eta}]}{\delta\overline{\eta}(x_{1}^{\prime})\dots\delta\overline{\eta}(x_{n}^{\prime})\delta\eta(x_{1})\dots\delta\eta(x_{n})}\Bigg{\rvert}_{\eta,\overline{\eta}=0}.$
(S6)

Eq. (S3) is completed with an initial condition for the effective action
reading as [33]

$\Gamma^{T\to+\infty}[\eta,\overline{\eta}]=\mathcal{S}^{T\to+\infty}[\eta,\overline{\eta}].$
(S7)

In this work, we _truncate_ the effective action at the two-particle level,
that is, we set to zero all the three- and more particle correlators. We also
ignore the flow of the self-energy, that is, we keep the one-particle vertex
fixed to $(\hat{G}_{0}^{T})^{-1}$.

The only vertex function we retain is the two particle vertex

$V^{\sigma_{1}\sigma_{2}\sigma_{3}\sigma_{4},T}(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3},n_{1},n_{2},n_{3})=\Gamma^{(4)T}.$
(S8)

which depends, in principle, on the four spins of the two incoming and two
outgoing electrons, and on three momenta and Matsubara frequencies. We
introduce a further approximation by neglecting the frequency dependence.

To disentangle the dependence of $V^{T}$ on the three momenta, we employ a
channel decomposition [70, 71]:

$\begin{split}V^{\sigma_{1}\sigma_{2}\sigma_{3}\sigma_{4}}(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3})&=U\left(\delta_{\sigma_{1}\sigma_{3}}\delta_{\sigma_{2}\sigma_{4}}-\delta_{\sigma_{1}\sigma_{4}}\delta_{\sigma_{2}\sigma_{3}}\right)\\\
&\phantom{=}+P^{\sigma_{1}\sigma_{2}\sigma_{3}\sigma_{4}}_{\mathbf{k}_{1},\mathbf{k}_{3}}(\mathbf{k}_{1}+\mathbf{k}_{2})+D^{\sigma_{1}\sigma_{2}\sigma_{3}\sigma_{4}}_{\mathbf{k}_{1},\mathbf{k}_{4}}(\mathbf{k}_{3}-\mathbf{k}_{1})-D^{\sigma_{1}\sigma_{2}\sigma_{4}\sigma_{3}}_{\mathbf{k}_{1},\mathbf{k}_{3}}(\mathbf{k}_{2}-\mathbf{k}_{3}),\end{split}$
(S9)

where, to simplify the notation, we have dropped the $T$-dependence. We refer
to $P$ and $D$ as particle-particle and particle-hole channels, respectively.
It is clear from Eq. (S7) that at $T\to+\infty$ one has $V=U$ and therefore
$P=D=0$. The flow equations read as

$\displaystyle\frac{d}{dT}P^{\sigma_{1}\sigma_{2}\sigma_{3}\sigma_{4}}_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{q})=-\sum_{\begin{subarray}{c}s_{1},s_{2}\\\
s_{3},s_{4}\end{subarray}}\int_{\mathbf{p}}V^{\sigma_{1}\sigma_{2}s_{3}s_{4}}(\mathbf{k},\mathbf{q}-\mathbf{k},\mathbf{p})\,\dot{\Pi}^{s_{1}s_{2}s_{3}s_{4}}_{P,\mathbf{p}}(\mathbf{q})\,V^{s_{1}s_{2}\sigma_{3}\sigma_{4}}(\mathbf{p},\mathbf{q}-\mathbf{p},\mathbf{k}^{\prime}),$
(S10a)
$\displaystyle\frac{d}{dT}D^{\sigma_{1}\sigma_{2}\sigma_{3}\sigma_{4}}_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{q})=2\sum_{\begin{subarray}{c}s_{1},s_{2}\\\
s_{3},s_{4}\end{subarray}}\int_{\mathbf{p}}V^{\sigma_{1}s_{2}\sigma_{3}s_{4}}(\mathbf{k},\mathbf{p}+\mathbf{q},\mathbf{k}+\mathbf{q})\,\dot{\Pi}^{s_{1}s_{2}s_{3}s_{4}}_{D,\mathbf{p}}(\mathbf{q})\,V^{s_{1}\sigma_{2}s_{3}\sigma_{4}}(\mathbf{p},\mathbf{k}^{\prime}+\mathbf{q},\mathbf{p}+\mathbf{q}),$
(S10b)

with

$\displaystyle\Pi_{P,\mathbf{p}}^{s_{1}s_{2}s_{3}s_{4}}(\mathbf{q})=\sum_{n}\left[G^{T}_{0,s_{3}s_{1}}(\mathbf{p},\nu_{n})G^{T}_{0,s_{4}s_{2}}(\mathbf{q}-\mathbf{p},-\nu_{n})\right],$
(S11a)
$\displaystyle\Pi_{D,\mathbf{p}}^{s_{1}s_{2}s_{3}s_{4}}(\mathbf{q})=\sum_{n}\left[G^{T}_{0,s_{4}s_{1}}(\mathbf{p},\nu_{n})G^{T}_{0,s_{3}s_{2}}(\mathbf{p}+\mathbf{q},\nu_{n})\right].$
(S11b)

The symbol $\dot{\Pi}$ is a shorthand for $d\Pi/dT$. To simplify the
treatment, we parametrize the spin dependence of $P$ and $D$ as

$\displaystyle
P^{\sigma_{1}\sigma_{2}\sigma_{3}\sigma_{4}}_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{q})=\frac{1}{2}\sum_{\alpha,\beta=0}^{3}P^{\alpha\beta}_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{q})\,t^{\alpha}_{\sigma_{1}\sigma_{2}}\left(t^{\beta}_{\sigma_{3}\sigma_{4}}\right)^{\dagger},$
(S12a) $\displaystyle
D^{\sigma_{1}\sigma_{2}\sigma_{3}\sigma_{4}}_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{q})=\frac{1}{2}\sum_{\alpha,\beta=0}^{3}D^{\alpha\beta}_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{q})\,\tau^{\alpha}_{\sigma_{1}\sigma_{3}}\tau^{\beta}_{\sigma_{2}\sigma_{4}},$
(S12b)

where $\tau^{1,2,3}$ are the Pauli matrices, $\tau^{0}=\mathbb{1}$, and
$t^{\alpha}=i\tau^{\alpha}\tau^{2}$. To treat the dependence of the channels
on $\mathbf{k}$ and $\mathbf{k}^{\prime}$, we resort to a _truncated unity_
approach [72, 73], that is, we expand them in a complete basis of form factors
$\\{f^{\ell}_{\mathbf{k}}\\}$:

$X^{\alpha\beta}_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{q})=\sum_{\ell,\ell^{\prime}}X^{\alpha\beta}_{\ell\ell^{\prime}}(\mathbf{q})f^{\ell}_{\mathbf{k}}(f^{\ell^{\prime}}_{\mathbf{k}^{\prime}})^{*},$
(S13)

with $X=P$ or $D$, and truncate the sum up to a given order. We choose form
factors that transform in the irreducible representations of the lattice
symmetry group $C_{4v}$. It is possible to choose $f^{\ell}_{\mathbf{k}}$ in
such a way that their real space representation takes nonzero values only in a
given shell of neighbors of a given reference site [74]. It is therefore
convenient to truncate the sum (S13) by including all form factors up to a
given shell. For our calculations, we truncate the sum to the first shell of
neighbors, that is, we only consider the form factors

$\displaystyle f^{0}_{\mathbf{k}}=1,$ (S14a) $\displaystyle
f^{1}_{\mathbf{k}}=\cos k_{x}+\cos k_{y},$ (S14b) $\displaystyle
f^{2}_{\mathbf{k}}=\cos k_{x}-\cos k_{y},$ (S14c) $\displaystyle
f^{3}_{\mathbf{k}}=\sqrt{2}\sin k_{x},$ (S14d) $\displaystyle
f^{4}_{\mathbf{k}}=\sqrt{2}\sin k_{y}.$ (S14e)

Inserting the identity

$\delta_{\mathbf{k},\mathbf{k}^{\prime}}=\sum_{\ell}f^{\ell}_{\mathbf{k}}(f^{\ell}_{\mathbf{k}^{\prime}})^{*}$
(S15)

into (S10), we obtain

$\frac{d}{dT}X^{\alpha\beta}_{\ell\ell^{\prime}}(\mathbf{q})=\frac{1}{2}\zeta_{X}\sum_{\alpha^{\prime},\beta^{\prime}}\sum_{m,m^{\prime}}X[V]_{\ell
m}^{\alpha\beta^{\prime}}(\mathbf{q})\,\dot{\Pi}_{X,mm^{\prime}}^{\beta^{\prime}\alpha^{\prime}}(\mathbf{q})\,X[V]_{m^{\prime}\ell^{\prime}}^{\alpha^{\prime}\beta}(\mathbf{q}),$
(S16)

with $X=P$ or $D$, $\zeta_{P}=-1$ and $\zeta_{D}=+2$. The symbols $P[V]$ and
$D[V]$ represent the projection of the full vertex (S9) onto the $P$ and $D$
channels and they read as

$\displaystyle
P[V]_{\ell\ell^{\prime}}^{\alpha\beta}(\mathbf{q})=\int_{\mathbf{k},\mathbf{k}^{\prime}}\frac{1}{2}\sum_{\begin{subarray}{c}\sigma_{1},\sigma_{2}\\\
\sigma_{3},\sigma_{4}\end{subarray}}V^{\sigma_{1}\sigma_{2}\sigma_{3}\sigma_{4}}(\mathbf{k},\mathbf{q}-\mathbf{k},\mathbf{k}^{\prime})(t_{\sigma_{1}\sigma_{2}}^{\alpha})^{\dagger}t^{\beta}_{\sigma_{3}\sigma_{4}}(f^{\ell}_{\mathbf{k}})^{*}f^{\ell^{\prime}}_{\mathbf{k}^{\prime}},$
(S17a) $\displaystyle
D[V]_{\ell\ell^{\prime}}^{\alpha\beta}(\mathbf{q})=\int_{\mathbf{k},\mathbf{k}^{\prime}}\frac{1}{2}\sum_{\begin{subarray}{c}\sigma_{1},\sigma_{2}\\\
\sigma_{3},\sigma_{4}\end{subarray}}V^{\sigma_{1}\sigma_{2}\sigma_{3}\sigma_{4}}(\mathbf{k},\mathbf{k}^{\prime}+\mathbf{q},\mathbf{k})\tau_{\sigma_{3}\sigma_{1}}^{\alpha}\tau^{\beta}_{\sigma_{4}\sigma_{2}}(f^{\ell}_{\mathbf{k}})^{*}f^{\ell^{\prime}}_{\mathbf{k}^{\prime}}.$
(S17b)

Similar expressions hold for the bubbles
$\dot{\Pi}_{X,\ell\ell^{\prime}}^{\beta\alpha}(\mathbf{q})$. For more details
on how to perform the vertex projections, see for example Refs. [72, 73].

We run a flow starting from a high-temperature value $T_{\mathrm{ini}}\simeq
10t$ at which $P$ and $D$ can be initialized by a second order perturbation
theory result:

$X_{\ell\ell^{\prime}}^{\alpha\beta\,T_{\mathrm{ini}}}(\mathbf{q})=\frac{1}{2}\zeta_{X}U^{2}\Pi_{X,00}^{\alpha\beta\,T_{\mathrm{ini}}}(\mathbf{q})\delta_{\ell,0}\delta_{\ell^{\prime},0},$
(S18)

to (approximately) account for the integration from $T=+\infty$ to
$T=T_{\mathrm{ini}}$. We stop the flow at a temperature $T_{c}$ where one of
$P$ or $D$ exceeds the value $25t$ or at a temperature
$T_{\mathrm{min}}=10^{-2}t$ where the momentum integration in the projections
of the bubbles (S11) onto form factors becomes numerically unstable.

###  SI.2 Renormalized mean-field theory

To continue the flow below the stopping temperature
$T_{s}\equiv\mathrm{max}(T_{c},T_{\mathrm{min}})$, and therefore detect the
leading superconducting instabilities, we resort to a combination of mean-
field theory with the functional renormalization group [34, 35, 37]. This
approximation accounts for continuing the flow in the superconducting regime,
that is, by introducing a pairing gap and anomalous vertices, by keeping only
particle-particle diagrams. The simplified flow equations can be formally
integrated and cast in the form of _renormalized_ mean-field equations. In
this way, one obtains the gap equation (valid for any $0\leq T\leq T_{s}$)

$\Delta_{\mathbf{k}}^{\sigma\sigma^{\prime}}=\sum_{s,s^{\prime}}\int_{\mathbf{k}^{\prime}}\widetilde{V}^{\sigma\sigma^{\prime}ss^{\prime}}_{\mathbf{k}\mathbf{k}^{\prime}}\left[T\sum_{n}F^{ss^{\prime}}(\mathbf{k}^{\prime},\nu_{n})\right],$

where $F^{\sigma\sigma^{\prime}}(\mathbf{k},\nu_{n})$ are the spin components
of the anomalous propagator $\hat{F}(\mathbf{k},\nu_{n})$, defined through

$\mathbf{G}_{\mathrm{BdG}}(\mathbf{k},i\nu_{n})=\left(\begin{array}[]{cc}\hat{G}(\mathbf{k},\nu_{n})&\hat{F}(\mathbf{k},\nu_{n})\\\
\hat{F}^{\dagger}(\mathbf{k},\nu_{n})&-[\hat{G}(-\mathbf{k},-\nu_{n})]^{\mathrm{T}}\end{array}\right)=\left(\begin{array}[]{cc}i\nu_{n}-\hat{\mathcal{H}}_{\mathbf{k}}&-\hat{\Delta}_{\mathbf{k}}\\\
-\hat{\Delta}_{\mathbf{k}}^{\dagger}&i\nu_{n}+[\hat{\mathcal{H}}_{-\mathbf{k}}]^{\mathrm{T}}\end{array}\right)^{-1},$
(S19)

with $\hat{\Delta}_{\mathbf{k}}$ the gap matrix (in spin space). The function
$\widetilde{V}^{\sigma_{1}\sigma_{2}\sigma_{3}\sigma_{4}}_{\mathbf{k}\mathbf{k}^{\prime}}$
is the two-particle-irreducible vertex function in the particle-particle
channel at zero center-of-mass momentum. It is determined by inverting a
Bethe-Salpeter equation at the stopping temperature $T_{s}$:

$V^{\sigma_{1}\sigma_{2}\sigma_{3}\sigma_{4},T_{s}}(\mathbf{k},-\mathbf{k},\mathbf{k}^{\prime})=\widetilde{V}^{\sigma_{1}\sigma_{2}\sigma_{3}\sigma_{4}}_{\mathbf{k}\mathbf{k}^{\prime}}+\sum_{\begin{subarray}{c}s_{1},s_{2}\\\
s_{3},s_{4}\end{subarray}}\int_{\mathbf{p}}\widetilde{V}^{\sigma_{1}\sigma_{2}s_{3}s_{4}}_{\mathbf{k}\mathbf{p}}\,\Pi^{s_{1}s_{2}s_{3}s_{4},T_{s}}_{P,\mathbf{p}}(\mathbf{q}=\mathbf{0})\,V^{s_{1}s_{2}\sigma_{3}\sigma_{4},T_{s}}(\mathbf{p},-\mathbf{p},\mathbf{k}^{\prime}).$
(S20)

Eq. (S20) can be inverted with an approach similar to the one presented in the
previous section: one can paramterize the spin dependence of $\widetilde{V}$,
$V^{T_{s}}$, and $\Pi^{T_{s}}_{P}$ by means of the $t^{\alpha}$ matrices (see
Eq. (S12a)) and expand the $\mathbf{k}$\- and $\mathbf{k}^{\prime}$-dependence
in form factors and truncate the expansion at a given order. Notice that the
order kept in this inversion can be also higher than the one used in the fRG
calculation of $V^{T_{s}}$. One obtains

$P[V^{T_{s}}]^{\alpha\beta}_{\ell\ell^{\prime}}(\mathbf{q}=\mathbf{0})=2\widetilde{V}^{\alpha\beta}_{\ell\ell^{\prime}}+\sum_{\alpha^{\prime},\beta^{\prime}}\sum_{m,m^{\prime}}\widetilde{V}^{\alpha\beta^{\prime}}_{\ell
m}\,\Pi^{\beta^{\prime}\alpha^{\prime},T_{s}}_{P,mm^{\prime}}(\mathbf{q}=\mathbf{0})\,P[V^{T_{s}}]^{\alpha^{\prime}\beta}_{m^{\prime}\ell^{\prime}}(\mathbf{q}=\mathbf{0}).$
(S21)

Defining a multi-index $\mathfrak{i}=(\alpha,\ell)$, $\widetilde{V}$,
$P[V^{T_{s}}]$, and $\Pi^{T_{s}}_{P}$ become matrices, and $\widetilde{V}$ can
be cast as

$\widetilde{V}_{\mathfrak{i}\mathfrak{i}^{\prime}}=\left\\{\left[[2\mathbf{P}[V^{T_{s}}](\mathbf{q}=\mathbf{0})]^{-1}-\bm{\Pi}_{P}^{T_{s}}(\mathbf{q}=\mathbf{0})\right]^{-1}\right\\}_{\mathfrak{i}\mathfrak{i}^{\prime}}.$
(S22)

By expanding also the gap in form factors and $t^{\alpha}$ matrices, the gap
equation takes the form

$\Delta^{\alpha}_{\ell}=\sum_{\beta,\ell^{\prime}}\widetilde{V}^{\alpha\beta}_{\ell\ell^{\prime}}\mathcal{F}_{\ell^{\prime}}^{\beta},$
(S23)

where we have defined

$\Delta_{\mathbf{k}}^{\sigma\sigma^{\prime}}=\sum_{\alpha=0}^{3}\sum_{\ell}\Delta^{\alpha}_{\ell}t^{\alpha}_{\sigma\sigma^{\prime}}f^{\ell}_{\mathbf{k}},$
(S24)

and

$\mathcal{F}_{\ell^{\prime}}^{\beta}=\frac{1}{2}T\sum_{n}\int_{\mathbf{k}}\Tr\left[F(\mathbf{k},\nu_{n})(t^{\alpha})^{\dagger}\right](f^{\ell}_{\mathbf{k}})^{*}.$
(S25)

Notice that Eq. (S23) still depends on the temperature via
$\mathcal{F}^{\alpha}_{\ell}$.

Eq. (S23) can be used to determine the gap function at a given temperature
$T\geq T_{s}$ but also to simply determine the leading superconducting state.
This can be achieved by _linearizing_ it with respect to
$\Delta^{\alpha}_{\ell}$. By noticing that

$\frac{\delta\mathcal{F}^{\alpha}_{\ell}}{\delta\Delta^{\beta}_{\ell^{\prime}}}\bigg{\rvert}_{\Delta^{\alpha}_{\ell}=0}=\Pi_{P,\ell\ell^{\prime}}^{\alpha\beta}(\mathbf{q}=\mathbf{0}),$
(S26)

we obtain

$\Delta^{\alpha}_{\ell}=\sum_{\beta,\gamma}\sum_{\ell^{\prime},m}\widetilde{V}^{\alpha\beta}_{\ell\ell^{\prime}}\Pi_{P,\ell^{\prime}m}^{\beta\gamma}(\mathbf{q}=\mathbf{0})\Delta^{\gamma}_{m}.$
(S27)

For temperatures close to the superconducting transition temperature
$T_{\mathrm{SC}}$, where the overall magnitude of the gap is small and the
linearization carried out above is justifiable, the gap is an eigenvector of
the matrix $\widetilde{V}\circ\Pi_{P}(\mathbf{q}=\mathbf{0})$ with eigenvalue
1. $T_{\mathrm{SC}}$ can therefore be determined as the temperature at which
the largest positive eigenvalue of the matrix
$\widetilde{V}\circ\Pi_{P}(\mathbf{q}=\mathbf{0})$ becomes 1. Since in a Fermi
liquid some components of $\Pi_{P}(\mathbf{q}=\mathbf{0})$ diverge like
$\log(T^{-1})$ (or $\log^{2}(T^{-1})$ at a van Hove singularity) as $T$
approaches 0, one will always find a finite value of $T_{\mathrm{SC}}$, or, in
other words, Eq. (S23) has always a nontrivial solution at $T=0$. The
numerical values of $\Delta^{\alpha}_{\ell}$, however, can be exponentially
small, so that their calculation becomes numerically involved. For this
reason, we analyze the linearized gap equation, Eq. (S27), at a low
temperature, $T_{\mathrm{low}}=10^{-6}t$, and determine the leading
superconducting eigenvalue $\lambda_{\mathrm{SC}}$, defined as the largest
positive eigenvalue of $\widetilde{V}\circ\Pi_{P}(\mathbf{q}=\mathbf{0})$. The
symmetry of the corresponding eigenvector determines the symmetry of the
leading superconducting instability.

##  SII Symmetry classification of the superconducting states

| $\mathbf{k}$ | $\vec{\tau}$ | $\hat{U}$
---|---|---|---
E |  $k_{x}\to k_{x}$ $k_{y}\to k_{y}$ |  $\tau^{1}\to\tau^{1}$ $\tau^{2}\to\tau^{2}$ $\tau^{3}\to\tau^{3}$ | $\mathbb{1}$
$\mathrm{C}_{4}^{\pm}$ |  $k_{x}\to\mp k_{y}$ $k_{y}\to\pm k_{x}$ |  $\tau^{1}\to\mp\tau^{2}$ $\tau^{2}\to\pm\tau^{1}$ $\tau^{3}\to\tau^{3}$ | $e^{\pm i\frac{\pi}{4}\tau^{3}}$
$\mathrm{C}_{2}$ |  $k_{x}\to-k_{x}$ $k_{y}\to-k_{y}$ |  $\tau^{1}\to-\tau^{1}$ $\tau^{2}\to-\tau^{2}$ $\tau^{3}\to\tau^{3}$ | $e^{i\frac{\pi}{2}\tau^{3}}$
$\Sigma_{v}^{x(y)}$ |  $k_{x}\to\pm k_{x}$ $k_{y}\to\mp k_{y}$ |  $\tau^{1}\to\mp\tau^{1}$ $\tau^{2}\to\pm\tau^{2}$ $\tau^{3}\to\tau^{3}$ |  $i\tau^{3}\mathcal{K}$ ($\mathbb{1}\mathcal{K}$)
$\Sigma_{d}^{1(2)}$ |  $k_{x}\to\pm k_{y}$ $k_{y}\to\pm k_{x}$ |  $\tau^{1}\to\mp\tau^{2}$ $\tau^{2}\to\mp\tau^{1}$ $\tau^{3}\to\tau^{3}$ | $e^{\pm i\frac{\pi}{4}\tau^{3}}\mathcal{K}$
Table SI: Discrete spin-lattice transformations that leave the Rashba-Hubbard Hamiltonian invariant. In the first column we show their effect on the lattice momentum $\mathbf{k}$, in the second one how they transform the Pauli matrices $\vec{\tau}$, and in the third one the operator that implements the spin part of the transformation on the electron operators. We remark that the transformation of the Pauli matrices is implemented as $\vec{\tau}\to\hat{U}^{\dagger}\vec{\tau}\hat{U}$. By $\mathcal{K}$ we denote the complex conjugation operator. | E | $\mathrm{C}_{4}^{\pm}$ | $\mathrm{C}_{2}$ | $\Sigma_{v}^{x(y)}$ | $\Sigma_{d}^{1(2)}$ | $\hat{\Delta}_{\mathbf{k}}$
---|---|---|---|---|---|---
$A_{1}$ | 1 | 1 | 1 | 1 | 1 | $g^{A_{1}}_{\mathbf{k}}t^{0}-g^{E,y}_{\mathbf{k}}t^{1}+g^{E,x}_{\mathbf{k}}t^{2}$
$A_{2}$ | 1 | 1 | 1 | -1 | -1 | $g^{A_{2}}_{\mathbf{k}}t^{0}+g^{E,x}_{\mathbf{k}}t^{1}+g^{E,y}_{\mathbf{k}}t^{2}$
$B_{1}$ | 1 | -1 | 1 | 1 | -1 | $g^{B_{1}}_{\mathbf{k}}t^{0}+g^{E,y}_{\mathbf{k}}t^{1}+g^{E,x}_{\mathbf{k}}t^{2}$
$B_{2}$ | 1 | -1 | 1 | -1 | 1 | $g^{B_{2}}_{\mathbf{k}}t^{0}+g^{E,x}_{\mathbf{k}}t^{1}-g^{E,y}_{\mathbf{k}}t^{2}$
$E$ | $\begin{pmatrix}1&0\\\ 0&1\end{pmatrix}$ | $\begin{pmatrix}0&\pm 1\\\ \mp 1&0\end{pmatrix}$ | $\begin{pmatrix}-1&0\\\ 0&-1\end{pmatrix}$ | $\begin{pmatrix}\pm 1&0\\\ 0&\mp 1\end{pmatrix}$ | $\begin{pmatrix}0&\pm 1\\\ \pm 1&0\end{pmatrix}$ | $\left(g^{E,x}_{\mathbf{k}}t^{3},\,g^{E,y}_{\mathbf{k}}t^{3}\right)$
Table SII: Transformation properties of the irreps under the action of the spin-lattice symmetry group. In the two-dimensional irrep $E$ the transformation matrices depend on the basis choice. In the last column, we show the form of the gap function in each irrep. The functions $g^{A_{1}}_{\mathbf{k}}$, $g^{A_{2}}_{\mathbf{k}}$, $g^{B_{1}}_{\mathbf{k}}$ and $g^{B_{2}}_{\mathbf{k}}$ are irreps of the sole $C_{4v}$, while $g^{E,x}_{\mathbf{k}}$ and $g^{E,y}_{\mathbf{k}}$ are a basis of the two-dimensional irrep $E$ of $C_{4v}$ chosen such that they transform under $C_{4v}$ as shown in the last line of the table. The lowest harmonic contributions for these functions are listed in Table SIII. $g^{A_{1}}_{\mathbf{k}}$ | 1, $\cos k_{x}+\cos k_{y}$, $2\cos k_{x}\cos k_{y}$, $\cos(2k_{x})+\cos(2k_{y})$, …
---|---
$g^{A_{2}}_{\mathbf{k}}$ |  $2\sqrt{2}\left(\cos k_{x}-\cos k_{y}\right)\sin k_{x}\sin k_{y}$, …
$g^{B_{1}}_{\mathbf{k}}$ |  $\cos k_{x}-\cos k_{y}$, $\cos(2k_{x})-\cos(2k_{y})$, …
$g^{B_{2}}_{\mathbf{k}}$ |  $2\sin k_{x}\sin k_{y}$, $2\sqrt{2}\left(\cos k_{x}+\cos k_{y}\right)\sin k_{x}\sin k_{y}$, …
$\left(g^{E,x}_{\mathbf{k}},g^{E,y}_{\mathbf{k}}\right)$ |  $\left(\sqrt{2}\sin k_{x},\sqrt{2}\sin k_{y}\right)$, $\left(2\cos k_{y}\sin k_{x},2\cos k_{x}\sin k_{y}\right)$, …
Table SIII: Lowest harmonics for the functions $g^{A_{1}}_{\mathbf{k}}$,
$g^{A_{2}}_{\mathbf{k}}$, $g^{B_{1}}_{\mathbf{k}}$, $g^{B_{2}}_{\mathbf{k}}$,
$g^{E,x}_{\mathbf{k}}$, $g^{E,y}_{\mathbf{k}}$ shown in Table SII. Note that
we have employed the normalization convention
$\int_{\mathbf{k}}|g_{\mathbf{k}}^{X}|^{2}=1$.

The symmetry classification of the superconducting symmetries relies on the
symmetry group of the Hamiltonian (1). Because of the Rashba spin-orbit
coupling, the lattice point group $C_{4v}$ does not leave the Hamiltonian
invariant. We have therefore to consider a different group that combines
$C_{4v}$ with spin rotations. This group can be simply obtained from the
square lattice point group by adding to each of its elements a spin
transformation that leaves the Rashba term $\sim(-\sin k_{y}\,\tau^{1}+\sin
k_{x}\,\tau^{2})$ invariant. In Table SI we list all the discrete combined
spin-lattice transformations that leave the Hamiltonian (1) invariant. These
symmetries form a discrete group with eight elements
$\\{\mathrm{E},\mathrm{C}_{4}^{+},\mathrm{C}_{4}^{-},\mathrm{C}_{2},\Sigma_{v}^{x},\Sigma_{v}^{y},\Sigma_{d}^{1},\Sigma_{d}^{2}\\}$.
Since this group is equivalent to the point group $C_{4v}$, its irreducible
representations (irreps) are the same: four are one-dimensional ($A_{1}$,
$A_{2}$, $B_{1}$, $B_{2}$), and one two-dimensional ($E$). We can therefore
classify the symmetry of the superconducting gap according to these irreps. We
decompose the gap into form factors and $t^{\alpha}$ matrices, according to
Eq. (S24). We notice that the action of the group on the matrices $t^{\alpha}$
is given by

$t^{\alpha}\to\hat{U}^{\dagger}t^{\alpha}\hat{U}^{*},$ (S28)

with the different transformation matrices $\hat{U}$ listed in Table SI. As a
consequence, when $\alpha=1,2,3$, the matrices $t^{\alpha}$ transform exactly
as the Pauli matrices (as shown in Table SI), while $t^{0}$ transforms always
trivially. We also remark that, in order to preserve the overall antisymmetry
of the pairing wavefunction, the singlet (triplet) components of the gap
function $\Delta^{\alpha}_{\ell}$, corresponding to $\alpha=0$
($\alpha=1,2,3$), must be symmetric (antisymmetric) for
$\mathbf{k}\to-\mathbf{k}$. For this reason, only the pairs of indices
$(\alpha,\ell)$ with $\alpha=0$ and $f^{\ell}_{\mathbf{k}}$ even or
$\alpha=1,2,3$ and $f^{\ell}_{\mathbf{k}}$ odd are allowed in
$\Delta^{\alpha}_{\ell}$. In Table SII, we list the irreps of the combined
lattice-spin discrete symmetry group, together with the typical form that the
gap function $\hat{\Delta}_{\mathbf{k}}$ takes in each one of them. In Table
SIII we show the lowest harmonics for the functional dependence of
$\hat{\Delta}_{\mathbf{k}}$ on $\mathbf{k}$. Note that equal spin triplet
pairing (described with the matrices $t^{1}$ and $t^{2}$) can only occur in
combination with singlet pairing (described by $t^{0}$) in a one-dimensional
representation. On the other hand, opposite spin triplet pairing ($t^{3}$) can
take place only in the two-dimensional irrep $E$.

##  SIII Real space renormalized mean-field

In Sec. SI.2, we solve the renormalized mean-field Hamiltonian in momentum
space due to a translational invariance. In this section, we give the details
of solving the renormalized mean-field Hamiltonian in real space. In the
absence of any inhomogeneities, the two approaches are analogues and give the
same results. However, inhomoegenieties arising either due to edges/surfaces
or disorder can only be captured in a real space picture. Here, we ignore any
disorder and focus on the edge properties and possible symmetry breakings
within the real space picture. We decouple $H_{U}$ in Eq. (1) only in the
Cooper channel and the resultant real space Hamiltonian is,

$H_{\rm eff}=H_{0}+\sum_{\begin{subarray}{c}i\mu\\\
\sigma\sigma^{\prime}\end{subarray}}\Delta_{\mu}^{\sigma\sigma^{\prime}}(i)c^{\dagger}_{i+\mu,\sigma^{\prime}}c^{\dagger}_{i,\sigma}+\textrm{H.c.},$
(S29)

where $\mu$ denotes the nearest neighbor bonds $\pm x,\pm y$ of the lattice
site $i$, $H_{0}$ is defined in Eq. (1) and
$\Delta_{\mu}^{\sigma\sigma^{\prime}}$ is the superconducting order parameter
given by

$\Delta_{\mu}^{\sigma\sigma^{\prime}}(i)=\sum_{s,s^{\prime}}\sum_{\nu}\widetilde{V}^{\sigma\sigma^{\prime}ss^{\prime}}_{\mu\nu}\langle
c_{i,s}c_{i+\nu,s^{\prime}}\rangle,$ (S30)

where the interaction strengths are

$\widetilde{V}^{\sigma\sigma^{\prime}ss^{\prime}}_{\mu\nu}=\frac{1}{2}\sum_{\alpha\beta}\sum_{\ell\ell^{\prime}}\widetilde{V}^{\alpha\beta}_{\ell\ell^{\prime}}t^{\alpha}_{\sigma\sigma^{\prime}}(t^{\beta}_{ss^{\prime}})^{\dagger}f^{\ell}(\mu)[f^{\ell^{\prime}}(\nu)]^{*},$
(S31)

with $f^{\ell}(\mu)$ being the Fourier transform of $f^{\ell}_{\mathbf{k}}$
and $\widetilde{V}^{\alpha\beta}_{\ell\ell^{\prime}}$ are obtained from the
fRG calculation. $H_{\rm MF}$ is diagnolized using the Bogoliubov-de Gennes
transformations [75],
$c_{i\sigma}=\sum_{n}(\gamma_{n}u_{i\sigma}^{n}-\sigma\gamma^{\dagger}_{n}{v_{i\sigma}^{n}}^{*})$,
where $\gamma^{\dagger}_{n\sigma}$ ($\gamma_{n\sigma}$) are the creation
(annihilation) operators of the Bogoliubov quasiparticle at a state $n$, and
$u_{i\sigma}^{n}$ and $v_{i\sigma}^{n}$ the eigenfunctions with eigenvalues
$E$. The resulting eigen-system is then solved self-consistently for the
superconducting order parameters $\Delta_{\mu}^{\sigma\sigma^{\prime}}$.
Moreover, the findings of the order paramters obtained in the momentum space
suggest that the dominant pairing is only on nearest neighbors, i.e., either
$A_{1}$ or $B_{1}$ phases, for most parts of the phase diagram except in a
narrow region of the phase diagram ($B_{2}$ phase red region in Fig. 1).
Hence, we restrict the order parameters only to nearest neighbors in real
space, as already expressed in Eq. (S29). This helps to reduce the numer of
independent self-consistent order parameters to 16 per lattice site. We then
project local $\Delta_{\mu}^{\sigma\sigma^{\prime}}$ on onsite equivalents of
different pairing components to analyse the symmetry of the order parameter.
As shown in Tables. SII and SIII, the only possible nearest neighbour singlet
components are $d$-wave and extended $s$-wave order parameters. The onsite
equivalents of the singlets are given as

$\Delta_{d/s}(i)=\frac{1}{8}\sum_{\mu}\epsilon^{d/s}_{\mu}\sum_{\sigma}\mathrm{sgn}(\sigma)\Delta_{\mu}^{\sigma\overline{\sigma}}(i),$
(S32)

where $\Delta_{d}$ is the $d$-wave order parameter, which in momentum space is
proportional to $(\cos k_{x}-\cos k_{y})$. $\Delta_{s}$ is the extended
$s$-wave order parameter, which in momentum space is proportional to $(\cos
k_{x}+\cos k_{y})$. The perfactor $\epsilon^{d/s}_{\mu}$ is
$\epsilon^{d/s}_{\mu}=1$ for $\mu=\pm x$ and $\epsilon^{d/s}_{\mu}=-1/1$ for
$\mu=\pm y$. The phase of the $d$-wave order parameter is characterized with
$\theta$ with the definition
$\Delta_{d}(i)=\left|\Delta_{d}(i)\right|\exp(i\theta)$. The triplet order
parameters are mainly $p$-wave with the onsite equivalent given as

$\Delta_{p_{\nu}}^{\sigma\sigma^{\prime}}(i)=\frac{1}{2i}\sum_{\mu}\mathrm{sgn}(\mu)\Delta_{\mu}^{\sigma\sigma^{\prime}},$
(S33)

where $\Delta_{p_{x}/p_{y}}^{\sigma\sigma^{\prime}}$ are the $p$-wave order
parameters, which in momentum space is proportional to $\sin k_{x}$ or $\sin
k_{y}$. Different spin components of the $p$-wave order parameters can be
combined to obtain the order parameters in the irreducible represenations as

$\displaystyle\Delta_{A_{1}}^{t}(i)=\frac{1}{4}\left(\sum_{\sigma}\mathrm{sgn}(\sigma)\Delta_{p_{y}}^{\sigma\sigma}-i\sum_{\sigma}\Delta_{p_{x}}^{\sigma\sigma}\right),$
(S34)
$\displaystyle\Delta_{B_{1}}^{t}(i)=\frac{1}{4}\left(\sum_{\sigma}\mathrm{sgn}(\sigma)\Delta_{p_{y}}^{\sigma\sigma}+i\sum_{\sigma}\Delta_{p_{x}}^{\sigma\sigma}\right).$
(S35)

Other mixed-spin triplet components are negligible compared to the equal spin
triplets.

The interaction strengths
$\widetilde{V}^{\sigma\sigma^{\prime}ss^{\prime}}_{\mu\nu}$ obtained from fRG
give rise to very small values of the superconducting order parameters in most
parts of the phase diagram. As a result, the corresponding coherence lengths
are large giving the necessity of working with extremely large system sizes to
avoid interference effects of the edges while studying the topological edge
states. Hence, we scale the interaction strengths obtained from fRG by 10 in
order to solve the gap equation Eq. (S30) in real space for all representative
points in the phase diagram investigated. We have verified that the
qualitative features do not change with reducing the scale factor. We first
match the self-consistent superconducting order parameters with periodic
boundary condition using the solutions of the momentum space and real space
gap equations. We then use open boundary conditions to investigate each of the
phases. The choice of the geometry of the boundary depends on the
characteristic edge states of different superconducting phases, as explained
in the main text.

##  SIV effective model around Van Hove singularity points

In this section we derive low-energy Hamiltonians near the two VHS points of
the square lattice and discuss the logarithmic divergences in the DOS near
these VHSs. The non-interacting Hamiltonian matrix of the Rashba model on the
square lattice is given by

$\displaystyle
h(\mathbf{k})=(\epsilon_{\mathbf{k}}-\mu)\tau_{0}+\lambda\vec{g}_{\mathbf{k}}\cdot\vec{\tau},$
(S36)

where $\epsilon_{\mathbf{k}}=-2t(\cos k_{x}+\cos k_{y})-4t^{\prime}\cos
k_{x}\cos k_{y}$ and $\vec{g}_{\mathbf{k}}=2t(-\sin k_{y},\sin k_{x},0)$. The
eigenvalues are $E_{\mathbf{k}}=(\epsilon_{\mathbf{k}}-\mu)\pm
2t\lambda\sqrt{\sin^{2}k_{x}+\sin^{2}k_{y}}$. The chemical potentials of the
two VHS points are $\mu_{\textrm{VHS}}=\pm 2(\sqrt{(t\pm
2t^{\prime})^{2}+t^{2}\lambda^{2}}-t)$. The low-energy Hamiltonian around the
VHS point ($k_{c}$, $0$) on the $\Gamma-X$ direction (VHS1) reads

$\displaystyle\epsilon(k_{c}+q_{x},q_{y})$ $\displaystyle=$
$\displaystyle-2\left[(\cos k_{c}-\sin k_{c}\lambda+1)+2\cos
k_{c}t^{\prime}\right]+q_{x}^{2}\left[t(\cos k_{c}-\sin k_{c}\lambda)+2\cos
k_{c}t^{\prime}\right]$ $\displaystyle+$ $\displaystyle q_{y}^{2}\left(2\cos
k_{c}t^{\prime}+\frac{|\lambda|t}{\sin k_{c}}+t\right),$ $\displaystyle=$
$\displaystyle-2\left[t(\cos k_{c}-\sin k_{c}\lambda+1)+2\cos
k_{c}t^{\prime}\right]-q_{x}^{2}\sin
k_{c}\frac{(t+2t^{\prime})^{2}+t^{2}\lambda^{2}}{\lambda
t}+q_{y}^{2}\frac{2\sin k_{c}\cos k_{c}t^{\prime}+t\sin k_{c}+|\lambda|t}{\sin
k_{c}},$

where $k_{c}$ is given by $\tan k_{c}=-\frac{\lambda t}{t+2t^{\prime}}$. The
dispersion can be further simplified to

$\displaystyle\epsilon(k_{c}+q_{x},q_{y})$ $\displaystyle=$ $\displaystyle
2(D_{1}-t)-q_{x}^{2}D_{1}+q_{y}^{2}\left[-\frac{2(t+2t^{\prime})t^{\prime}}{D_{1}}+D_{1}+t\right],$
(S38)

with $D_{1}=\sqrt{(t+2t^{\prime})^{2}+\lambda^{2}t^{2}}$.

The low-energy Hamiltonian around the VHS point ($\pi$, $k^{\prime}_{c}$) on
the $X-M$ direction (VHS2) reads

$\displaystyle\epsilon(\pi+q_{x},k^{\prime}_{c}+q_{y})$ $\displaystyle=$
$\displaystyle 2\left[t(-\cos k^{\prime}_{c}-\sin k_{c}\lambda+1)+2\cos
k^{\prime}_{c}t^{\prime}\right]+q_{y}^{2}\left[t(\cos k^{\prime}_{c}+\sin
k^{\prime}_{c}\lambda)-2\cos k^{\prime}_{c}t^{\prime}\right]$ (S39)
$\displaystyle-$ $\displaystyle q_{x}^{2}\left(2\cos
k^{\prime}_{c}t^{\prime}+\frac{|\lambda|t}{\sin k^{\prime}_{c}}+t\right)$
$\displaystyle=$ $\displaystyle 2\left[t(-\cos k^{\prime}_{c}-\sin
k_{c}\lambda+1)+2\cos k^{\prime}_{c}t^{\prime}\right]+q_{y}^{2}\sin
k_{c}\frac{(t-2t^{\prime})^{2}+t^{2}\lambda^{2}}{\lambda t}$ $\displaystyle-$
$\displaystyle q_{x}^{2}\left(2\cos
k^{\prime}_{c}t^{\prime}+\frac{|\lambda|t}{\sin k^{\prime}_{c}}+t\right),$

where $k^{\prime}_{c}$ is given by $\tan k^{\prime}_{c}=\frac{\lambda
t}{t-2t^{\prime}}$. The dispersion can be further simplified to

$\displaystyle\epsilon(\pi+q_{x},q_{y})$ $\displaystyle=$
$\displaystyle-2(D_{2}-t)+q_{y}^{2}D_{2}-q_{x}^{2}\left[\frac{2(t-2t^{\prime})t^{\prime}}{D_{2}}+D_{2}+t\right],$
(S40)

with $D_{2}=\sqrt{(t-2t^{\prime})^{2}+\lambda^{2}t^{2}}$. The two VHSs are
derived from the saddle at the X/Y point with Rashba spin-orbit coupling. In
the absence of the next-nearest hopping $t^{\prime}$, the dispersion around
these two VHSs are similar and the corresponding effective masses along $x$
and $y$ are the same, leading to the same divergent density of states (see
blue curves in Fig. S1). With increasing $|t^{\prime}|$ (negative
$t^{\prime}$), the VHS1 on the $\Gamma-X$ line rapidly moves away from the X
point, while the VHS2 on the $X-M$ line slightly move closer to the X point.
As the coefficient of the logarithmically divergent DOS is proportional to
$\sqrt{m_{x}m_{y}}$, where $m_{x/y}$ is the effective mass along $x/y$
direction, this will make the DOS around the two VHSs different, which can be
seen in Fig. S1. The saddle points on the $X-M$ lines (VHS2) are close to the
X point and this will generate strong ferromagnetic spin fluctuations, which
promoting spin-triplet pairing. This explains the dominant spin-triplet
pairing near the VHS2 filling ($|t^{\prime}|>0.3$) in our fRG calculations.

Figure S1:  Band structures (a) and density of states (DOS) (b) around two the
VHSs of the square lattice as a function of the second neighbor hopping
$t^{\prime}$.

##  SV topological edge and corner states

In this section we study the topological edge and corner states of the $A_{1}$
and $A_{1}+iB_{1}$ superconducting phases. With the spinor
$\Psi_{\mathbf{k}}=(c^{\dagger}_{\mathbf{k}\uparrow},c^{\dagger}_{\mathbf{k}\downarrow},c_{\bm{-k}\uparrow},c_{\mathbf{k}\downarrow})$
the BdG Hamiltonian can be written as
$\mathcal{H}_{\textrm{BdG}}=\sum_{\mathbf{k}}\Psi^{\dagger}_{\mathbf{k}}H_{\textrm{BdG}}(\mathbf{k})\Psi_{\mathbf{k}}$
with the Hamiltonian matrix

$\displaystyle
H_{\textrm{BdG}}(\mathbf{k})=\left(\begin{array}[]{cc}h(\mathbf{k})&\Delta(\mathbf{k})\\\
\Delta^{\dagger}(\mathbf{k})&-h^{\star}(\mathbf{k})\\\ \end{array}\right),$
(S43)

where the pairing term is
$\Delta(\mathbf{k})=[\Delta_{S}(\mathbf{k})+\Delta_{T}(\mathbf{k})\vec{d}\cdot\vec{\tau}]i\tau_{2}$,
with $\Delta_{S}$($\Delta_{T}$) being the spin-singlet (spin-triplet) pairing
component. For the helical $p$-wave ($A_{1}$ pairing state), the pairing
matrix is

$\displaystyle\Delta^{A_{1}}_{T}\vec{d}=\Delta_{p}(-\sin k_{y},\sin k_{x}),$
(S44)

with $\vec{d}$ being parallel to $\vec{g}$
($\vec{d}_{\mathbf{k}}\parallel\vec{g}_{\mathbf{k}}$). The singlet pairing
component is

$\displaystyle\Delta_{S}(\mathbf{k})=\Delta_{0}+\Delta_{s}(\mathbf{k})+\Delta_{d}(\mathbf{k})=\Delta_{0}+\Delta_{s}(\cos
k_{x}+\cos k_{y})+\Delta_{d}(\cos k_{x}-\cos k_{y}).$ (S45)

The BdG Hamiltonian can be written as,

$\displaystyle
H_{\textrm{BdG}}(\mathbf{k})=\epsilon_{\mathbf{k}}\eta_{z}\tau_{0}+g_{x}(\mathbf{k})\eta_{0}\tau_{x}+g_{y}(\mathbf{k})\eta_{z}\tau_{y}-\Delta_{p}[d_{x}(\mathbf{k})\eta_{x}\tau_{z}+d_{y}(\mathbf{k})\eta_{y}\tau_{z}]-\Delta_{S}(\mathbf{k})\eta_{y}\tau_{y},$
(S46)

where $\vec{\eta}$ and $\vec{\tau}$ are Pauli matrices in the Nambu and spin
space, respectively. This Hamiltonian has time-reversal symmetry, particle-
hole symmetry and chiral symmetry, satisfying
$\mathcal{T}H_{\textrm{B}dG}(\mathbf{k})\mathcal{T}^{-1}=H_{\textrm{BdG}}(-\mathbf{k})$
with $\mathcal{T}=i\tau_{2}\mathcal{K}$,
$\mathcal{C}H_{\textrm{BdG}}(\mathbf{k})\mathcal{C}^{-1}=-H_{\textrm{BdG}}(-\mathbf{k})$
with $\mathcal{C}=\eta_{1}\mathcal{K}$ and
$\mathcal{S}H_{\textrm{BdG}}(\mathbf{k})\mathcal{S}^{-1}=-H_{\textrm{BdG}}(\mathbf{k})$
with $\mathcal{S}=\mathcal{T}\mathcal{C}=i\eta_{1}\tau_{2}$. Owing to
$\mathcal{T}^{2}=-1$, $\mathcal{C}^{2}=+1$ and $\mathcal{S}^{2}=1$, the system
belongs to the DIII class of topological superconductors [3] and is
characterized by a $\mathbb{Z}_{2}$ invariant in two dimensions [76].

To study the edge states, we expand the above model around the $\Gamma$ point.
(The expansion around the M point is briefly discussed at the end.) The
effective model around the $\Gamma$ point reads,

$\displaystyle
H^{\Gamma}_{\textrm{BdG}}(\bm{q})=[M_{0}+B(q^{2}_{x}+q^{2}_{y})]\eta_{z}\tau_{0}-Aq_{y}\eta_{0}\tau_{x}+Aq_{x}\eta_{z}\tau_{y}-\Delta_{p}(-q_{y}\eta_{x}\tau_{z}+q_{x}\eta_{y}\tau_{z})-\Delta_{S}(\bm{q})\eta_{y}\tau_{y}.$
(S47)

Here, the parameters are given by

$\displaystyle M_{0}=-4t-4t^{\prime}-\mu,\qquad B=t+2t^{\prime},\qquad
A=t\lambda,$
$\displaystyle\textrm{and}\quad\Delta_{S}=\Delta_{0}+2\Delta_{s}+\frac{1}{2}\Delta_{s}\left(q^{2}_{x}+q^{2}_{y}\right)+\frac{1}{2}\Delta_{d}\left(q^{2}_{x}-q^{2}_{y}\right).$
(S48)

We first omit the singlet pairing and study the topological edge states for
the helical $p$-wave pairing. Considering open boundary conditions along the
$x$ direction, we let $q_{x}\rightarrow-i\partial_{x}$ in the Hamiltonian
(S47) and obtain

$\displaystyle H^{\Gamma}_{\textrm{BdG}}(-i\partial_{x},q_{y})$
$\displaystyle=$ $\displaystyle
H^{\Gamma}_{0}(-i\partial_{x},q_{y})+H^{\Gamma}_{1}(-i\partial_{x},q_{y}),$
$\displaystyle H^{\Gamma}_{0}(-i\partial_{x},q_{y})$ $\displaystyle=$
$\displaystyle(M_{0}-B\partial^{2}_{x})\eta_{z}\tau_{0}-iA\eta_{z}\tau_{y}\partial_{x}+\Delta_{p}\eta_{y}\tau_{0}\partial_{x},$
$\displaystyle H^{\Gamma}_{1}(-i\partial_{x},q_{y})$ $\displaystyle=$
$\displaystyle-Aq_{y}\eta_{0}\tau_{x}+\Delta_{p}q_{y}\eta_{x}\tau_{z},$ (S49)

where we keep the dominant linear $q_{y}$ terms but neglect the insignificant
$q^{2}_{y}$ terms. In the following, we solve the Hamiltonian
$H^{\Gamma}_{0}(-i\partial_{x},q_{y})$ and treat
$H^{\Gamma}_{1}(-i\partial_{x},q_{y})$ perturbatively. As the system belongs
to DIII class, zero-energy localized edge modes will appear in the
topologically nontrivial regime. The eigenvalue equation for these zero modes
reads $H^{\Gamma}_{0}(-i\partial_{x},q_{y})\psi(x)=0$ with the boundary
condition $\psi(0)=\psi(+\infty)=0$. With the wavefunction ansatz
$\psi(x)=e^{-\alpha x}\phi$, the eigenvalue equation is further simplified to

$\displaystyle\left[(M_{0}-B\alpha^{2})\eta_{z}\tau_{0}+iA\alpha\eta_{z}\tau_{y}-i\Delta_{p}\alpha\eta_{y}\tau_{0}\right]\phi=0.$
(S50)

By multiplying $\eta_{z}\tau_{0}$ from the right hand side, the above equation
can be brought into the form

$\displaystyle[(M_{0}-B\alpha^{2})+iA\alpha\tau_{y}-\Delta_{p}\alpha\eta_{x}]\phi=0.$
(S51)

Therefore, the four-component wavefunction $\phi$ is an eigenstate of both
$\sigma_{y}$ and $\eta_{x}$, and the corresponding eigenbasis is

$\displaystyle\phi_{1}$ $\displaystyle=$
$\displaystyle|\eta_{x}=1\rangle\otimes|\tau_{y}=1\rangle,$
$\displaystyle\phi_{2}$ $\displaystyle=$
$\displaystyle|\eta_{x}=1\rangle\otimes|\tau_{y}=-1\rangle,$
$\displaystyle\phi_{3}$ $\displaystyle=$
$\displaystyle|\eta_{x}=-1\rangle\otimes|\tau_{y}=1\rangle,$
$\displaystyle\phi_{4}$ $\displaystyle=$
$\displaystyle|\eta_{x}=-1\rangle\otimes|\tau_{y}=-1\rangle.$ (S52)

For the four $\phi_{i}$ of the eigenbasis we obtain the equation
$B\alpha^{2}+(-iAf^{i}_{\sigma}+\Delta_{p}f^{i}_{\eta})\alpha-M_{0}=0$, with
$f^{i}_{\sigma/\eta}=\pm$, which has the following roots

$\displaystyle\alpha_{1,2}=\frac{(-\Delta_{p}f^{i}_{\eta}+iAf^{i}_{\sigma})\pm\sqrt{(-\Delta_{p}f^{i}_{\eta}+iAf^{i}_{\sigma})^{2}+4BM_{0}}}{2B}.$
(S53)

The boundary condition $\psi(+\infty)=0$ implies $\text{Re}(\alpha)>0$. With
$\alpha_{1}+\alpha_{2}=(-\Delta_{p}f^{i}_{\eta}+iAf^{i}_{\sigma})/B$ and
taking $A,B>0$, only $\phi_{3,4}$ satisfies the above constrain with
$f^{i}_{\eta}=-1$. As the pairing potential $\Delta_{p}$ is small, we omit
$A^{2}/B^{2}$ and $A/B^{2}$ terms in the following calculations. With
including the boundary condition $\psi(0)=0$, the two wavefunctions can be
written as

$\displaystyle\psi_{3/4}(x)=\mathcal{N}_{0}e^{-\beta_{0}x}e^{\mp
i\beta_{1}x}\sin(\beta_{2}x)\phi_{3/4},$ (S54)

with $\beta_{0}=\Delta_{p}/2B$, $\beta_{1}=A/2B$,
$\beta_{2}=\sqrt{-4BM_{0}+A^{2}}/2B$ and
$\mathcal{N}_{0}=\sqrt{\frac{4\beta_{0}(\beta^{2}_{0}+\beta^{2}_{2})}{\beta^{2}_{2}}}$.
In the basis of $\psi_{3,4}(x)$, the elements of
$H^{\Gamma}_{1}(-i\partial_{x},q_{y})$ can be written as

$\displaystyle\langle\psi_{3}|H^{\Gamma}_{1}(-i\partial_{x},q_{y})|\psi_{3}\rangle$
$\displaystyle=$
$\displaystyle\langle\psi_{4}|H^{\Gamma}_{1}(-i\partial_{x},q_{y})|\psi_{4}\rangle=0,$
$\displaystyle\langle\psi_{3}|H^{\Gamma}_{1}(-i\partial_{x},q_{y})|\psi_{4}\rangle$
$\displaystyle=$ $\displaystyle\mathcal{N}_{1}(A-i\Delta_{p})q_{y},$
$\displaystyle\langle\psi_{4}|H^{\Gamma}_{1}(-i\partial_{x},q_{y})|\psi_{3}\rangle$
$\displaystyle=$ $\displaystyle\mathcal{N}_{1}(A+i\Delta_{p})q_{y},$ (S55)

with the factor
$\mathcal{N}_{1}=\frac{\Delta_{p}}{A}(1-\frac{A^{2}}{4BM_{0}})$. Here, we used
$\tau_{x}|\tau_{y}=\pm 1\rangle=\pm i|\tau_{y}=\mp 1\rangle$ and
$\tau_{z}|\tau_{y}=\pm 1\rangle=|\tau_{y}=\mp 1\rangle$. Therefore, the
effective Hamiltonian of the helical edge states reads

$\displaystyle\tilde{H}^{\Gamma}_{1}(q_{y})=\mathcal{N}_{1}(Aq_{y}s_{x}-\Delta_{p}q_{y}s_{y}),$
(S56)

with $\bm{s}$ being the Pauli matrices in the $\psi_{3,4}(x)$ space. The bulk
topological invariant can be directly calculated using Ref. [58]. For the used
parameter setting $t=1,t^{\prime}=-0.5,\mu=-1.6$, there are four Fermi
surfaces: the pockets around $\Gamma$ and $X/Y$ points are attributed to the
positive helicity (+) band and the pocket around the M point is attributed to
the negative helicity (-) band. As the pairing potential for the pockets
around $\Gamma$ and M changes sign, the system is $\mathbb{Z}_{2}$
topologically nontrivial according to
$N_{2D}=\Pi_{s}[\text{sgn}(\Delta_{s})]^{m_{s}}$ [58] and the Dirac cone in
the edge states should be located at the projection of the M point, namely
$\overline{X}/\overline{Y}$ in the 1D edge Brillouin zone.

Next, we consider the spin-singlet pairing perturbation
$H^{\Gamma}_{2}=\Delta_{S}(\bm{q})\eta_{y}\tau_{y}$ with a real gap. It can be
easily shown that this perturbation cannot gap out the edge states. Then we
consider the $iB_{1}$ pairing perturbation
$H^{\Gamma}_{3}(\bm{q})=\frac{1}{2}\Delta_{d}(q^{2}_{x}-q^{2}_{y})\eta_{x}\tau_{y}$,
which breaks time-reversal symmetry. With open boundary conditions in the $x$
direction, we replace $q_{x}$ by the momentum operator and obtain
$H^{\Gamma}_{3}(\bm{q})=-\frac{1}{2}\Delta_{d}(\partial^{2}_{x}+q^{2}_{y})\eta_{x}\tau_{y}$.
Projecting this into the basis of $\psi_{3/4}(x)$ gives

$\displaystyle\tilde{H}^{\Gamma}_{3}(q_{y})=-\frac{\mathcal{N}_{2}\Delta_{d}}{2}s_{z},$
(S57)

with $\mathcal{N}_{2}=\frac{A^{2}-4BM_{0}}{2B}$. Thus, the time-reversal
breaking $iB_{1}$ pairing introduces a mass term for the edge states.
Moreover, we find that the mass term changes sign between the $x$-terminated
and $y$-terminated edges, due to the $d$-wave pairing nature. Hence, the edge
Hamiltonians on the $x$\- and $y$-terminated edges are given by

$\displaystyle\tilde{H}^{\Gamma}_{x}(q_{y})=\mathcal{N}_{1}(Aq_{y}s_{x}-\Delta_{p}q_{y}s_{y})-\frac{\mathcal{N}_{2}\Delta_{d}}{2}s_{z},$
(S58)
$\displaystyle\tilde{H}^{\Gamma}_{y}(q_{x})=\mathcal{N}_{1}(Aq_{x}s_{x}-\Delta_{p}q_{x}s_{y})+\frac{\mathcal{N}_{2}\Delta_{d}}{2}s_{z},$
(S59)

respectively. When these two edges meet at a corner, the mass term must go
through zero. As a consequence, a Majorana zero-energy state appears at the
corner, indicating higher-order topological superconductivity.

Expanding the BdG Hamiltonian (S46) around the M point, we obtain an
expression of the same form as Eq. (S47) but with the modified parameters

$\displaystyle M_{0}$ $\displaystyle=$ $\displaystyle 4t-4t^{\prime}-\mu,$
$\displaystyle B$ $\displaystyle=$ $\displaystyle-t+2t^{\prime},$
$\displaystyle A$ $\displaystyle=$ $\displaystyle-t\lambda.$ (S60)

With our parameter choices for $t$, $t^{\prime}$ and $\lambda$, $B<0$ and
$A<0$ but $AB>0$ still holds. Therefore, the wavefunctions and the
corresponding edge Hamilitoians remain similar. Therefore, Majorana corner
modes still appear.

##  SVI Disorder effects on the superconducting states

Here we study the disorder effects of the three different kinds of topological
superconducting states obtained in the phase diagram shown in Fig. 1 of the
main text. We introduce non-magnetic chemical potential disorder by adding a
term $H_{V}=\sum_{i}V_{i}n_{i}$ to the Hamiltonian with the disordered
Hamiltonian given by $H_{\rm{dis}}=H_{\rm{eff}}+H_{V}$ where $H_{\rm{eff}}$ is
given in Eq. (S29) and $V_{i}$ is a site-dependent non-magnetic impurity
potential drawn from a random distribution, such that $V_{i}\in[-V/2,V/2]$
uniformly, also commonly known as Anderson disorder. With the new method of
fRG+MFT developed in this paper, we can now investigate the disoder effects of
the different superconducting phases. Using the same method as in Sec. SIII,
we self-consistently solve for all the order parameters with periodic boundary
conditions in a $30\times 30$ lattice. In the presence of disorder, all the
order parameters vary at every site. Hence to quantify the average
superconducting properties, we define average of the order parameter
magnitudes as $\langle\Delta_{s/d}\rangle=1/N\sum_{i}|\Delta_{s/d}|$ and
$\langle\Delta^{t}_{A_{1}/B_{1}}\rangle=1/N\sum_{i}|\Delta^{t}_{A_{1}/B_{1}}|$
with the site-dependent order parameters given by the definitions in Eq. (S32)
and Eq. (S33). In Fig. S2, we show the disorder-dependence of average order
parameter magnitudes for three representative points of the phase diagram
marked by three crosses in Fig. 1 of the main text with different dominant
superconducting instabilities. For $V=0$ the dominant instabilities are
clearly visible in (a) and (b): for $B_{1}$ pairing state it is the singlet
$\Delta_{d}$ and for $A_{1}$ pairing state it is the $\Delta^{t}_{A_{1}}$. In
(c), the coexistence of $B_{1}$ and $A_{1}$ pairing states is also apparent
with all the order parameters being comparable in magnitude. With increasing
disorder, there is a reduction in the dominant order parameters. However, the
dominant orders survive with considerable magnitudes even upto $V=2.5$ for all
the three cases making them quite robust to disorder. Since the computational
cost increases in the presence of disorder, the disorder effect on the whole
phase diagram is left for a future work. We have also verified that the edge
properties of all three topological state do not change atleast for $V\leq
2.5$. The results in this section are only presented for a single
configuration of disorder. We do not expect the disorder-dependence of the
order parameters to change after disorder averaging.

Figure S2: Dependence of average order parameter magnitudes as a function of
disorder strength $V$ for three representative points of the phase diagram,
(a) $\mu=-0.5$ and $t^{\prime}=-0.5$ (cross I in Fig. 1 of the main text), (b)
$\mu=-1.6$ and $t^{\prime}=-0.5$ (cross III in Fig. 1 of the main text), and
(c) $\mu=-1.2$ and $t^{\prime}=-0.5$ (cross II in Fig. 1 of the main text).
$V$ is expressed in units of t.