Higher-order topological superconductivity from repulsive interactions in kagome and honeycomb systems

We discuss a pairing mechanism in interacting two-dimensional multipartite lattices that intrinsically leads to a second order topological superconducting state with a spatially modulated gap. When the chemical potential is close to Dirac points, oppositely moving electrons on the Fermi surface undergo an interference phenomenon in which the Berry phase converts a repulsive electron-electron interaction into an effective attraction. The topology of the superconducting phase manifests as gapped edge modes in the quasiparticle spectrum and Majorana Kramers pairs at the corners. We present symmetry arguments which constrain the possible form of the electron-electron interactions in these systems and classify the possible superconducting phases which result. Exact diagonalization of the Bogoliubov-de Gennes Hamiltonian confirms the existence of gapped edge states and Majorana corner states, which strongly depend on the spatial structure of the gap. Possible applications to vanadium-based superconducting kagome metals AV$_3$Sb$_3$ (A=K,Rb,Cs) are discussed.


I. INTRODUCTION
Since the discovery of high-temperature superconductors, investigations of superconducting instabilities driven purely by the repulsive Coulomb interaction have led to the discovery of numerous unconventional superconducting phases [1][2][3][4][5][6][7][8]. Among these, topological superconductors -hosting robust, zero-energy edge modes -have taken a key role in the pursuit of platforms for quantum computing [9][10][11]. Recently, a new classification of topological superconductors has emerged which relies on crystalline symmetries in addition to non-spatial ones [12][13][14][15][16][17][18][19][20][21][22][23]. These "higher-order" topological superconductors host non-trivial boundary states with multiple dimensionalities: in two dimensions, the second-order topological superconducting phase hosts a gapped spectrum of propagating modes along the edge as well as zero-energy modes localized at the corners of the system in addition to a gapped bulk quasiparticle spectrum.
In the search for candidate materials hosting these phases, it is crucial to investigate the possible mechanisms that might give rise to them in physical systems, based on a microscopic analysis of the electron-electron interactions. In this work we describe an approach to interactions in two dimensional materials with Dirac points that naturally leads to a second order topological p + iτ p superconducting phase, arising past a critical doping for arbitrarily weak values of the repulsive Coulomb interaction. Our analysis is performed in the weak coupling regime, in which Stoner instabilities and correlated insulating states are absent, but the effects we describe dominate in materials where the atomic orbitals are strongly localized, which distinguishes our theory from those previously used to study graphene.
The Coulomb repulsion is converted into an attraction via destructive interference for certain scattering channels near the Dirac points which originate from the Berry phase, causing the formation of Cooper pairs to lower the correlation energy. Thus, our theory relies on Dirac Fermi surfaces, a scenario distinct from tuning the Fermi level to a nesting density or van Hove singularity, commonly studied as an explanation for superconductivity [24][25][26][27][28][29][30]. In contrast to previously studied mechanisms for higher order topological superconductivity, we also do not require proximitization to an existing superconductor [31][32][33][34][35][36][37][38], negative-U Hubbard interactions [39], or bosonic fluctuations [40].
This pairing mechanism was previously studied in the context of an artificial honeycomb lattice [41]; here we classify the resulting superconducting states and demonstrate their topological properties. Our results follow from general symmetry considerations, but we will present specific results for kagome and honeycomb lattices. In Section II we discuss Hubbard-like models of these systems, presenting the most general set of interactions consistent with the lattice symmetries and enumerating the possible superconducting states which result. In Section III we discuss how the superconducting instabilities are affected by screening effects. In Section IV we perform exact diagonalization of the Bogoliubov-de Gennes Hamiltonian describing the superconducting states, confirming the existence of edge and corner modes for samples whose edge respects the point group symmetry of the lattice. We conclude with a discussion of candidate materials, and suggest the mechanism can explain recently observed superconductivity in vanadium-based kagome metals AV 3 Sb 3 (A=K,Rb,Cs), where superconductivity appears to be unconventional but correlations are weak [42][43][44][45][46][47][48].

II. SUPERCONDUCTING INSTABILITIES IN DIRAC MATERIALS
We investigate the pairing instabilities of a twodimensional lattice with hexagonal symmetry hosting Dirac points at the K points, which we describe via an extended Hubbard model where c † r = (c † r,↑ , c † r,↓ ) are two component spinors and the spin inner product is implied in the single particle and interacting terms in the Hamiltonian (1).
We choose lattice vectors a 1 = (a, 0), a 2 = ( a 2 , √ 3a 2 ), and specialize to lattices possessing C 6v symmetry, which consists of twofold and threefold rotations in the plane and mirror reflections about the x and y axes. We only consider models without spin-orbit interaction, thus we may choose a basis of Wannier orbitals so that T (r − r ) and U (r 1 , r 2 , r 3 , r 4 ) are real. The Hamiltonian is invariant under SU(2) spin rotations and possesses a time reversal symmetry satisfying T 2 = 1, which is identical to complex conjugation in the coordinate representation. For the remainder of this paper we denote the time reversal operation which leaves spin unaffected by T .
The resulting band structure features Dirac points at the K and K points, which we distinguish by a valley index τ = ±. The eigenstates of the single particle Hamiltonian may be classified by their eigenvalues e 2πi 3 α under threefold rotations, where α = ±. We shall refer to the α degree of freedom as pseudospin. At the Dirac points there are four degenerate energy eigenstates |τ, α which transform into one another under time reversal, twofold rotations and mirror symmetry. We shall introduce the Pauli operators τ µ , α µ acting on the valley and pseudospin degrees of freedom. The degeneracy of the α eigenstates at each valley is protected by either M x or R π T . We obtain the following representation of crystal symmetries: where φ I is a phase which generally depends on the lattice and the chosen center of inversion. Doping slightly above the valleys ±K = ±( 4π 3a , 0), the Fermi surface consists of two circular pockets surrounding the K points with Fermi momenta k F |K|, with the single particle energy eigenstates formed from linear combinations of the pseudospin basis states (i.e. eigenstates of threefold rotations). By obtaining the band structure of (1) and projecting the interactions onto states near the Dirac points, we obtain a quantum field theory describing interaction processes close to the Fermi level which accounts for the τ , α and spin (s) degrees of freedom and may be formulated using an 8-component local field operator ψ † τ αs (r). The most general Hamiltonian density consistent with the symmetries (2) is of the form where we have expanded the interaction in the adjoint basis consisting of products J a = τ µ α ν , and V ab are constants which may be obtained by projecting the extended Hubbard interactions (1) onto the pseudospin and valley eigenstates. As we show in the Appendix, only operators which are even under time reversal symmetry are permitted in the interactions, which are {J 1 , J 2 , . . . , The Hamiltonian must transform as a scalar under the spatial symmetries (2) which implies that V is a diagonal matrix. There are only five independent coupling constants g 1 , g 2 , g 3 , g 4 , g 5 , which we relate to the interaction matrix V ab by In addition to a density-density interaction g 1 , the interactions involve valley-conserving (J 2 = τ z α z ) and valleymixing (J 5 = τ x , J 6 = τ y ) mass operators. The remaining six operators are conserved SU(2) valley currents which can be expressed as products of the valley operators iτ x , iτ y , τ z and the electric current j µ = (j x , j y ) = (τ z α x , τ z α y ). The relation between the coupling constants g i appearing in the field theory (3) and the extended Hubbard model parameters U (r 1 , r 2 , r 3 , r 4 ) for the kagome and honeycomb lattices are given in Table I. We denote the three sublattices of the kagome lattice {A, B, C}, and the two sublattices of the honeycomb lattice by {A, B}. We list only interactions involving nearest neighbors, and denote the interactions in which the lattice coordinates r 1 , r 2 , r 3 , r 4 exist on the σ 1 , σ 2 , σ 3 , σ 4 sublattices by U σ1σ2σ3σ4 . In the kagome lattice, each site has two nearest neighbors in each of the other sublattices, and in the honeycomb lattice, each site has three nearest neighbors in the other sublattice.
The eigenstates of the Dirac Hamiltonian in the upper band are given by f † k,τ,s Evaluating the Born amplitudes from (3) we find that the scattering vertex Γ τ1τ2τ3τ4 (θ) in the Cooper channel is given by Kagome Honeycomb g1/Ω 1 6 (2UAAAA + 8UABAB + 8UAABA + 4UACBC + UAABB) where θ is the scattering angle. BCS theory predicts that superconducting pairing occurs when the scattering amplitude between Cooper pairs is negative, i.e. for attractive scattering. For pairing of electrons within the same valley τ this corresponds to Γ τ τ τ τ < 0, while pairing of electrons in opposite valleys requires the symmetrized The critical temperature is given where λ = Γ τ τ τ τ for intravalley pairing or λ = Γ +−+− ± Γ +−−+ for intervalley pairing, and ν 0 = k F /(2πv) is the density of states per spin per valley at the Fermi level. In order to investigate the physical features of the condensate we introduce the mean field Hamiltonian where f † k is the creation operator for a Dirac fermion in the upper band and ∆(k) is the superconducting gap matrix (spin and valley indices are implied). Accounting for spin, valley and angular momentum structure, there are eight possible superconducting phases, which we list in Table  II along with the corresponding scattering amplitude.
The topological properties of the gaps are associated with their Altland-Zirnbauer class [49][50][51], which classifies mean field Hamiltonians based on time reversal, charge conjugation and chiral symmetries. For the spin singlet pairing phases, ∆ ∝ s 0 , accounting for SU(2) spin rotation symmetry we find that the charge conjugation symmetry satisfies C 2 = −1. In the case of spin triplet pairing, the spontaneous violation of SU(2) spin rotation symmetry fixes a direction d = (d x , d y , d z ), with the gap ∆ ∝ d µ s µ . As we show in the Appendix, magnetization of the condensate is energetically disfavored, which implies that d can be chosen to be purely real. The mean-field Hamiltonian is invariant under spin rotation about the d-direction. Accounting for this U(1) symmetry [88], we may decompose the mean field Hamiltonian in two spin blocks, with each separately possessing a charge conjugation symmetry satisfying C 2 = +1. We therefore conclude that the phases with spin singlet pairing are either in class CI or C, while those with spin triplet pairing are either in class BDI or D depending on whether the time-reversal symmetry of the normal state Hamiltonian survives in the condensed phase. Consulting the periodic table of topological invariants [51], we find that first order topological superconductivity is possible for the phases in which time reversal symmetry is spontaneously broken, i.e. those in class C or D, and posses a Chern number which counts the number of chiral modes propagating along the boundary and whose parity is equal to the parity of Majorana modes localized at a vortex. However, both the intravalley p-wave spin triplet, and intervalley s-wave spin triplet states are in class BDI. For the intravalley spin triplet phase, the gap is p + ip in one valley and p − ip in the other, preserving time reversal symmetry, and we refer to this as the p + iτ p phase. In these cases lowest order topology is absent while second order topology is possible [14,16,20]. These phases are indicated by the existence of Majorana corner states of definite spin protected by crystalline symmetries. In addition to the spinless time reversal symmetry operator satisfying T 2 = +1, our system possesses a spinful time reversal symmetry T = T e iπSy satisfying (T ) 2 = −1 which interchanges opposite spin blocks; thus corner modes occur in Kramers pairs. In Section IV we present exact diagonalization results which demonstrate that the intravalley p + iτ p spin triplet phase indeed possesses a second order topology.
We now apply our previous results to analyze the possibility of superconducting instabilities in the kagome and honeycomb lattices. The leading Hubbard interactions do not involve tunneling between sites, thus we only account for interactions U (r 1 , r 2 , r 3 , r 4 ) in (1) with r 1 = r 3 , r 2 = r 4 . Denoting the parameter λ in (6) by λ 1 for intervalley s-wave spin triplet pairing, λ 2 for intervalley p-wave spin triplet pairing and λ 3 for intravalley p-wave spin triplet pairing, we find for the kagome lattice where Ω is the area of the unit cell and for the honeycomb lattice.  In the limit U ABAB → 0 the scattering amplitudes for all three spin triplet instabilities vanish for both lattices, implying that the system approaches a critical point for the onset of triplet superconductivity, including for the topologically nontrivial phases. The vanishing of the scattering amplitude in the p + ip and p + iτ p channels is a consequence of the Dirac nature of the scattering states, and can be understood as follows: the components of the wavefunction in the pseudospin and valley basis have momentum-dependent phase factors which are a manifestation of the Berry phase surrounding the Dirac points. Despite our model only containing local interactions, upon projecting the Dirac wavefunction onto the interactions in (3) we find that the Berry phase gives rise to p + ip scattering amplitudes which are essential to obtain a nontrivial topology. The pseudospin and valley dependence of the operators J a allows some g i to contribute to the partial wave amplitudes λ with negative signs, as can be seen in Table II. This causes destructive interference between the contributions to the scattering amplitude with distinct pseudospin and valley structurei.e. the g i . In the scattering channels which give rise to the topological phases, the contributions from different g i sum to zero, leading to a total cancellation of the on-site Coulomb repulsion.
An additional negative contribution to the scattering amplitude which can overcome the weak nearest neighbor repulsion would result in superconductivity. One interesting scenario would be a weak electron-phonon interaction, which in our case could drive the instability to an unconventional spin triplet state, rather than the spin singlet state which normally results. In the next section, however, we will show that even in the absence of additional attractive interactions, a negative contribution to the nearest neighbor repulsion arises as a result of overscreening. All three triplet amplitudes become negative if the bare nearest neighbor repulsion is sufficiently weak, and the density is sufficiently high.

III. SCREENING EFFECTS
We now consider how g i are modified by screening processes. We show these effects here analytically in the Random Phase Approximation (RPA), equivalent to taking the limit of large N where N = 4 is the degeneracy of the Dirac points [53][54][55][56]. This analysis was performed for the special case of an artificial semiconductor honeycomb lattice in [41]; we extend this analysis to encompass generic Dirac materials satisfying C 6v and time-reversal symmetry, which may be described by the Hamiltonian (3). The RPA results in replacing the interaction constants V ab in (3) with screened frequency and momentum dependent couplings, which are determined by the RPA equation where Π cd (ω, q) is a generalized polarization operator. These represent the susceptibility of the system to a perturbation δH = a φ a (r, t)J a , and are calculated in the Appendix. Accounting for the frequency and momentum dependence of the screened scattering vertex, solution of the Eliashberg equations gives the critical temperature where λ is the RPA renormalised λ of Table II obtained from interactions between electrons on the Fermi surface and Λ is a frequency cutoff determined by scattering processes away from the Fermi surface; explicit calculation (discussed in the Appendix) shows Λ ≈ E F . Scattering on the Fermi surface corresponds to the range of momentum and frequency transfer q < 2k F , ω = 0. In this range the polarization operators are constant, and λ is simply given by the previous expressions λ (Table  II) evaluated with screened couplings g i at ω = 0, q = 0. We have studied the screening of the interactions g 1 and g 2 in a previous paper [41], and found that g 1 is reduced while g 2 is enhanced as the chemical potential is increased. We can explain the qualitative effects of static and homogeneous screening on the interactions on general physical grounds. The coupling g 1 is screened by the long wavelength density-density response, which as usual weakens g 1 with increasing chemical potential. The intra and intervalley mass couplings g 2 and g 4 are screened via the polarization operator associated with the response of the system to an external perturbation ∝ α z τ z , τ x , τ y which gap the Dirac spectrum. This perturbation lowers the energy of negative energy states and increases the energy of the positive energy states. The response of the occupied negative energy states below the Dirac point is divergent and may be fully subtracted by regularizationthe result of which is that g 2 and g 4 should be replaced by reduced interaction constants which are different from those provided in Table I; thus within the Dirac model the only physically significant effect of screening is the increase in energy of the occupied positive energy states which grows stronger with increasing chemical potential. The interactions g 3 , g 5 are valley current interactions and the polarization operators responsible for their screening also contain a UV divergent part which may be subtracted by replacing the bare interaction constants with reduced values, and after this subtraction are unscreened due to a Ward identity. Explicitly, we find in RPA where i = 2, 4 and j = 3, 5. This results in replacing the nearest-neighbor coupling with a screened value, U ABAB → U ABAB in the relations in Table I and (8), (9). The screened nearest neighbor interaction is given by for the kagome lattice and for the honeycomb lattice. For sufficiently large ν 0 , we find that U ABAB < 0, resulting in superconductivity, and the dominant instability occurs in the intravalley, p-wave, spin triplet channel. This occurs for µ > µ * with the critical value for the kagome lattice and for the honeycomb lattice. The Dirac theory is only applicable up to a cutoff beyond which the electronic dispersion becomes nonlinear. Our analysis predicts a superconducting instability when the critical doping µ * lies within the regime where Dirac theory is justified. In the limit g 3 → 0 or g 2 → 0 respectively for the kagome and honeycomb lattices, the critical doping becomes infinite. We therefore require significant pseudospin dependent couplings, which are suppressed when the atomic orbitals are delocalized.
So far, we have neglected the long ranged part of the Coulomb interaction, which provides weak couplings between sites with large separation. This may be accounted for by replacing g 1 with its screened value at q → 0, g 1 → 1/N ν 0 , and simply results in a shift of the critical doping to higher values.

IV. TOPOLOGICAL EDGE AND CORNER MODES
Our RPA results show that the dominant superconducting instability for both the honeycomb and Kagome lattice is the intravalley p-wave spin triplet phase, with . This gap possesses a number of unconventional features that distinguish it from the superconducting phases more commonly studied in either high-T c superconductors or Dirac materials. The gap exists at both K points and has a p + iτ p wave angular structure, which points to the presence of a π-Berry phase in two-particle scattering from which superconducting pairing originates. However, the winding of the gap is opposite in the two valleys, with the superconducting order parameter at the two K points being related to each other by complex conjugation. As a result, the usual topological features associated with chiral p-wave superconductors are strongly modified. We may view the phase as consisting of four copies of a topological p + ip superconductor [10], one per spin and valley species, with the gaps at opposite valleys possessing topological invariants of opposite sign. This results in a topologically trivial superconductor according to the Altland-Zirnbauer classification. However, in the presence of crystalline symmetries, a second order topology emerges which is characterized by anomalous corner modes.
In the intravalley paired state, Cooper pairs carry a finite quasimomentum, and as a result the gap is spatially modulated, forming a pair density wave which is commensurate with the lattice with a periodicity of three unit cells [57][58][59][60][61][62][63][64]. The condensate possesses a U(1) order parameter φ associated with the spatial position of the pair density wave and does not couple to the electromagnetic field; the system also exhibits Goldstone modes associated with fluctuations of φ, which are physically the sliding modes of the pair density wave. The spatial modulation of the gap spontaneously breaks the crystalline symmetries of the normal state except at special values φ = nπ 2 , n ∈ Z; at these values the gap may be classified according to a second order topological invariant which predicts the number of Majorana corner modes of a finite sample of the superconductor which preserves the crystalline symmetries.
We study the topology of the p + iτ p intravalley spin triplet phase via the bulk boundary correspondence for higher order topological superconductors, and identify the anomalous boundary physics in hexagonal superconducting flakes. By choosing a spin quantization axis which is aligned so that d ŷ, we find that pairing occurs between electrons with both spins aligned either along the z or −z axis. Thus the pairing term may be decomposed into a sum of two spin sectors. We derive the topological features of the superconducting state by analyzing the spinless Hamiltonian corresponding to pairing within a single spin sector. In order to make numerical diagonalization easier we employ a simplified form of the gap function in which we only include nearest neighbor pairing and hopping terms. In the Appendix, we derive the real space form of the gap from the p + iτ p momentum space structure.
The spinless mean field lattice Hamiltonian which results is where ∆(r , r) = −∆(r, r ), and ∆(r, r ) = ∆ cos(K · (r + r ) + φ) , for the kagome lattice, and for the honeycomb lattice, and we distinguish the parameter ∆ from the bulk gap ∆. We present exact diagonalization results for the kagome and honeycomb lattices in Figs. 1, 2 respectively. We have calculated the spectrum for all values of φ, however the topological classification must be performed at the values φ = 0, π 2 for which the gap obeys a twofold rotation symmetry. We have performed calculations for a variety of edge geometries, and in each case we preserve the twofold rotational symmetry of the bulk. In a number of cases we observe a gapless edge spectrum for all values of φ. This occurs for the geometries shown in the last two rows of Fig. 1 and the last row of Fig. 2. These surfaces states are remnants of the chiral Majorana edge modes expected from the Read-Green model [10] which are the most easily identifiable manifestions of the chiral p-wave nature of the gap. For other edge terminations, however, the edge spectrum is gapped, indicating that the gapless edge modes are not topologically protected, which is consistent with the fact that, since our gap lies in class BDI, it cannot possess a first order topological invariant. For the other edge terminations, we find that at φ = π 2 , the edge states are gapped, and six zero energy modes appear within each spin sector which are localized at each of the corners of the flake, forming the three lowest modes within the quasiparticle spectrum of a single spin sector highlighted in red in Figs. 1, 2c,f,i. This occurs for both the honeycomb and kagome lattices. At φ = 0, we either find no Majorana zero modes when the edge spectrum is gapped or a gapless edge spectrum. This allows us to conclude that both the honeycomb and kagome lattices are in a second order topological phase for φ = π 2 and trivial for φ = 0. We note that, since our results are for a spinless model obtained after performing a spin decomposition of the mean-field Hamiltonian, the Majorana zero modes appear as Kramers pairs. Since pairing occurs between electrons with parallel spin, the Majorana modes individually possess definite spin, and are protected against hybridization by a spinful time reversal symmetry T = e iπSy T .
Since the boundary modes in a higher order topological phase require crystalline symmetries, they are not protected against local perturbations. The robustness of the bulk-boundary correspondence in the presence of crystalline symmetry breaking therefore provides an important, general question for systems exhibiting higher order topology. In our examples, for general values of the pair density wave order parameter φ, all the point group symmetries of the normal state are spontaneously broken by the pairing term except for mirror reflection about the x axis, which is preserved due to the fact that the phase K · (r + r ) + φ of the pair density wave does not depend on the y coordinate. However, at discrete values φ = πn 2 , the superconductor recovers a twofold rotational symmetry. The parameter φ allows continuous tuning of the system between a second order topological phase at φ = π(n + 1 2 ) and a trivial phase at φ = nπ. We observe, despite the breaking of crystalline symmetries, that the anomalous corner modes persist and remain exponentially close to zero energy in a finite region of values around |φ − π(n + 1 2 )| < φ * up to a critical value φ * at which the edge gap closes. The corner modes are protected by both the bulk and edge gap, and thus persist even when the point group symmetry is weakly broken [13,14].
The anomalous corner modes coexist with a gapped edge spectrum. This originates from the p + iτ p nature of the gap, which provides a 1D boundary superconductor for each copy of the Moore-Read model existing in each spin and valley species. The manifestation of the bulk second order topology is shadowed by the behavior of the boundary superconductor: for certain edge terminations, the edge modes remain anomalous due to symmetries along the edge which forbid their hybridization. In particular, for the bottom geometries of Fig. 1 and 2, the system possesses a mirror symmetry at φ = πn 2 that prohibits the counterpropagating chiral Majorana modes from gapping out. In this case, the edge states are anomalous and associated with a bulk crystalline symmetry protected topological phase protected by mirror symmetry. For other edge terminations, mixing of the chiral modes between opposite valleys along the same edge occurs, and the edge spectrum is gapped. Interestingly, the emergence of anomalous corner modes in this situation is not general but depends on the value of φ.
We may make an interesting comparison to the p + ip intervalley spin triplet phase, which is subleading in our analysis in the regime where the onsite couplings dominate over the nearest neighbor couplings, but may arise as the dominant instability for models with more complicated lattice interactions. In this case, pairing occurs between the electron pockets surrounding opposite valleys, and the order parameter spontaneously chooses a winding number corresponding to either p + ip or p − ip pairing. Time reversal symmetry is broken, which puts the superconductor in class D, rather than BDI as in the intravalley case. This phase possesses a Chern number which leads to a gapless, anomalous edge spectrum. In a forthcoming study, we show that this phase may become dominant when spin orbit coupling is present, and analyze the topology of this phase in detail.

V. DISCUSSION
Our results show that the first order topological p + ip and second order topological p + iτ p superconducting states emerge naturally in the presence of Dirac points near the chemical potential. While in a previous study we explored the pseudospin pairing mechanism giving rise to a p + iτ p instability in a semiconductor-based artificial honeycomb lattice within the RPA [41], in this work we have extended our analysis of the pairing mechanism to generic lattices satisfying C 6v symmetry and shown that this instability is not specific to any lattice structure but rather emerges due to the universal properties of the Dirac fermion excitations in these materials. We have illustrated that both the Berry phase and screening effects underlying the pairing mechanism are active in kagome lattices, and presented exact diagonalization results which demonstrate the striking similarity between the boundary physics in kagome and honeycomb lattices, which result from the fact that they possess identical topological properties.
In addition to the p + iτ p phase, the intervalley p + ip and s-wave phase are possible, but are subleading in our analysis. Unlike the p-wave states, the s-wave state is protected against the effects of weak disorder by a generalized Anderson's theorem, so this state might posses a higher T c than the topological phases in a sufficiently disordered system. In a forthcoming study we also suggest that this phase may arise due to spin-orbit coupling and study its properties in detail.
We call attention to a number of relevant features of this mechanism: (i) the mechanism can be applied for weak repulsive interactions, so strong correlations are not necessary, (ii) superconductivity appears beyond a critical doping beyond the Dirac points, at which the screening effects become strong enough to give rise to a net attraction, (iii) the electronic origin of the pairing means the cutoff which determines T c is the Fermi energy rather than the Debye frequency, (iv) superconductivity is more pronounced in systems with localized orbitals.
Superconductivity has recently been discovered in quasi two dimensional vanadium-based kagome metals AV 3 Sb 3 (A=K,Rb,Cs) [42,45]. The agreement between DFT calculations of the bandstructure and ARPES suggests that these materials are weakly correlated, an interpretation also supported by the results of DFT/DMFT calculations [42,65,74]. The Fermi surface consists of two pockets surrounding the Dirac points, a pockets surrounding the Γ point and an additional Fermi contour [66]. Superconductivity exists alongside density wave order which sets in at a higher temperature and appears to compete with superconductivity [67][68][69][70][71][72][73]. STM and Josephson STS measurements suggest a spatially modulated superconducting gap [47] and zero bias peaks inside magnetic field-induced vortex cores are suggestive of Majorana bound states [48].
Given that phonons appear unable to account for the measured T c [47,74] and yet the materials appear to be weakly correlated, we argue that our mechanism might provide an explanation for superconductivity in kagome metals, in contrast to existing theoretical proposals which attribute superconductivity to the effects of nesting or competing density wave order [75][76][77]. In this interpretation, superconductivity originates from the Dirac-like Fermi surfaces and is not directly related to the observed charge density wave order [67], which is potentially due to van Hove singularities or nested portions of the band structure near the Fermi level. Future theoretical work should investigate the interplay between the pairing mechanism and the presence of density wave order, which breaks rotational symmetry [68]. Future experiments could look for evidence of the edge or corner modes present in the p + ip and p + iτ p phases [22,23], probe the spin structure of the gap through NMR, or further investigate the real space structure of the superconducting gap.
Additional kagome systems have recently been discovered to host superconductivity, and might also be possible candidates for the pseudospin mechanism, including bilayer kagome systems [78], ferromagnetic kagome metals [79], and lanthanum-based materials [80]. Superconductivity has also been seen in Dirac surface states of doped topological insulators [82][83][84][85][86], which might also be explained by this theory. In materials such as transition metal dichalcogenides like MoS 2 and MoTe 2 , and fewlayer stanene [87], the presence of spin orbit interactions gaps the Dirac points and/or introduces a valley-Zeeman field, but the Dirac physics persists nonetheless. In a forthcoming paper, we have examined the effects of spinorbit coupling on the possible superconducting states, and find that the pseudospin mechanism results in topological superconductivity in these systems as well. As a general guiding principle, promising candidate materials for this mechanism are those with strongly localized orbitalswhich lead to larger values of the couplings g i and lower values of the critical doping µ * -and larger Dirac-like Fermi surfaces -which enhance the screening effects.

VI. ACKNOWLEDGEMENTS
The authors thank P.

A1. SYMMETRY RELATIONS FOR THE EFFECTIVE FIELD THEORY
We may explicitly derive the interaction parameters V ab in Eq. 3 of the main text by projecting the extended Hubbard interactions onto the valley and pseudospin eigenstates u τ α (r). We find The Hubbard interactions may be expressed in terms of the Wannier orbitals φ r (x) (with x being a 3D coordinate vector) via where V (x − x ) is the Coulomb interaction, and we have chosen the Wannier orbitals to be purely real. The interactions thus satisfy U (r 1 , r 2 , r 3 , r 4 ) = U (r 3 , r 2 , r 1 , r 4 ). We may use the representation of time reversal symmetry (Eq. (2) in the main text) T = ΩK with Ω = τ x α x to derive a relation Applying this to the factors u τ1α1 (r 1 ) and u τ3α3 (r 3 ) in (A1) we find We may express this in the form Since the adjoint basis elements J a are Hermitian and transform with definite sign under time reversal, (Ω † J a Ω) T = (Ω † J a Ω) * = T −1 J a T = ±J a , we find that V ab = 0 for any operator J a which is odd under time reversal. We may derive further symmetry constaints using the relations where Λ is a symmetry operation acting on coordinates and U Λ is its unitary representation. Applying a simultaneous transformation r i → Λr i in (A1) we find Since there always exists a symmetry operation under which J a and J b transform with opposite signs unless a = b, we find that V ab is a diagonal matrix.

A2. LANDAU-GINZBURG ANALYSIS OF THE SUPERCONDUCTING GAP SYMMETRY
The generic mean field Hamiltonian to account for all pairing possibilities, in the upper band, is The gap takes any of the forms presented in Table II, which contains the valley, spatial and spin structure. The spin structure employs standard notation of d µ -vectors; The analysis presented in the main text states that the gap functions can be chosen as simply one of the possible degenerate sets; for the gap structure discussed in Section IV it is important that we can choose the d-vector along one axis ∆ ∝ d y s y , and that a non-unitary choice of ∆ k is energetically penalized. We explicitly show that this is the case below. We compute the free energy to quartic order from (A8), and find where Γ are the Cooper channel Born amplitudes, explicitly given in (5), and the trace includes summation over momentum and frequency. Each gap structure entering Table II has a critical temperature set by the coupling λ, where λ = Γ τ τ τ τ for intravalley pairing, while λ = Γ +−+− ± Γ +−−+ for intervalley pairing. We restrict our attention to work within each degenerate set of gaps separately. In most cases, the valley structure is trivially traced out -the free energy reduces to separate cases based on the overall spin degeneracy d · s and degenerate chiral phase factors e ±iθ k . Using this logic, there are only five distinct cases to consider Case (iii) enumerates the two physically distinct valley structures that share the same spin structure and chirality, but themselves are not degenerate, i.e. have different T c . The differing valley structures of these gaps does not affect the contribution to F that is quartic in ∆; the valley structure enters through Γ and affects the value of T c but not the stability analysis. The brackets in (v) accounts for the four physically distinct gap structures which are degenerate, i.e. have the same value of λ and hence have same critical temperature, c.f. Table II. We combine these three gap functions into a vector (v 0 , v −2τ , u τ ±1 ), which will use below in the free energy analysis. Evaluating the free energy to quartic order in the gaps, with coefficients a i > 0, b > 0, i = 1, 2, 3. We summarize the main conclusions: • Case (i) is a s-wave singlet paired state. The free energy straightforwardly establishes stability of this phase.
• Case (ii) corresponds to an intervalley p + ip state, in which electrons from opposite valleys pair with a definite chirality. From F 2 [d 0 ± s 0 e ±iθ k ], we find that coexistence of opposite chiralities is energetically penalized so that the system spontaneously chooses one chirality, e.g. d 0 + s 0 e +iθ k , therefore breaking time reversal. This state exhibit first order topology with gapless edge states.
• Case (iii) concerns two triplet paired states for which the free energy has the same form; from F 3 [d · s] we find that non-unitary paring is penalised, i.e. d is purely real (or purely imaginary); the d vector spontaneously chooses a direction, breaking spin SU (2) symmetry, but the condensate has vanishing magnetization. The two superconducting states included in this case are: intravalley p + iτ p, in which electrons undergo p + ip pairing in one valley and p − ip pairing in the other, and intervalley s-wave, in which electrons from opposite valleys undergo triplet s-wave pairing. Both states respect time-reversal symmetry.
• Case (iv) examines spin-triplet p-wave intervalley pairing. From F 4 [d ± · s e ±iθ k ], similar to case (ii), the resulting state spontaneously breaks time reversal symmetry, resulting in a first order topological p + ip state e.g. d + e +iθ k . As in case (iii), we again find the d-vector is purely real (or imaginary), and hence this state breaks time reversal symmetry but does not support a magnetization.

A3. SCREENED INTERACTIONS
In this appendix we discuss the screened electron-electron interactions. We begin by presenting general expressions for the polarization operators, and making some general comments about regulating UV divergences in effective quantum field theories. We then explicitly calculate the intra and inter valley current-current susceptibilities. We present the formulae for the screened couplings, and plot the frequency dependence of the resulting scattering amplitudes for Cooper pairs. We conclude by discussing alternative regularizations for computing the susceptibilities.

A. Preliminaries
The RPA equations for the screened interactions are where Π αγ is the polarization operator, where the vertices J µ = {J 1 , . . . , J 10 } = {τ 0 α 0 , τ z α z , τ 0 α x , τ 0 α y , τ x α 0 , τ y α 0 , τ x α x , τ x α y , τ y α x , τ y α y }, as per Section ?? of the main text. We choose units in which v = 1, and upon performing the frequency integral by residues, manipulations presented in [A1] result in The first term -the interband polarization operator, denoted Π + -contributes when µ = 0, while the second termthe intraband polarization operator, denoted Π − -only contributes when µ = 0. Explicitly, we write Π µν + (ω, q) = Tr The polarization operator Π 11 is the standard charge polarization operator for graphene [A2-A4], which describes the screening of the density-density interaction α 0 τ 0 ⊗ α 0 τ 0 . The function Π 22 is the pseudospin polarization operator which describes the screening behaviour of the pseudospin-pseudospin interaction α z τ z ⊗ α z τ z which was first calculated in [A1] (denoted there as Π zz;00 due to different notation and a difference choice of basis for the pseudospin states). The remaining functions describe the screening of the intra and intervalley chiral currents. The current-current susceptibility for graphene has been discussed in [A5-A7], though to the best of our knowledge the intervalley current-current susceptibility has not been previously investigated. Below, it will be convenient for us to derive expressions for both functions, and present them in a form somewhat more general than those already in the literature.
All polarization operators may all be written in terms of four basic functions, for which we will introduce simple notation in what follows. Use the spacetime indices J µ = (τ 0 α 0 , τ 0 α i ) where i = x, y, and denote Π 22 as Π zz . Considering J µ = α i τ 0 and J ν = α j τ 0 , ie the chiral current screening Π 33 , Π 34 , Π 44 , we then write this function as Π ij and decompose into longitudinal and transverse parts: The chiral current polarization operator obeys an important identity due to electromagnetic gauge invariance. In the Dirac theory, the electromagnetic current operator J = α so that the vector potential of electromagnetism A appears in the Hamiltonian as the perturbation eJ · A. Gauge invariance, or equivalently the conservation of charge, can be shown to imply the Ward identity, From this we find Combining these, which gives us an expression for the longitudinal part of the current-current susceptibility in terms of the density-density response Comparing our below results with the expression for Π 00 cite [A1] shows our expressions satisfy this relation. We shall also derive an expression for the transverse part of the susceptibility in the dimensional regularization scheme, which therefore allows the current-current screening operator to be written entirely in terms of the density and pseudospin responses. As is common quantum field theories, the polarization operators are formally divergent quantities and require regularization. However, the physical origins of the divergences and the effects of regularization differ between the polarization operators. The dimensional regularization scheme has the effect of simply setting all UV contributions to zero, leaving just the effects of those degrees of freedom in the Dirac effective theory near the K-points. Since actual materials are UV-completed by a lattice, placing the theory on a lattice is a more physical regularization scheme, and gives the physical relation between the numerical values of the couplings computed by band structure methods and the couplings in the effective theory.
Firstly, the charge susceptibility Π 00 (ω, q) obeys the exact compressibility sum rule, Π 00 (ω = 0, q → 0) = −ν 0 where ν 0 is the density of states at the Fermi level, which means that cutoff-dependent quantities never appear in any sensible regularization scheme. For small frequencies and momenta, the regularized polarization operator must be the same as that calculated in dimensional regularization.
Second, the current susceptibilities Π ij (ω, q) obey the Ward identity discussed above, as a result of the fact that α is the current operator in the Dirac theory, which also excludes cutoff-dependent contributions. However, away from the K-points, the current operator is no longer given by α, which will result in UV-dependent contributions to Π ij (ω, q). In the lattice regularization scheme, this means that one can make changes to the lattice -for instance by modifying the dispersion near the Γ point -that will result in a different value for Π ij (ω, q), but these UV-contributions originate from physics away from the Dirac point where the Ward identity applies. Similar to Π 00 (ω, q), for small frequencies and momenta, the regularized Π ij (ω, q) must be the same in lattice and dimensional regularization.
The pseudospin susceptibility Π zz (ω, q), by contrast, obeys no such constraint and receives a constant and negative cutoff dependent contribution in a lattice regularization. We may distinguish between two types of cutoff dependencecontributions which originate from physics near the Dirac point, and contributions which originate from physics away from the Dirac point. While only the latter affects Π ij (ω, q), as we explained above, both types of UV-contributions appear in Π zz (ω, q); one way to see this is to use a hard cutoff Λ in the Dirac theory, where (restoring the Dirac velocity v) one finds Π 00 (ω = 0, q → 0) = −Λ/(2πv) + ν 0 . Thus by solely modifying physics near the Dirac pointfor instance by changing the velocity -these UV contributions which originate from physics near the Dirac point are changed. As a result, these contributions which appear in a lattice regularization are important even at small frequencies and momenta, and are described in the main text.
When working with an effective field theory, one always performs calculations in terms of UV-independent quantities. However, if one wishes to perform a microscopic calculation of the couplings g i in (??) -by computing the wavefunctions of a material for e.g. through band structure diagonalization and then taking matrix elements of the Coulomb interaction -the couplings which appear in the effective quantum field theory will be related to those values through lattice regularization, in which Π ij (ω, q) and Π zz (ω, q) receive contributions from physics away from the Dirac point, and Π zz (ω, q) receives an additional constant contribution from physics near the Dirac point.

B. Intravalley current-current polarization operator
We now turn to a calculation of Π ij (ω, q). We begin by calculating the interband part through dimensional regularisation, using the formulae Focus first on the denominator in Eq. (A16). Using the Schwinger-Feynman parametrization, where we now use relativistic notation l µ = (E, k), p µ = (ω, q), and l 2 = E 2 − k 2 , p 2 = ω 2 − q 2 . Shifting l → l − xp (bear in mind that this will affect the numerator as well), and Wick rotating E → iE, the expression becomes The corresponding numerator (before Wick rotation) is Now we substitute J µ , J ν = α i , α j with i, j = x, y and perform the pseudospin trace using the identity where ∆ = x(x − 1)p 2 . Using rotational symmetry to replace (l i ) 2 → 1 3 l 2 , we arrive at where we account for the valley and spin trace through a factor N = 4. Integrating using Eq. (A23), we find the imaginary part Similarly using Eq. (A22), the real part is calculated from Eq. (A29) to be The functions Π ij ⊥ (Ω, q) and Π ij (Ω, q) are respectively the transverse and longitudinal components of the polarization operator. From above, we see the their interband parts are and Calculating the density response J µ = J ν = α 0 through the formalism above, one may see that the above expression for Π ij satisfies the Ward identity (A20).

C. Intervalley polarization operator
Taking J µ = τ ± α i and J ν = τ ∓ α j with τ ± = 1 2 (τ x ± iτ y ) and performing the trace changes the sign of the terms containing factors of τ z . If we denote this intervalley polarization operator Π ij II then one finds the simple result In other words, the transverse and longitudinal response are reversed between intra and intervalley. The same trace relations show that Π 22 = Π 55 = Π 66 = Π zz , ie the screening of the operators τ ± is the same as that of τ z α z .
As is shown by examination of the formulae in the previous subsection, the scattering amplitudes have a roughly step-like behavior. Below, we plot the p-wave inter and intra valley on-shell couplings -λ =1 intra = Γ =1 τ τ τ τ (|k| − |p|, |k − p|) as well as λ =1 inter = Γ +−+− (|k| − |p|, |k − p|) − Γ +−−+ (|k| − |p|, |k − p|) -for kagome and honeycomb to illustrate this phenomenon. Treating the scattering amplitude as a step function in frequency space, the frequency dependent gap equation can therefore be straightforwardly solved as in the Anderson-Morel treatment of the electron-phonon problem [A11], the result being a renormalized coupling appearing in the exponential form of T c and an effective frequency cutoff, c.f. the treatment in S3 of [A1].