Interaction-driven first-order and higher-order topological superconductivity

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 employ 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 + i B_1$ pairing. While the $A_1$ and $B_1$ states have first order topology with helical and flat-band Majorana edge states, respectively, the $A_1 + i B_1$ pairing exhibits second-order topology with Majorana corner modes. We investigate the disorder stability of the bulk superconducting states, analyze interaction-induced instabilites of the edge states, and discuss implications for experimental systems.

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].A variety of heterostructures [7] and candidate materials [8][9][10][11], where topological superconductivity is expected to occur, have been investigated.However, despite tremendous efforts, the ideal topological superconductor suitable for the envisioned applications has yet to be found.Two major obstacles in this research field are small superconducting gaps and omnipresent disorder.Since the topology only protects against perturbations smaller than the superconducting gap, it is of paramount importance to find intrinsic topological superconductors with larger gaps compared to proximty-induced superconductors [12,13].
Finding intrinsic topological superconductivity needs superconductivity to be originating from electronelectron repulsive interactions since electron-phonon interactions commonly only give conventional nontopological superconductivity [14,15].One possible path to high-T c topological superconductivity is to search for platforms using high-T c cuprate superconductors, for example as proposed in twisted bilayer cuprates [16].However, the lack of tunability in material control of cuprates has made it hard to experimentally realize twisted bilayers of cuprates.An alternative strategy investigated in recent times is to use a large density of states at the Fermi level [17].This strategy has been envisaged in material candidates where the Fermi level lies close to a van Hove singularity [18] or near flat bands [19].But a major hindrance towards such possibility of superconductivity is that a large density of states is also associated with particle-hole orders, like ferromagnetism, due to the Stoner criterion [20] or spin/charge density waves due to nesting [21,22].However, doping away from a ferro-magnetic phase to its critical point can result in triplet topological superconductivity, due to ferromagnetic fluctuations [23].But the resulting triplet superconductivity is highly fine-tuned, exists only in a narrow region of parameter space, and has an exponentially small T c , due to the large distance from the singular density of states [24].
Another ingredient which is often believed to be crucial for topological superconductivity is spin-orbit coupling (SOC).SOC is known to split one van Hove singularity into two and also changes the topological nature of the Fermi surfaces.Hence, the presence of spinorbit coupling gives the possibility of having high density of states in a large parameter space and additionally gives non-trivial topological band structures, providing a possible route to intrinsic topological superconductivity with large T c .With recent progress in experimental techniques, it has now become possible to fabricate 2D van der Waals materials [25][26][27] with high tunability both in doping or filling and spin-orbit coupling [28][29][30][31].
Motivated by this, we theoretically analyze the conditions under which topological superconductivity emerges in the Rashba-Hubbard model [32][33][34][35][36][37][38][39] in 2D, to capture the simultaneous role of electronic correlations driven by Hubbard interaction, SOC within Rashba model, and singular density of states at van Hove singularities in 2D.Theoretically investigating interactions and singular density of states is challenging, since methods like meanfield or random phase approximation calculations cannot capture the mutual interference between particle-particle and particle-hole instabilities at singularities [40], and Quantum Monte Carlo applied to the Rashba-Hubbard model is plagued by the sign problem.We therefore employ the functional renormalization group (fRG) [41][42][43], which treats all instability channels on an equal footing, and augment it with a mean-field theory (MFT) in or-der to capture the nature of the superconductivity deep inside the phase.
Specifically, we find that while magnetism is suppressed by SOC, magnetic fluctuations near the van Hove singularities lead to A 1 and B 1 superconducting states, with dominant p-and d-wave like pairing character, respectively, in a large parameter regime.Both of these pairing states exhibit nontrivial first order topology [3] and host helical and flat-band Majorana edge states, respectively.Remarkably, we also find near the phase boundary between these two states a time-reversal breaking A 1 +iB 1 pairing state with higher-order topology and Majorana corner modes [44][45][46][47][48][49][50][51][52].Furthermore, we extend the fRG+MFT method to real-space, which allows us to deduce the topological and edge properties of the superconducting states, as well as the (in-)stabilities of the edge states against residual interactions, and the stability of these states to disorder.Our fRG+MFT approach in real space 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 [53][54][55].We further find that the Majorana corner modes in the A 1 + iB 1 pairing state remain robust to disorder, defying the usual expectation that the sign-changing B 1 pairing, responsible for the corner modes, is sensitive to disorder.Model and methods.-Westart from the Rashba-Hubbard Hamiltonian on the square lattice given by H = H 0 +H U , with where c j,σ (c † j,σ ) is the annihilation (creation) operator of an electron at site j with spin projection σ, ⃗ τ are the Pauli matrices, r j is the lattice coordinate of site j, n j,σ = c † j,σ c j,σ is the spin-resolved particle number operator, µ is the chemical potential.In the following, we choose the hopping amplitudes such that t ⟨j,j ′ ⟩ = −t when j and j ′ are nearest neighbors, t ⟨⟨j,j ′ ⟩⟩ = −t ′ when j and j ′ are second neighbors, and zero otherwise.We set t as the energy unit.We also consider the Rashba SOC to be nonzero only for nearest neighbors λ ⟨j,j ′ ⟩ = λt.In the rest of the paper, we take λ = 0.3.Results for different values of λ are shown in the Supplemental Material (SM) [56].
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) [57][58][59][60].This results in a renormalized su- perconducting gap equation of the form where ν n are fermionic Matsubara frequencies, T is the temperature, ∆ σσ ′ k is the SC order parameter, and F ss ′ (k ′ , ν n ) is the spin-resolved anomalous propagator, describing the propagation of a hole that gets reflected into a particle or vice-versa.The function V σσ ′ ss ′ kk ′ describes an effective interaction, computed by means of the fRG (see SM for details [56]).To determine the symmetry (and not the size) of the SC gap, we linearize Eq. (2) in ∆ σσ ′ 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.-InFig. 1, we show the superconducting phase diagram as a function of chemical potential µ ≤ 0 and second-neighbor hopping t ′ ≤ 0. We also show the different topologies of the two Fermi surfaces, 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 marked differently: VHS1 (left line) has a larger density of states than VHS2 (right line) and is also located further away from the (π, 0) point, see SM [56].In the orange region of Fig. 1, labeled as iAF, the fRG flow signals an (incommensurate) antiferromagnetic instability.In the region enclosed by the dashed-dotted black lines the superconducting transition temperature exceeds T = 10 −6 t, which is the lowest temperature accessible in our fRG computation.
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) [56].The resulting gap function can to leading order be expressed as where t µ = iτ µ τ 2 , and α s and α 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 SOC.Higher-order harmonics are considered in our calculations, but the effective attractions in the higher-harmonic channels are found to be very small.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 swave in the singlet component and a helical p-wave in the triplet: For larger −µ and small −t ′ , 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 a dominant singlet component in the d xy symmetry channel, and a subdominant triplet part.We note that near the boundary between differently colored regions in Fig. 1 there is the possibility of a coexistence phase that combines the two order parameters (see discussion below).Hence, in general there will be two transitions as one goes from, e.g., deep inside the B 1 state to the A 1 state.Finally, we note that SOC disfavors ferromagnetic phases expected near van Hove singularities [24].Instead, spin fluctuations drive the formation of the SC phases in a large parameter regime even at van Hove fillings, see SM [56].
Having identified different superconducting phases in Fig. 1, we then focus on each of them individually by investigating the characteristic points marked by red crosses in Fig. 1.Additionally, we extend the fRG+MFT method to real-space to investigate the disorder stability of the phases and to determine possible interactioninduced edge instabilities [56].Furthermore, we scale the interaction strengths by 10 in order to render the investigation of the edge properties on finite size systems computationally amenable.B 1 pairing state.-InFig. 2, we investigate the properties of the B 1 phase for the parameters marked by cross I (µ = −0.5 and t ′ = −0.5) in Fig. 1.In this regime, the gap functions show (quasi-) nodes along the diagonal k x = ±k y , see SM [56].This nodal structure suggests that a (11) edge is pair-breaking [61,62] and will form a flat band of topological zero-energy states [63][64][65][66]  protected by time-reversal symmetry and translational symmetry along (11).The large degeneracy of these zero-energy states makes them thermodynamically unstable and prone to symmetry breaking.To investigate the edge properties, we consider an open boundary geometry (see Fig. 2(a)) with one (11) edge.As shown in (a), at low temperatures, the (11) edge hosts oscillations in the phase θ of the d-wave superconducting order parameter, called phase crystals breaking both timereversal symmetry and translational symmetry along the (11) edge [53][54][55].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 triplet superconducting order parameter originating in the bulk due to SOC.We calculate the spatially averaged density of states where N is the total number of lattice sites, and u n iσ and v n iσ 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. 2(b), N (E) shows a large zero-bias peak for T /T c > 0.17, showing the presence of a flat band of zero-energy states, which does not change with increasing temperature.Due to the formation of the phase crystal at T /T c ≈ 0.17, the zerobias peak gets suppressed for lower temperatures since the phase crystals Doppler shifts the zero-energy states to finite energies.With lowering temperature, the shift increases.A 1 pairing state.-Wenow investigate the topological edge states of the A 1 phase for the parameters marked by cross III (µ = −1.6 and t ′ = −0.5) in Fig. 1.This A 1 pairing is a time-reversal symmetric fully gapped superconductor and belongs to the DIII class, characterized by a 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 = Π s [sgn(∆ s )] ms , where m s is the number of timereversal invariant points enclosed by the sth Fermi sur- face [67].Due to this non-trivial topology, we find helical edge states with open boundaries, shown in Sec.SVII of the SM [56].A 1 + iB 1 pairing state.-Thetransition from the B 1 superconducting state to the A 1 state with increasing |µ| 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 (µ = −1.2 and t ′ = −0.5) in Fig. 1, which is near the phase boundary.Interestingly, we find a timereversal symmetry breaking A 1 + iB 1 superconducting phase as the lowest energy state.To verify that the relative phase of π/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 ϕ.As shown in the inset of Fig. 3(a), ϕ = ±π/2 gives the largest condensation energy and consequently A 1 ± 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 secondorder topological superconductor [44], see SM [56].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. 3(a).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 ϱ s , where n ′ are the four zero-energy states.In Fig. 3(b), we see that ϱ 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.
To analyze the stability of the A 1 +iB 1 state to perturbations, we study the disorder effects on this phase.We introduce non-magnetic chemical potential disorder by adding a term H V = i V i n i to the Hamiltonian, with V i being a non-magnetic impurity potential drawn from a random distribution, such that V i ∈ [−V /2, V /2] uniformly (i.e., Anderson disorder), and perform a fully selfconsistent calculation.As shown in the inset of Fig. 3(c), the bulk gap E g reduces with increasing disorder strength V , but remains finite for realistic disorder strengths of V ≤ 2.0.Notably, the average order parameters, shown in Fig. S4(c) of the SM [56], show more robust behavior.The disorder-robust behavior of A 1 + iB 1 is remarkable since the broken time-reversal symmetry makes Anderson's theorem [68] not applicable, and might be related to the presence of interactions [69][70][71][72].We also look at the corner states.Fig. 3(c) and (d) show the eigenvalues and the spin projected wavefunction of the lowest energy states for V = 1.5.The Majorana corner states persist even in the presence of disorder, showing the stability of the higher-order topological superconductor.This finding is also striking from the following point of view.It is known that topological nature of first-order topological phases make Majoranas survive moderately strong disorder [73][74][75].However, the corner modes in the higherorder A 1 + iB 1 arise due to the sign change of B 1 pairing at adjacent edges.Now, it is commonly believed that B 1 pairing is sensitive even to non-magnetic disorder [76].Hence, the persistence of Majorana corner modes in the presence of disorder is highly non-intuitive and opens a new way of understanding sensitivity of higher-order topological superconductors.It will also be interesting to investigate in the future the effects of other models of disorder.
Discussion.-We have shown that the combined effects of van Hove singularities, Rashba SOC, and repulsive Hubbard interactions give rise to a cascade of topological su-perconducting 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, which is also disorder-robust.
It is remarkable that we find interaction-driven topological superconducting phases in the vicinity of van Hove singularities due to the presence of SOC.We emphasize that the topological superconducting phases emerge robustly in the Rashba-Hubbard model, independent of the size of the Rashba SOC and survive the inclusion of further neighbor hoppings.The size of the Rashba SOC however decides the doping window where topological superconductivity is obtained.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 spinlattice symmetry groups.While we have only considered effects of a Hubbard onsite interaction, longer ranged interactions can bring in additional interesting aspects, which will be investigated in a future work.For example, on kagome lattice systems, nearest-neighbor interactions play a crucial role due to the sublattice interference effects [77,78].
Our results provide a guide to understand and design topological superconductivity in heavy-atom superlattices [79] and van der Waals materials [25][26][27].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 [80] or EuO/KTO [81] it is possible to tune the carrier density, and therefore the Fermi surface, by electric gating.With the recent progress in experimental techniques, the strength of the SOC can also be highly tuned.It can either be tuned by applying an electric field [28][29][30][31] or by changing the geometry of a superlattice.For example, SOC in the CeCoIn 5 /YbCoIn 5 superlattices can be adjusted by the width of the YbCoIn 5 blocks [28,82].In CeCoIn 5 /YbCoIn 5 superlattices [28,82], there are already prospective signatures of topological superconductivity below T c ≃ 2 K, a high value compared to topological superconductivity proposed in semiconductorsuperconductor nanowire devices [12,13].Other promising candidates can be van der Waals heavy-atom materials [83][84][85].All these experimental developements and our findings taken together may open up a route to intrinsic topological superconductors being used for the design of quantum information devices.Moreover, the A 1 + iB 1 superconductor may show interesting diode [86,87] and piezoelectric effects [88], since it breaks both inversion and time-reversal symmetry.
We are thankful to A. Greco  In order to determine the system's instabilities and the leading superconducting symmetry, we run a truncated functional renormalization group flow [1][2][3][4].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 [5].By rescaling the Grassmann fields according to where T is the temperature, we obtain the action where n labels fermionic Matsubara frequencies and Ĥk is the Fourier transform of H 0 in Eq. ( 1), reading as The temperature therefore takes the role of a renormalization group scale and one can derive an exact evolution equation for the effective action functional [6,7] Here, we have defined and where x = (k, n, σ) and the symbol [•] T denotes matrix transposition and does not have to be confused with the temperature T .Expanding the effective action functional Γ T [η, η] on both sides of (S4), one can derive an infinite hierarchy of flow equations for the n particle vertices Eq. ( S4) is completed with an initial condition for the effective action reading as [7] In this work, we truncate the effective action at the two-particle level, that is, we set to zero all the three-and higher particle correlators.We also ignore the flow of the self-energy, that is, we keep the one-particle vertex fixed to ( ĜT 0 ) −1 .The only vertex function we retain is the two particle vertex 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 [1,2]: 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. (S8) that at T → +∞ one has V = U and therefore P = D = 0.The flow equations read as with The symbol Π is a shorthand for dΠ/dT .To simplify the treatment, we parametrize the spin dependence of P and D as where τ 1,2,3 are the Pauli matrices, τ 0 = 1, and t α = iτ α τ 2 .To treat the dependence of the channels on k and k ′ , we resort to a truncated unity approach [4,8], that is, we expand them in a complete basis of form factors {f ℓ k }: 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 ℓ 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 [9].It is therefore convenient to truncate the sum (S14) 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 Inserting the identity into (S11), we obtain Similar expressions hold for the bubbles Πβα X,ℓℓ ′ (q).For more details on how to perform the vertex projections, see for example Refs.[4,8].
We run a flow starting from a high-temperature value T ini ≃ 10t at which P and D can be initialized by a second order perturbation theory result: to (approximately) account for the integration from T = +∞ to T = T 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 min = 10 −2 t where the (adaptive) momentum integration in the projections of the bubbles (S12) onto form factors becomes numerically unstable.

B. Renormalized mean-field theory
To continue the flow below the stopping temperature T s ≡ max(T c , T min ), and therefore detect the leading superconducting instabilities, we resort to a combination of mean-field theory with the functional renormalization group [10][11][12].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 ≤ T ≤ T s ) where F σσ ′ (k, ν n ) are the spin components of the anomalous propagator F (k, ν n ), defined through with ∆k the gap matrix (in spin space).The function V σ1σ2σ3σ4 kk ′ 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 : Eq. (S22) can be inverted with an approach similar to the one presented in the previous section: one can paramterize the spin dependence of V , V Ts , and Π Ts P by means of the t α matrices (see Eq. (S13a)) and expand the k-and k ′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 Ts .One obtains Defining a multi-index i = (α, ℓ), V , P [V Ts ], and Π Ts P become matrices, and V can be cast as By expanding also the gap in form factors and t α matrices, the gap equation takes the form where we have defined and Notice that Eq. (S25) still depends on the temperature via F α ℓ .Eq. (S25) can be used to determine the gap function at a given temperature T ≥ T s but also to simply determine the leading superconducting state.This can be achieved by linearizing it with respect to ∆ α ℓ .By noticing that we obtain For temperatures close to the superconducting transition temperature T 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 V • Π P (q = 0) with eigenvalue 1. T SC can therefore be determined as the temperature at which the largest positive eigenvalue of the matrix V • Π P (q = 0) becomes 1.Since in a Fermi liquid some components of Π P (q = 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 SC , or, in other words, Eq. (S25) has always a nontrivial solution at T = 0.The numerical values of ∆ α ℓ , however, can be exponentially small, so that their calculation becomes numerically involved.For this reason, we analyze the linearized gap equation, Eq. (S29), at a low temperature, T low = 10 −6 t, and determine the leading superconducting eigenvalue λ SC , defined as the largest positive eigenvalue of V • Π P (q = 0).The symmetry of the corresponding eigenvector determines the symmetry of the leading superconducting instability.

C. Approximations
In this Section, we briefly discuss the robustness of our fRG+MFT method against the approximations we have made.First of all, we have employed a so-called one loop approximation, that is, we have neglected all three particles and higher order correlators.This approximation is reasonable in the limit of weak-to-intermediate interactions, where the fRG is known to give reasonable answers on the leading order tendencies [7,13].This is why we have chosen the intermediate value of U = 3t for the Hubbard interaction in the main text.
Another point is the neglect of the electron self-energy and of frequency dependencies in the vertex function.Such an approximation has been previously shown to give reasonable results at weak or intermediate coupling when compared with a full treatment of the frequencies and/or self-energies [11,14,15].Therefore, we do not expect this approximation to significantly impact the results presented in the main text.
The reader might wonder whether the obtained values for T c can be considered as good estimates for the transition temperatures in the Rashba-Hubbard model.Unfortunately, while the above mentioned approximations do not impact the results qualitatively, they do affect them quantitatively.Another important aspect to be kept in mind is that the one-loop truncation predicts a mean-field-like superconducting transition, while in two spatial dimensions this should belong to the Berezinskii-Kosterlitz-Thouless (BKT) universality class.Since BKT transition temperatures are typically smaller than their mean-field analogues, we can interpret our obtained values for T c as overestimates of the actual transition temperatures our model displays.However, we stress that, within the sets of approximations employed, the linearization of the gap equation does not represent a further approximation for the calculation of T c .It may only slightly affect the relative strengths of the different superconducting order parameters deep in the phase, but since the values of the transition temperatures are generally very low, we expect this difference to be minimal.

SII. SYMMETRY CLASSIFICATION OF THE SUPERCONDUCTING STATES
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 ∼ (− sin k y τ 1 + sin k x τ 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 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 α matrices, according to Eq. (S26).We notice that the action of the group on the matrices t α is given by with the different transformation matrices Û listed in Table SI.As a consequence, when α = 1, 2, 3, the matrices t α 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 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 k, in the second one how they transform the Pauli matrices ⃗ τ , 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 ⃗ τ → Û † ⃗ τ Û .By K we denote the complex conjugation operator.
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 are irreps of the sole C4v, while g E,x k and g E,y k are a basis of the twodimensional irrep E of C4v chosen such that they transform under C4v as shown in the last line of the table.The lowest harmonic contributions for these functions are listed in Table SIII.
the gap function ∆ α ℓ , corresponding to α = 0 (α = 1, 2, 3), must be symmetric (antisymmetric) for k → −k.For this reason, only the pairs of indices (α, ℓ) with α = 0 and f ℓ k even or α = 1, 2, 3 and f ℓ k odd are allowed in ∆ α ℓ .In Table SII, we list the irreps of the combined lattice-spin discrete symmetry group, together with the typical form that the gap function ∆k takes in each one of them.In Table SIII we show the lowest harmonics for the functional dependence of ∆k on 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.
k shown in Table SII.Note that we have employed the normalization convention k |g X k | 2 = 1.

SIII. REAL SPACE RENORMALIZED MEAN-FIELD
In Sec.SI B, 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, where µ denotes the nearest neighbor bonds ±x, ±y of the lattice site i, H 0 is defined in Eq. ( 1) and ∆ σσ ′ µ is the superconducting order parameter given by where the interaction strengths are with f ℓ (µ) being the Fourier transform of f ℓ k and V αβ ℓℓ ′ are obtained from the fRG calculation.H MF is diagnolized using the Bogoliubov-de Gennes transformations [16], ), where γ † nσ (γ nσ ) are the creation (annihilation) operators of the Bogoliubov quasiparticle at a state n, and u n iσ and v n iσ the eigenfunctions with eigenvalues E. The resulting eigen-system is then solved self-consistently for the superconducting order parameters ∆ σσ ′ µ .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. (S31).This helps to reduce the numer of independent self-consistent order parameters to 16 per lattice site.We then project local ∆ σσ ′ µ 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 where ∆ d is the d-wave order parameter, which in momentum space is proportional to (cos k x − cos k y ).∆ s is the extended s-wave order parameter, which in momentum space is proportional to (cos k x + cos k y ).The perfactor ϵ where ∆ σσ ′ px/py 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 Other mixed-spin triplet components are negligible compared to the equal spin triplets.
The interaction strengths V σσ ′ ss ′ µν 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. (S32) 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 h where k c is given by tan k c = − λt t+2t ′ .The dispersion can be further simplified to with The low-energy Hamiltonian around the VHS point (π, where k ′ c is given by tan k ′ c = λt t−2t ′ .The dispersion can be further simplified to 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 ′ , 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 ′ | (negative t ′ ), the VHS1 on the Γ − 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 √ 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 ′ | > 0.3) in our fRG calculations.

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 Ψ k = (c † k↑ , c † k↓ , c −k↑ , c k↓ ) the BdG Hamiltonian can be written as where the pairing term is ∆ , with ∆ S (∆ T ) being the spin-singlet (spin-triplet) pairing component.For the helical p-wave (A 1 pairing state), the pairing matrix is (S45) The BdG Hamiltonian can be written as, ) where ⃗ η and ⃗ τ are Pauli matrices in the Nambu and spin space, respectively.This Hamiltonian has time-reversal symmetry, particle-hole symmetry and chiral symmetry, satisfying +1 and S 2 = 1, the system belongs to the DIII class of topological superconductors [17] and is characterized by a Z 2 invariant in two dimensions [18].
To study the edge states, we expand the above model around the Γ point.(The expansion around the M point is briefly discussed at the end.)The effective model around the Γ point reads, H Γ BdG (q) = [M 0 + B(q 2 x + q 2 y )]η z τ 0 − Aq y η 0 τ x + Aq x η z τ y − ∆ p (−q y η x τ z + q x η y τ z ) − ∆ S (q)η y τ y .(S47) 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 HΓ x (q y ) = N 1 (Aq y s x − ∆ p q y s y ) − N 2 ∆ d 2 s z , (S58) 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 ) With our parameter choices for t, t ′ and λ, 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. MAGNETIC FLUCTUATIONS AND SINGLET TO TRIPLET RATIO
In the upper left panel of Fig. S2, we plot the parameter η, quantifying the strength of spin-triplet in the leading superconducting state.Values of η that are close to zero imply a singlet-dominated SC state, while η ∼ 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 ′ .This property can be understood by analyzing the spin fluctuations in the model, which are responsible for pairing [20,21].In the insets of panels (c) and (d) of Fig. S2, 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.We have defined them as where S α q,ω is the Fourier transform of the α-component of the spin operator.Note that if in Eq. (S62a) we had chosen to define χ trans from the y-component of the spin, the insets of panel (c) of Fig. S2 would be rotated by 90 degrees, according to the symmetries of the Rashba term.
In the main panels, we show the ratio between the susceptibility evaluated at q = 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 q = 0, viceversa, when it is close to zero the peak is far away from q = 0.This is the case of the region where a magnetic instability is found, due to a divergence of the transverse susceptibility for q near the antiferromagnetic wave-vector (π, π).The comparison between panels (a), (c) and (d) of Fig. S2 is very instructive.Indeed, we note that spin fluctuations peaked around (π, π) tend to favor a B 1 singlet-dominated pairing, while ferromagnetic fluctuations enhance a triplet-dominated A 1 state.
In panel (b) of Fig. S2, we also report the value of the leading superconducting eigenvalue λ SC , which provides information on the relative size of the superconducting transition temperature.Note that λ 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 ′ = −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 λ 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. S2) 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 ′ |, 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.

SVII. HELICAL EDGE-STATES OF THE A1 SUPERCONDUCTING STATE
To study the topological edge states of the A 1 we take the parameters marked by cross III (µ = −1.6 and t ′ = −0.5) in Fig. 1 of the main text and consider a nanoribbon geometry with open boundary conditions along x and periodic boundary conditions along y.Here, we obtain the self-consistent order parameters with periodic boundary conditions along both x and y.Using this self-consistent solution we employ the nanoribbon geometry with no further selfconsistency, in order to reach a large system size.We have, however, verified for a smaller system size that all the edge features are insensitive to self-consistency.In the inset of Fig. S3(a), 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 for the lowest positive eigenstate with E ≈ 0 in the main panel of Fig. S3(a).The wavefunction is localized on the open (10) edges with oscillating weights going away from the edge.The localization length is less than half of the length along the 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 spin density of states A s (k y , E, i x ) = n,σ sgn(σ)|u n ix,ky,σ | 2 δ(E − E n ), where u n ix,ky,σ is the Fourier transform of u n i,σ 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. S3(b), we show A s (k y , E, i x ) averaged over i x = 0, 10 from the left (10) edge.We find a linearly dispersing Dirac edge spectrum inside the gap at k y = π with the helical nature reflected by positive values for positive velocities and negative for negative velocities.

SVIII. 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 = i V i n i to the Hamiltonian with the disordered Hamiltonian given by H dis = H eff + H V where H eff is given in Eq. (S31) and V i is a site-dependent non-magnetic impurity potential drawn from a random distribution, such that V i ∈ [−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 × 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 | with the site-dependent order parameters given by the definitions in Eq. (S34) and Eq.(S35).In Fig. S4, 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 ∆ d and for A 1 pairing state it is the ∆ t A1 .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 ≤ 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.

SIX. PHASE DIAGRAMS FOR DIFFERENT CHOICES OF λ AND U
In this section, we present phase diagrams obtained for different choices of the spin orbit coupling λ and Hubbard interaction U .All plots are scans in the chemical potential µ at a fixed t ′ = −0.35t.
By inspecting Figs.S5, S6, and S7, we notice that increasing λ tends to suppress the B 1 SC channel, while it favors the A 1 one.
Comparing Figs.S6 and S8, one can convince himself that increasing U from 3t to 4t does not change the leading SC symmetry channel, but it only enhances their strength.

SX. GAP FUNCTIONS
In this section we analyze the momentum dependence of the solutions of the gap equation (S20), where, as done for the calculations in real space presented in the main text, we have rescaled the effective interactions V by a factor of 10.A clearer picture is provided by transforming these eigenvectors from the spin to the band basis  where η, η ′ = ± and U k is the transformation matrix that diagonalizes the quadratic part of Hamiltonian (1), with e iϕ k = (− sin k y + i sin k x )/ sin 2 k y + sin 2 k x .In the general case, interband pairing terms will appear, that is, In Fig. S9, we show the momentum dependence of the gap functions calculated in three points of the phase diagram where the B 1 , A 1 , and A 1 + iB 1 phases emerge.From panel (a) of Fig. S9, we see that in the B 1 phase, the triplet component of the gap function removes nodal lines along the diagonals of the singlet component.Some nodal points persists, however, in the vicinity of the Fermi surfaces, making the overall gap in the Bogoliubov bands rather small.In the A 1 phase, instead, the intraband gaps exhibit nodal points only at (0, 0), (0, π), (π, 0), and (π, π), making the quasi-particle bands fully gapped, unless one of the Fermi surfaces passes through one of the above mentioned points.Note that in the A 1 state the interband gaps vanish.Finally, in the A 1 + iB 1 phase, we observe that the intraband gap functions are large on the normal Fermi surfaces, rendering this state fully gapped.

Figure 1 .
Figure 1.Phase diagram of the Rashba-Hubbard model on the square lattice as a function of chemical potential µ and second-neighbor hopping t ′ for U = 3t.The dashed lines indicate the values of the chemical potential µvH 1 (t ′ ), µvH 2 (t ′ ) where the van Hove singularities VHS1 and VHS2 occur.The dashed-dotted lines enclose the region where Tc exceeds 10 −6 t.

Figure 2 .
Figure 2. (a) Phase crystal forming at the (11) edge.The sine of the phase θ of the d-wave superconducting order parameter is plotted in real space at T = 0 with the bulk value of sin θ being subtracted.At the (11) edge there are clear modulations visible in sin θ.(b) Density of states N (E) around zero energy for different temperatures T /Tc.The inset shows the density of states at T /Tc = 0.23 in a wider range.

Figure 3 .
Figure 3. (a) Eigenvalues with open boundary conditions along both x and y, showing four Majorana states.Upper inset: Normalized condenstation energy ∆F as a function of the relative phase ϕ between A1 and B1 order parameters.Lower inset: zoomed view of the four Majorana states.(b) Spin projected wave function of the four Majorana states.(c,d) Same as (a,b) but with disorder V = 1.5.Upper inset of (c): Bulk energy gap Eg as a function of disorder strength V .Here, T = 0.

2 iν n − Ĥk − 1 .
0) and ⃗ τ are the Pauli matrices.We note that upon rescaling (S1), only the quadratic part of the action remains temperature-dependent, with bare propagator ĜT 0 (k, ν n ) = T 1 (S3) * These two authors contributed equally ) with X = P or D, ζ P = −1 and ζ D = +2.The symbols P [V ] and D[V ] represent the projection of the full vertex (S10) onto the P and D channels and they read as µ = ±x and ϵ d/s µ = −1/1 for µ = ±y.The phase of the d-wave order parameter is characterized with θ with the definition ∆ d (i) = |∆ d (i)| exp(iθ).The triplet order parameters are mainly p-wave with the onsite equivalent given as

Figure S1 .
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 ′ .

Figure S2 .
Figure S2.Triplet-to-singlet ratio and ferromagnetic fluctuations.Panel (a): triplet strength η of the superconducting state as a function of chemical potential µ and second-neighbor hopping t ′ .Panel (b): leading superconducting eigenvalue λ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.

Figure
Figure S3.(a) Wave-function amplitude of the lowest positive eigenstate near zero-energy for a nanoribbon geometry.The inset shows the eigenvalues.(b) Momentum resolved spin density of states averaged over a width near the left (10) edge of the nanoribbon, ⟨As(ky, E, ix)⟩i x , showing the helical nature of the edge states.Here, T = 0. Maximum values of the colorbars are cut for better visibility.

Figure S4 .
Figure S4.Dependence of average order parameter magnitudes as a function of disorder strength V for three representative points of the phase diagram, (a) µ = −0.5 and t ′ = −0.5 (cross I in Fig. 1 of the main text), (b) µ = −1.6 and t ′ = −0.5 (cross III in Fig. 1 of the main text), and (c) µ = −1.2 and t ′ = −0.5 (cross II in Fig. 1 of the main text).V is expressed in units of t.

Figure S5 .
Figure S5.Superconducting eigenvalues in each of the symmetry channels as a function of the chemical potential µ for fixed t ′ = −0.35t,λ = 0.2, U = 3t.

Figure S6 .
Figure S6.Superconducting eigenvalues in each of the symmetry channels as a function of the chemical potential µ for fixed t ′ = −0.35t,λ = 0.3, U = 3t.This plots corresponds to a cut along the t ′ = −0.35tline of Fig. 1.

Figure S7 .
Figure S7.Superconducting eigenvalues in each of the symmetry channels as a function of the chemical potential µ for fixed t ′ = −0.35t,λ = 0.4, U = 3t.

Figure S8 .
Figure S8.Superconducting eigenvalues in each of the symmetry channels as a function of the chemical potential µ for fixed t ′ = −0.35t,λ = 0.3, U = 4t.
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.D.C. acknowledges financial support from Kungl.Vetenskapsakademien and C.F. Liljewalchs stipendiestiftelse Foundation.