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

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–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–11]. Recently, a new classification of topological superconductors has emerged which relies on crystalline symmetries in addition to non-spatial ones [12–24]. 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\tau 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 [25–31]. In contrast to previously studied mechanisms for higher order topological superconductivity, we also do not require proximitization to an existing superconductor [32–41], negative-U Hubbard interactions [42], or bosonic fluctuations [43].

This pairing mechanism was previously studied in the context of an artificial honeycomb lattice [44]; 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 2 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 3 we discuss how the superconducting instabilities are affected by screening effects. In section 4 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₃Sb₅ (A = K, Rb, Cs), where superconductivity appears to be unconventional but correlations are weak [45–51].

We investigate the pairing instabilities of a two-dimensional lattice with hexagonal symmetry hosting Dirac points at the K points, which we describe via an extended Hubbard model [52–54]

$$\begin{align} H & = \sum_{\boldsymbol{r},\boldsymbol{r}^{^{\prime}}}{T(\boldsymbol{r}-\boldsymbol{r}^{^{\prime}}) c^\dagger_{\boldsymbol{r}} c_{\boldsymbol{r}^{^{\prime}}}} \nonumber \\ & \quad+ \frac{1}{2}\sum_{\boldsymbol{r}_1,\boldsymbol{r}_2,\boldsymbol{r}_3,\boldsymbol{r}_4}{U(\boldsymbol{r}_1,\boldsymbol{r}_2,\boldsymbol{r}_3,\boldsymbol{r}_4)(c^\dagger_{\boldsymbol{r}_3} c_{\boldsymbol{r}_1})(c^\dagger_{\boldsymbol{r}_4} c_{\boldsymbol{r}_2})} \end{align} \tag{ 1 }$$

where $c^\dagger_{\boldsymbol{r}} = (c^\dagger_{\boldsymbol{r},\uparrow},c^\dagger_{\boldsymbol{r},\downarrow})$ 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 $\boldsymbol{a}_1 = (a,0)$, $\boldsymbol{a}_2 = (\frac{a}{2},\frac{\sqrt{3}a}{2})$, and specialize to lattices possessing $C_{6\upsilon}$ 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(\boldsymbol{r}-\boldsymbol{r}^{^{\prime}})$ and $U(\boldsymbol{r}_1,\boldsymbol{r}_2,\boldsymbol{r}_3,\boldsymbol{r}_4)$ are real. The Hamiltonian is invariant under SU(2) spin rotations and possesses a time reversal symmetry satisfying $\mathcal{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 $\mathcal{T}$.

The resulting band structure features Dirac points at the K and K' points, which we distinguish by a valley index $\tau = \pm$. The eigenstates of the single particle Hamiltonian may be classified by their eigenvalues $e^{\frac{2\pi i}{3} \alpha}$ under threefold rotations, where $\alpha = \pm$. We shall refer to the α degree of freedom as pseudospin. At the Dirac points there are four degenerate energy eigenstates $|\tau,\alpha\rangle$ which transform into one another under time reversal, twofold rotations and mirror symmetry. We shall introduce the Pauli operators $\tau_\mu,\alpha_\mu$ 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_{\pi}\mathcal{T}$. We obtain the following representation of crystal symmetries:

$$\begin{align} & U_{\mathcal{R}_{\frac{2\pi}{3}}} = e^{\frac{2\pi i}{3} \alpha_z} , \ \ U_{M_x} = \alpha_x , \nonumber \\ & U_{\mathcal{R}_\pi} = \cos \phi_I \tau_x + \sin \phi_I \tau_y , \ \ \mathcal{T} = \tau_x\alpha_x \mathcal{K} \end{align} \tag{ 2 }$$

where φ_I is a phase which generally depends on the lattice and the chosen center of inversion.

Doping slightly above the valleys $\pm\boldsymbol{K} = \pm (\frac{4\pi}{3a},0)$, the Fermi surface consists of two circular pockets surrounding the K points with Fermi momenta $k_F \ll |\boldsymbol{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 eight-component local field operator $\psi^\dagger_{\tau\alpha s}(\boldsymbol{r})$. The most general Hamiltonian density consistent with the symmetries (2) is of the form

$$\begin{align} \mathcal{H} = \psi^\dagger (-i\upsilon \tau_z\boldsymbol{\alpha}\cdot\nabla)\psi +\frac{1}{2} \sum_{a,b}{\mathcal{V}_{ab} \left(\psi^\dagger J^{\,a} \psi\right)\left(\psi^\dagger J^{\,b}\psi\right)} \end{align} \tag{ 3 }$$

where we have expanded the interaction in the adjoint basis consisting of products $J^{\,a} = \tau_\mu \alpha_\nu$, and $\mathcal{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,\dots, J^{10}\} = \{\tau_0\alpha_0, \tau_z\alpha_z, \tau_0\alpha_x,\, \tau_0\alpha_y,\, \tau_x, \tau_y,\, \tau_x\alpha_x, \tau_x\alpha_y, \tau_y\alpha_x, \tau_y\alpha_y\}$. The Hamiltonian must transform as a scalar under the spatial symmetries (2) which implies that $\mathcal{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 $\mathcal{V}_{ab}$ by

$$\begin{align} \mathcal{V}_{11} = g_1 , \ \mathcal{V}_{22} = g_2 , \ \mathcal{V}_{33} = \mathcal{V}_{44} = \frac{g_3}{2} , \nonumber \\ \mathcal{V}_{55} = \mathcal{V}_{66} = \frac{g_4}{2} \ ,\ \mathcal{V}_{77} = \mathcal{V}_{88} = \mathcal{V}_{99} = \mathcal{V}_{10,10} = \frac{g_5}{4} . \end{align} \tag{ 4 }$$

In addition to a density–density interaction g₁, the interactions involve valley-conserving ($J^2 = \tau_z\alpha_z$) and valley-mixing ($J^5 = \tau_x, J^6 = \tau_y$) mass operators. The remaining six operators are conserved SU(2) valley currents which can be expressed as products of the valley operators $i\tau_x,i\tau_y,\tau_z$ and the electric current $j_\mu = (j_x,j_y) = (\tau_z\alpha_x,\tau_z\alpha_y)$.

The relation between the coupling constants g_i appearing in the field theory (3) and the extended Hubbard model parameters $U(\boldsymbol{r}_1,\boldsymbol{r}_2,\boldsymbol{r}_3,\boldsymbol{r}_4)$ for the kagome and honeycomb lattices are given in table 1. 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 $\boldsymbol{r}_1,\boldsymbol{r}_2,\boldsymbol{r}_3,\boldsymbol{r}_4$ exist on the $\sigma_1,\sigma_2,\sigma_3,\sigma_4$ sublattices by $U_{\sigma_1\sigma_2\sigma_3\sigma_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.

Table 1. The relations between coupling constants (3) and the extended Hubbard model parameters (1) for the Kagome and honeycomb lattice with couplings involving nearest neighbors. We denote $U_{\sigma_1\sigma_2\sigma_3\sigma_4} = U(\boldsymbol{r}_1,\boldsymbol{r}_2,\boldsymbol{r}_3,\boldsymbol{r}_4)$ where $\boldsymbol{r}_1,\boldsymbol{r}_2,\boldsymbol{r}_3,\boldsymbol{r}_4$ are sites in sublattices $\sigma_1,\sigma_2,\sigma_3,\sigma_4$ separated at most by a nearest neighbor bond.

	Kagome	Honeycomb
$g_1/\Omega$	$\dfrac{1}{6}\left(2U_{AAAA} + 8U_{ABAB}+ 8U_{AABA} + 4U_{ACBC} + U_{AABB}\right)$	$\dfrac{1}{2}\left(U_{AAAA} + 3~U_{ABAB}\right)$
$g_2/\Omega$	$\dfrac{1}{2} U_{AABB}$	$\dfrac{1}{2}\left(U_{AAAA} - 3~U_{ABAB}\right)$
$g_3/\Omega$	$\dfrac{1}{3}\left( U_{AAAA} - 2U_{ABAB} + 4U_{AABA} - 4U_{ACBC} + 2U_{AABB}\right)$	$3~U_{ABBA}$
$g_4/\Omega$	$\dfrac{1}{3}\left( U_{AAAA} - 2U_{ABAB} + 3U_{AABA} - 3U_{ACBC} + 2U_{AABB} \right)$	$3~U_{ABBA}$
$g_5/\Omega$	$\dfrac{1}{3}\left( U_{AAAA} + U_{ABAB} + 4U_{AABA} + 8U_{ACBC} + 8U_{AABB} \right)$	UAAAA

The eigenstates of the Dirac Hamiltonian in the upper band are given by $f^{\,\dagger}_{\boldsymbol{k},\tau,s} = (\psi^\dagger_{\tau \tau s}(\boldsymbol{k}) + \tau e^{i\tau \theta_{\boldsymbol{k}}} \psi^\dagger_{\tau -\tau s}(\boldsymbol{k}))/\sqrt{2}$. Evaluating the Born amplitudes from (3) we find that the scattering vertex $\Gamma_{\tau_1\tau_2\tau_3\tau_4}(\theta)$ in the Cooper channel is given by

$$\begin{align} &\Gamma_{\tau\tau\tau\tau}(\theta) = \frac{g_1+g_2}{2} e^{-i\tau \theta}\cos\theta + \frac{g_1-g_2-g_3}{2} e^{-i\tau\theta}, \nonumber \\ &\Gamma_{+-+-}(\theta) = \Gamma_{-+-+}(\theta) = \frac{g_1 -g_2}{2}\cos\theta\nonumber \\ & \qquad \qquad \qquad \qquad \qquad + \frac{g_1+g_2+g_3}{2}, \nonumber \\ &\Gamma_{+-+}(\theta) = \Gamma_{-++-}(\theta) = -\frac{g_4}{2}\cos\theta + \frac{g_4+g_5}{2}, \end{align} \tag{ 5 }$$

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 $\Gamma^\ell_{\tau\tau\tau\tau} \lt 0$, while pairing of electrons in opposite valleys requires the symmetrized amplitudes $\Gamma^\ell_{+-+-} + \Gamma^\ell_{+-+} \lt 0$ or $\Gamma^\ell_{+-+-} - \Gamma^\ell_{-++-} \lt 0$. The critical temperature is given

$$\begin{align} T_c \sim E_F e^{-1/(\nu_0\lambda)} \end{align} \tag{ 6 }$$

where $\lambda = \Gamma^\ell_{\tau\tau\tau\tau}$ for intravalley pairing or $\lambda = \Gamma^{\ell}_{+-+-} \pm \Gamma^{\ell}_{+-+}$ for intervalley pairing, and $\nu_0 = k_F/(2\pi \nu)$ 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

$$\begin{align} H_\textrm{BdG}& = \sum_{\boldsymbol{k}}\varepsilon_{\boldsymbol{k}} f^{\,\dagger}_{\boldsymbol{k}}f_{\boldsymbol{k}} + \frac{1}{2}f^{\,\dagger}_{\boldsymbol{k}} \left(\Delta_{\boldsymbol{k}} i s_y i \tau_y\right) f^{\,\dagger}_{-\boldsymbol{k}} + \textrm{h.c.} \end{align} \tag{ 7 }$$

where $f^{\,\dagger}_{\boldsymbol{k}}$ is the creation operator for a Dirac fermion in the upper band and $\Delta(\boldsymbol{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 2 along with the corresponding scattering amplitude.

Table 2. Superconducting phases for 2D Dirac fermions in lattices with $C_{6\upsilon}$ symmetry. First column: gap structure, second column: order parameter, third column: coupling constant which determines T_c via equation (6), in terms of coupling constants in the Hamiltonian equation (3), fourth column: Altland–Zirnbauer class.

Gap structure	Δ	λ	AZ class
Intervalley, spin singlet, s-wave	τz	$\dfrac{1}{2}\left(g_1+g_2+g_3+g_4+g_5\right)$	CI
Intervalley, spin singlet, p-wave	$e^{\pm i \theta_{\boldsymbol{k}}}$	$\dfrac{1}{4}\left(g_1-g_2+g_4\right)$	C
Intervalley, spin triplet, s-wave	$d^\mu s_\mu$	$\dfrac{1}{2}\left(g_1+g_2+g_3-g_4-g_5\right)$	BDI
Intervalley, spin triplet, p-wave	$e^{\pm i \theta_{\boldsymbol{k}}} \tau_z (d^\mu s_\mu)$	$\dfrac{1}{4}\left(g_1 - g_2 - g_4\right)$	D
Intravalley, spin singlet, s-wave	$e^{i\tau_z\phi} (i\tau_y)$	$\dfrac{1}{4}\left(g_1+g_2\right)$	CI
Intravalley, spin singlet, d-wave	$e^{i\tau_z(\phi - 2\theta_{\boldsymbol{k}})} (i\tau_y)$	$\dfrac{1}{4}\left(g_1+g_2\right)$	CI
Intravalley, spin singlet, mixed s/d-wave	$e^{\pm i\theta_{\boldsymbol{k}}} e^{i\tau_z(\phi+\theta_{\boldsymbol{k}})}(i\tau_y)$	$\dfrac{1}{4}\left(g_1+g_2\right)$	C
Intravalley, spin triplet, p-wave	$e^{i\tau_z(\phi-\theta_{\boldsymbol{k}})} (d^\mu s_\mu)(i\tau_y)$	$\dfrac{1}{2}\left(g_1-g_2-g_3\right)$	BDI

The topological properties of the gaps are associated with their Altland–Zirnbauer class [55–57], which classifies mean field Hamiltonians based on time reversal, charge conjugation and chiral symmetries. For the spin singlet pairing phases, $\Delta \propto s^0$, accounting for SU(2) spin rotation symmetry we find that the charge conjugation symmetry satisfies $\mathcal{C}^2 = -1$. In the case of spin triplet pairing, the spontaneous violation of SU(2) spin rotation symmetry fixes a direction $\boldsymbol{d} = (d^{\,x}, d^{\,y}, d^{\,z})$, with the gap $\Delta \propto d^\mu s_\mu$. 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 ⁶ , we may decompose the mean field Hamiltonian in two spin blocks, with each separately possessing a charge conjugation symmetry satisfying $\mathcal{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 [57], 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 possess 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\tau 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 $\mathcal{T}^2 = +1$, our system possesses a spinful time reversal symmetry $\mathcal{T}^{\prime} = \mathcal{T} e^{i\pi S_y}$ satisfying $(\mathcal{T}^{\prime})^2 = -1$ which interchanges opposite spin blocks; thus corner modes occur in Kramers pairs. In section 4 we present exact diagonalization results which demonstrate that the intravalley $p+i\tau 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(\boldsymbol{r}_1,\boldsymbol{r}_2,\boldsymbol{r}_3,\boldsymbol{r}_4)$ in (1) with $\boldsymbol{r}_1 = \boldsymbol{r}_3,\boldsymbol{r}_2 = \boldsymbol{r}_4$. Denoting the parameter λ in (6) by λ₁ for intervalley s-wave spin triplet pairing, λ₂ for intervalley p-wave spin triplet pairing and λ₃ for intravalley p-wave spin triplet pairing, we find for the kagome lattice

$$\begin{align} \lambda_1 = \lambda_2 = \frac{\lambda_3}{2} = \frac{U_{ABAB}\Omega}{2} \end{align} \tag{ 8 }$$

where Ω is the area of the unit cell and

$$\begin{align} \lambda_1 = 0, \ \lambda_2 = \frac{\lambda_3}{2} = \frac{3~U_{ABAB}\Omega}{4} \end{align} \tag{ 9 }$$

for the honeycomb lattice.

In the limit $U_{ABAB}\rightarrow 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\tau 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 2. This causes destructive interference between the contributions to the scattering amplitude with distinct pseudospin and valley structure—i.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.

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 [59–62]. This analysis was performed for the special case of an artificial semiconductor honeycomb lattice in [44]; we extend this analysis to encompass generic Dirac materials satisfying $C_{6\upsilon}$ and time-reversal symmetry, which may be described by the Hamiltonian (3). The RPA results in replacing the interaction constants $\mathcal{V}_{ab}$ in (3) with screened frequency and momentum dependent couplings, which are determined by the RPA equation

$$\begin{align} \widetilde{\mathcal{V}}_{ab}(\omega,\boldsymbol{q}) = \mathcal{V}_{ab} + \sum_{cd}{\mathcal{V}_{ac}\Pi^{cd}(\omega,\boldsymbol{q}) \widetilde{\mathcal{V}}_{db}(\omega,\boldsymbol{q})}, \end{align} \tag{ 10 }$$

where $\Pi^{cd}(\omega,\boldsymbol{q})$ is a generalized polarization operator. These represent the susceptibility of the system to a perturbation $\delta \mathcal{H} = \sum_a{\phi_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

$$\begin{align} T_c = \Lambda e^{-1/(\nu_0~\widetilde{\lambda})} \end{align} \tag{ 11 }$$

where $\widetilde{\lambda}$ is the RPA renormalised λ of table 2 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 $\Lambda \approx E_F$.

Scattering on the Fermi surface corresponds to the range of momentum and frequency transfer $q \lt 2k_F$, ω = 0. In this range the polarization operators are constant, and $\widetilde{\lambda}$ is simply given by the previous expressions λ (table 2) evaluated with screened couplings $\widetilde{g}_i$ at $\omega = 0,q = 0$. We have studied the screening of the interactions g₁ and g₂ in a previous paper [44], and found that g₁ is reduced while g₂ 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₁ is screened by the long wavelength density–density response, which as usual weakens g₁ with increasing chemical potential. The intra and intervalley mass couplings g₂ and g₄ are screened via the polarization operator associated with the response of the system to an external perturbation $\propto \alpha_z\tau_z, \tau_{x}, \tau_{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 regularization—the result of which is that g₂ and g₄ should be replaced by reduced interaction constants which are different from those provided in table 1; 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

$$\begin{align} \widetilde{g}_1 & = \frac{g_1}{1 + N\nu_0~g_1} \ , \nonumber\\ \widetilde{g}_i & = \frac{g_i}{1 - N\nu_0~g_i} \ , \nonumber\\ \widetilde{g}_j & = g_j \end{align} \tag{ 12 }$$

where $i = 2,4$ and $j = 3,5$. This results in replacing the nearest-neighbor coupling with a screened value, $U_{ABAB}\rightarrow \widetilde{U}_{ABAB}$ in the relations in table 1, equations (8) and (9). The screened nearest neighbor interaction is given by

$$\begin{align} \widetilde{U}_{ABAB} = \frac{\widetilde{g}_1-\widetilde{g}_3}{2\Omega} = \frac{1}{2\Omega}\left[ \frac{g_1}{1 + N\nu_0~g_1} - g_3\right] \end{align} \tag{ 13 }$$

for the kagome lattice and

$$\begin{align} \widetilde{U}_{ABAB} = \frac{\widetilde{g}_1-\widetilde{g}_2}{3\Omega} = \frac{1}{3\Omega}\left[\frac{g_1}{1 + N\nu_0~g_1} - \frac{g_2}{1 - N \nu_0~g_2}\right] \end{align} \tag{ 14 }$$

for the honeycomb lattice. For sufficiently large ν₀, we find that $\widetilde{U}_{ABAB} \lt0$, resulting in superconductivity, and the dominant instability occurs in the intravalley, p-wave, spin triplet channel. This occurs for $\mu\gt \mu^*$ with the critical value

$$\begin{align} \mu^* = 2\pi \upsilon^2\left[ \frac{g_3^{-1}-g_1^{-1}}{N}\right] \end{align} \tag{ 15 }$$

for the kagome lattice and

$$\begin{align} \mu^* = 2\pi \upsilon^2\left[\frac{g_2^{-1} - g_1^{-1}}{2N}\right] \end{align} \tag{ 16 }$$

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 $\mu^*$ lies within the regime where Dirac theory is justified. In the limit $g_3\rightarrow 0$ or $g_2\rightarrow 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 $\widetilde{g}_1$ with its screened value at $q\rightarrow 0$, $\widetilde{g}_1 \rightarrow 1/N\nu_0$, and simply results in a shift of the critical doping to higher values.

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 gap structure $\Delta_{\boldsymbol{k}} = |\Delta| e^{i\tau_z(\phi-\theta_{\boldsymbol{k}})} (d^\mu s_\mu)(i\tau_y)$. 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\tau 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 [63–70]. 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 $\phi = \frac{n\pi}{2}$, $n\in \mathbb{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\tau 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 $\boldsymbol{d}\parallel \hat{\boldsymbol{y}}$, 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\tau p$ momentum space structure. The spinless mean field lattice Hamiltonian which results is

$$\begin{align} H_{\textrm{BdG}} = \sum_{\langle \boldsymbol{r},\boldsymbol{r}^{^{\prime}}\rangle}{-tc^\dagger_{\boldsymbol{r}} c_{\boldsymbol{r}^{^{\prime}}} + \frac{1}{2}\left[\Delta(\boldsymbol{r},\boldsymbol{r}^{^{\prime}}) c^\dagger_{\boldsymbol{r}}c^\dagger_{\boldsymbol{r}^{^{\prime}}} +\textrm{H.c.}\right]} , \end{align} \tag{ 17 }$$

where $\Delta(\boldsymbol{r}^{^{\prime}},\boldsymbol{r}) = -\Delta(\boldsymbol{r},\boldsymbol{r}^{^{\prime}})$, and

$$\begin{align} \Delta(\boldsymbol{r},\boldsymbol{r}^{^{\prime}}) = \Delta^{^{\prime}}_{}\cos(\boldsymbol{K}\cdot(\boldsymbol{r}+\boldsymbol{r}^{^{\prime}}) + \phi) \ , \ \begin{cases}\boldsymbol{r} \in B, \ \boldsymbol{r}^{^{\prime}}\in C \\ \boldsymbol{r}\in B, \ \boldsymbol{r}^{^{\prime}}\in A \\ \boldsymbol{r}\in A, \ \boldsymbol{r}^{^{\prime}}\in C \end{cases} \end{align} \tag{ 18 }$$

for the kagome lattice, and

$$\begin{align} \Delta(\boldsymbol{r},\boldsymbol{r}^{^{\prime}}) = \Delta^{^{\prime}}_{}\sin(\boldsymbol{K}\cdot(\boldsymbol{r}+\boldsymbol{r}^{^{\prime}}) + \phi), \ \boldsymbol{r}\in B, \boldsymbol{r}^{^{\prime}}\in A \end{align} \tag{ 19 }$$

for the honeycomb lattice, and we distinguish the parameter $\Delta^{^{\prime}}_{}$ from the bulk gap Δ.

We present exact diagonalization results for the kagome and honeycomb lattices in figures 1 and 2 respectively. We have calculated the spectrum for all values of φ, however the topological classification must be performed at the values $\phi = 0, \frac{\pi}{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 figure 1 and the last row of figure 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 $\phi = \frac{\pi}{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 figures 1, 2(c), (f) and (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 $\phi = \frac{\pi}{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 $\mathcal{T}^{^{\prime}} = e^{i\pi S_y} \mathcal{T}$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 $\boldsymbol{K}\cdot(\boldsymbol{r}+\boldsymbol{r}^{^{\prime}}) + \phi$ of the pair density wave does not depend on the y coordinate. However, at discrete values $\phi = \frac{\pi n}{2}$, the superconductor recovers a twofold rotational symmetry. The parameter φ allows continuous tuning of the system between a second order topological phase at $\phi = \pi(n+\frac{1}{2})$ and a trivial phase at $\phi = n\pi$. 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 $|\phi - \pi(n+\frac{1}{2})| \lt \phi^*$ up to a critical value $\phi^*$ 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\tau 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 figures 1 and 2, the system possesses a mirror symmetry at $\phi = \frac{\pi 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 separate study, we show that this phase may become dominant when spin orbit coupling is present, and analyze the topology of this phase in detail [71].

Our results show that the first order topological p + ip and second order topological $p+i\tau 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\tau p$ instability in a semiconductor-based artificial honeycomb lattice within the RPA [44], in this work we have extended our analysis of the pairing mechanism to generic lattices satisfying $C_{6\upsilon}$ 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\tau p$ phase, the intervalley p + ip and s-wave phase are possible, but are subleading in our analysis. The intervalley p + ip is a first order topological supercondutor with protected edge modes; the mechanism we describe here therefore can also result in first order topological superconductivity. 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 separate study, we also suggest that these two phases may arise due to spin–orbit coupling and study their properties in detail [71].

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-2D vanadium-based kagome metals AV₃Sb₅ (A = K, Rb, Cs) [45, 48]. 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 [45, 72, 73]. The Fermi surface contains pockets surrounding the Dirac points, a pocket surrounding the Γ point, and nearly-hexagonal Fermi contours [74]. Superconductivity exists alongside density wave order which sets in at a higher temperature and appears to compete with superconductivity [75–81]. STM and Josephson STS measurements suggest a spatially modulated superconducting gap [50] and zero bias peaks inside magnetic field-induced vortex cores are suggestive of Majorana bound states [51]. These materials also show signs of a phase transition between two distinct superconducting states as a function of doping [82].

Given that phonons appear unable to account for the measured T_c [50, 73] and yet the materials appear to be weakly correlated, we propose that our mechanism might provide a partial 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 [83–85]. In this interpretation, superconductivity originates from the Dirac-like Fermi surfaces and is not directly related to the observed charge density wave (CDW) order [75], which is potentially due to the nested portions of the band structure near the Fermi level. Due to the presence of multiple Fermi surfaces as well as Fermi surface reconstruction due to CDW order, our theory does not fully incorporate all the possible interactions. Given the apparent phase transition between two superconducting phases as doping is varied, we suggest that our theory may describe one of the observed phases, with different effects responsible for the other. Our theory can also not explain the observed periodicity of the pair density wave, however it is possible that this discrepancy is related to the reconstruction of the Fermi surface by the CDW, which would modify the spatial period predicted by our theory. Future experiments could look for evidence of the edge or corner modes present in the p + ip and $p+i\tau p$ phases [23, 24], probe the spin structure of the gap through NMR or further investigate the real space structure of the superconducting gap as a function of doping/pressure.

Additional kagome systems have recently been discovered to host superconductivity, and might also be possible candidates for the pseudospin mechanism, including bilayer kagome systems [86], ferromagnetic kagome metals [87], and lanthanum-based materials [88]. Superconductivity has also been seen in Dirac surface states of doped topological insulators [89–93], which might also be explained by this theory. In materials such as transition metal dichalcogenides like MoS₂ and MoTe₂, and few-layer stanene [94], 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 separate study, we have examined the effects of spin–orbit coupling on the possible superconducting states, and find that the pseudospin mechanism results in topological superconductivity in these systems as well [71]. Our mechanism may also play a role in the observed superconducting phases in Moiré systems such as twisted bilayer/trilayer graphene [95–98]. These systems exhibit honeycomb/kagome superlattices hosting Dirac points, however, while our theory is expected to be applicable, the interaction structure is significantly more complicated than for the monolayer honeycomb or kagome materials we have investigated. As a general guiding principle, promising candidate materials for the pseudospin mechanism are those with strongly localized orbitals—which lead to larger values of the couplings g_i and lower values of the critical doping $\mu^*$—and larger Dirac-like Fermi surfaces—which enhance the screening effects.

The authors thank B Ortiz, P Brouwer, M Scheurer and P Kotetes for helpful discussions. M G acknowledges support by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program under Grant Agreement No. 856526, and from the Deutsche Forschungsgemeinschaft (DFG) Project Grant 277101999 within the CRC network TR 183 'Entangled States of Matter' (subproject A03 and C01), and from the Danish National Research Foundation, the Danish Council for Independent Research | Natural Sciences. H D Scammell acknowledges funding from ARC Centre of Excellence FLEET. T L acknowledges support from the DFG within the CRC network TR 183.

All data that support the findings of this study are included within the article (and any supplementary files).

Appendix. Symmetry relations for the effective field theory

We may explicitly derive the interaction parameters $\mathcal{V}_{ab}$ in equation (3) of the main text by projecting the extended Hubbard interactions onto the valley and pseudospin eigenstates $u_{\tau\alpha}(\boldsymbol{r})$. We find

$$\begin{align} & \sum{\mathcal{V}_{ab} J^{\,a}_{\tau_1\tau_3;\alpha_1\alpha_3} J^{\,b}_{\tau_2\tau_4;\alpha_2\alpha_4}}\nonumber\\ & \;\; = \sum_{\boldsymbol{r}_i}U(\boldsymbol{r}_1,\boldsymbol{r}_2,\boldsymbol{r}_3,\boldsymbol{r}_4) u^*_{\tau_1\alpha_1}(\boldsymbol{r}_1) u_{\tau_3\alpha_3}(\boldsymbol{r}_3) \nonumber\\ & \qquad \times u^*_{\tau_2\alpha_2}(\boldsymbol{r}_2) u_{\tau_4\alpha_4}(\boldsymbol{r}_4). \end{align} \tag{ A1 }$$

The Hubbard interactions may be expressed in terms of the Wannier orbitals $\phi_{\boldsymbol{r}}(x)$ (with x being a 3D coordinate vector) via

$$\begin{align} U(\boldsymbol{r}_1,\boldsymbol{r}_2,\boldsymbol{r}_3,\boldsymbol{r}_4)& = \int V(x-x^{^{\prime}}) \phi_{\boldsymbol{r}_1}(x) \phi_{\boldsymbol{r}_3}(x) \phi_{\boldsymbol{r}_2}(x^{^{\prime}})\nonumber\\ & \quad \times \phi_{\boldsymbol{r}_4}(x^{^{\prime}}) d^3 x d^3 x^{^{\prime}}, \end{align} \tag{ A2 }$$

where $V(x-x^{^{\prime}})$ is the Coulomb interaction, and we have chosen the Wannier orbitals to be purely real. The interactions thus satisfy $U(\boldsymbol{r}_1,\boldsymbol{r}_2,\boldsymbol{r}_3,\boldsymbol{r}_4) = U(\boldsymbol{r}_3,\boldsymbol{r}_2,\boldsymbol{r}_1,\boldsymbol{r}_4)$.

We may use the representation of time reversal symmetry (equation (2) in the main text) $\mathcal{T} = \Omega \mathcal{K}$ with $\Omega = \tau_x\alpha_x$ to derive a relation

$$\begin{align} u^*_{\tau\alpha}(\boldsymbol{r}) = \sum_{\tau^{^{\prime}},\alpha^{^{\prime}}}{u_{\tau^{^{\prime}}\alpha^{^{\prime}}}(\boldsymbol{r}) \Omega_{\tau^{^{\prime}}\tau;\alpha^{^{\prime}}\alpha}}. \end{align} \tag{ A3 }$$

Applying this to the factors $u_{\tau_1\alpha_1}(\boldsymbol{r}_1)$ and $u_{\tau_3\alpha_3}(\boldsymbol{r}_3)$ in (A1) we find

$$\begin{align} & \sum{\mathcal{V}_{ab} J^{\,a}_{\tau_1\tau_3;\alpha_1\alpha_3} J^{\,b}_{\tau_2\tau_4;\alpha_2\alpha_4}} \nonumber\\ & \quad = \Omega_{\tau^{^{\prime}}_1\tau_1;\alpha^{^{\prime}}_1\alpha_1} \Omega^*_{\tau^{^{\prime}}_3\tau_3;\alpha^{^{\prime}}_3\alpha_3}\nonumber\\ & \qquad \times \sum{\mathcal{V}_{ab} J^{\,a}_{\tau^{^{\prime}}_3\tau^{^{\prime}}_1;\alpha^{^{\prime}}_3\alpha^{^{\prime}}_1} J^{\,b}_{\tau_2\tau_4;\alpha_2\alpha_4}}. \end{align} \tag{ A4 }$$

We may express this in the form

$$\begin{align} \sum_{ab}{}\mathcal{V}_{ab} J^{\,a}\otimes J^{\,b} = \sum_{ab}{}\mathcal{V}_{ab} (\Omega^\dagger J^{\,a} \Omega)^T \otimes J^{\,b}. \end{align} \tag{ A5 }$$

Since the adjoint basis elements J^a are Hermitian and transform with definite sign under time reversal, $(\Omega^\dagger J^{\,a} \Omega)^T = (\Omega^\dagger J^{\,a} \Omega)^* = \mathcal{T}^{-1} J^{\,a} \mathcal{T} = \pm J^{\,a}$, we find that $\mathcal{V}_{ab} = 0$ for any operator J^a which is odd under time reversal.

We may derive further symmetry constraints using the relations

$$\begin{align} u_{\tau\alpha}(\Lambda^{-1}\boldsymbol{r}) = u_{\tau^{^{\prime}}\alpha^{^{\prime}}}(\boldsymbol{r}) U^\Lambda_{\tau^{^{\prime}}\tau;\alpha^{^{\prime}}\alpha} \end{align} \tag{ A6 }$$

where Λ is a symmetry operation acting on coordinates and $U^\Lambda$ is its unitary representation. Applying a simultaneous transformation $\boldsymbol{r}_i \rightarrow \Lambda \boldsymbol{r}_i$ in (A1) we find

$$\begin{align} \sum_{ab}{\mathcal{V}_{ab} J^{\,a} \otimes J^{\,b} } = \sum_{ab}{\mathcal{V}_{ab} ((U^\Lambda)^\dagger J^{\,a} U^\Lambda)\otimes ((U^\Lambda)^\dagger J^{\,b} U^\Lambda)}. \end{align} \tag{ A7 }$$

Since there always exists a symmetry operation under which J^a and J^b transform with opposite signs unless a = b, we find that $\mathcal{V}_{ab}$ is a diagonal matrix.

Appendix. Landau–Ginzburg analysis of the superconducting gap symmetry

The generic mean field Hamiltonian to account for all pairing possibilities, in the upper band, is

$$\begin{align} H_{BdG}& = \sum_{\boldsymbol{k}, s, \tau}\varepsilon_{\boldsymbol{k}} f^{\,\dagger}_{\boldsymbol{k} s \tau}f_{\boldsymbol{k} s \tau}\nonumber\\[5pt] & \quad + \frac{1}{2}\sum_{\boldsymbol{k},s, \tau,s^{^{\prime}}, \tau^{^{\prime}}} f^{\,\dagger}_{\boldsymbol{k} s \tau} \left(\Delta_{\boldsymbol{k}} i s_y i \tau_y\right)_{s \tau, s^{^{\prime}} \tau^{^{\prime}}} f^{\,\dagger}_{-\boldsymbol{k} s^{^{\prime}} \tau^{^{\prime}}} + \textrm{h.c.} \nonumber\\[5pt] & \quad-\frac{1}{2}\sum_{\boldsymbol{k},\boldsymbol{p}}\Delta_{\boldsymbol{k}}^\dagger \Gamma^{-1}(\boldsymbol{k} ,\boldsymbol{p})\Delta_{\boldsymbol{p}}. \end{align} \tag{ A8 }$$

The gap takes any of the forms presented in table 2, which contains the valley, spatial and spin structure. The spin structure employs standard notation of d^µ -vectors;

$$\begin{align} & d^{\,x} = \frac{1}{2}\left(|{\uparrow\uparrow}\rangle-|{\downarrow\downarrow}\rangle\right), \ \ d^{\,y} = \frac{1}{2i}\left(|{\uparrow\uparrow}\rangle+|{\downarrow\downarrow}\rangle\right), \nonumber\\[5pt] & d^{\,z} = \frac{1}{2}\left(|{\uparrow\downarrow}\rangle+|{\downarrow\uparrow}\rangle\right), \ \ d^{\,0} = \frac{1}{2}\left(|{\uparrow\downarrow}\rangle-|{\downarrow\uparrow}\rangle\right). \end{align} \tag{ A9 }$$

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 4 it is important that we can choose the d-vector along one axis $\Delta\propto d_y s_y$, and that a non-unitary choice of $\Delta_{\boldsymbol{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

$$\begin{align} {\cal F}& = \frac{1}{2}\sum_{\boldsymbol{k},\boldsymbol{p}}\Delta_{\boldsymbol{k}}^\dagger \Gamma^{-1}(\boldsymbol{k} ,\boldsymbol{p})\Delta_{\boldsymbol{p}}+\frac{1}{4}\textrm{Tr} \left[\frac{\Delta_{\boldsymbol{k}} \Delta_{\boldsymbol{k}}^\dagger}{\omega^2+\epsilon_k^2}\right]\nonumber\\[5pt] & \quad + \frac{1}{8}\textrm{Tr} \left[\frac{\Delta_{\boldsymbol{k}} \Delta_{\boldsymbol{k}}^\dagger\Delta_{\boldsymbol{k}} \Delta_{\boldsymbol{k}}^\dagger}{(\omega^2+\epsilon_k^2)^2}\right] \end{align} \tag{ A10 }$$

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 2 has a critical temperature set by the coupling λ, where $\lambda = \Gamma^\ell_{\tau\tau\tau\tau}$ for intravalley pairing, while $\lambda = \Gamma^{\ell}_{+-+-} \pm \Gamma^{\ell}_{+-+}$ 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 $\boldsymbol{d}\cdot \boldsymbol{s}$ and degenerate chiral phase factors $e^{\pm i\theta_{\boldsymbol{k}}}$. Using this logic, there are only five distinct cases to consider

$$\begin{align} & (\textrm{i}) \ \ d^{\,0}~s_0 \ \tau_z \nonumber\\[5pt] & (\textrm{ii}) \ \ d^{\,0}_{\pm} s_0 \ e^{\pm i \theta_{\boldsymbol{k}}} \tau_0 \nonumber\\[5pt] & (\textrm{iii}) \ \ \boldsymbol{d}\cdot \boldsymbol{s} ,\ \ \boldsymbol{d}\cdot \boldsymbol{s}\ e^{i\tau_z(\phi-\theta_{\boldsymbol{k}})} \ (i\tau_y)\nonumber \\[5pt] & (\textrm{iv}) \ \ \boldsymbol{d}_\pm\cdot \boldsymbol{s} \ e^{\pm i \theta_{\boldsymbol{k}}} \ \tau_z \nonumber \\[5pt] & (\textrm{v}) \ \ d^{\,0}~s_0~\left(1, e^{\pm i \theta_{\boldsymbol{k}}} , e^{-2i\tau_z \theta_{\boldsymbol{k}}}e^{i\tau_z \theta_{\boldsymbol{k}}} \right) e^{i\tau_z\phi}(i\tau_y)\nonumber\\[5pt] & \qquad \equiv (\upsilon_0, \upsilon_{-2\tau}, u_{\tau\pm1}). \end{align} \tag{ A11 }$$

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 $\mathcal{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, cf table 2. We combine these three gap functions into a vector $(\upsilon_0,\upsilon_{-2\tau},u_{\tau\pm 1})$, which will use below in the free energy analysis.

Evaluating the free energy to quartic order in the gaps,

$$\begin{align} {\cal F}_1[d^{\,0}~s_0]& = a_1(T-T_{c_1})|d_0|^2 + b |d_0|^4 ,\\ {\cal F}_2[ d^{\,0}_{\pm} s_0 \ e^{\pm i \theta_{\boldsymbol{k}}}]& = a_2(T-T_{c_2})(|d_{0+}|^2+|d_{0-}|^2) + b \left(|d_{0+}|^4 + 4 |d_{0+}|^2 |d_{0-}|^2 + |d_{0-}|^4\right),\\ {\cal F}_3[ \boldsymbol{d}\cdot \boldsymbol{s} ]& = a_3(T-T_{c_3})(\boldsymbol{d}^*\cdot \boldsymbol{d}) + b \left\{(\boldsymbol{d}^*\cdot \boldsymbol{d})^2+|\boldsymbol{d}^*\times \boldsymbol{d}|^2\right\} ,\\ {\cal F}_4[ \boldsymbol{d}_\pm\cdot \boldsymbol{s} \ e^{\pm i \theta_{\boldsymbol{k}}} ]& = a_4(T-T_{c_4})(\boldsymbol{d}_+^*\cdot \boldsymbol{d}_+ + \boldsymbol{d}_-^*\cdot \boldsymbol{d}_-) + b \Big\{(\boldsymbol{d}_+^*\cdot \boldsymbol{d}_+)^2+(\boldsymbol{d}_-^*\cdot \boldsymbol{d}_-)^2+4(\boldsymbol{d}_+^*\cdot \boldsymbol{d}_+)(\boldsymbol{d}_-^*\cdot \boldsymbol{d}_-)\\ & \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad +|\boldsymbol{d}_+^*\times \boldsymbol{d}_+|^2 + |\boldsymbol{d}_-^*\times \boldsymbol{d}_-|^2\\ & \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad + 4(\boldsymbol{d}_+^*\times \boldsymbol{d}_+)(\boldsymbol{d}_-\times \boldsymbol{d}_-^*)\Big\}\\ {\cal F}_5[(\upsilon_0, \upsilon_{-2\tau}, u_{\tau\pm1})]& = a_5(T-T_{c_5})\left(|\upsilon_0|^2 + |\upsilon_{-2\tau}|^2+|u_{\tau+1}|^2+|u_{\tau-1}|^2\right)\\ & \quad + b\Big\{|\upsilon_0|^4 + |\upsilon_{-2\tau}|^4+|u_{\tau+1}|^4+|u_{\tau-1}|^4 + 4|\upsilon_0|^2\Big( |\upsilon_{-2\tau}|^2+ |u_{\tau+1}|^2+ |u_{\tau-1}|^2\Big)\\ & \qquad\quad +4|\upsilon_{-2\tau}|^2\Big(|u_{\tau+1}|^2+ |u_{\tau-1}|^2\Big) + 4|u_{\tau+1}|^2|u_{\tau-1}|^2 \Big\} \end{align}$$

with coefficients $a_i\gt0, b\gt0, 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 ${\cal F}_2[ d^{\,0}_{\pm} s_0 e^{\pm i \theta_{\boldsymbol{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 \theta_{\boldsymbol{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 ${\cal F}_3[ \boldsymbol{d}\cdot \boldsymbol{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\tau 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 ${\cal F}_4[ \boldsymbol{d}_\pm\cdot \boldsymbol{s} \ e^{\pm i \theta_{\boldsymbol{k}}} ]$, similar to case (ii), the resulting state spontaneously breaks time reversal symmetry, resulting in a first order topological p + ip state e.g. $\boldsymbol{d}_+e^{+ i \theta_{\boldsymbol{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.
• Case (v) looks at a degenerate manifold of s-wave, $d+i\tau d$, and a mixed state with s-wave in one valley and a $d+i d$ pairing in the opposite valley, the latter of which breaks time reversal. The intravalley d-wave state is a $d+i\tau d$ superconductor—a time reversal invariant combination of d + id pairing in one valley, and d − id pairing in the other—and exhibits second order topology with gapless corner states.

Appendix. 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.

Preliminaries

The RPA equations for the screened interactions are

$$\begin{align} \widetilde{V}_{\mu \nu}(\omega, \boldsymbol{q}) = V_{\mu \nu}(\boldsymbol{q}) + V_{\mu \alpha}(\boldsymbol{q}) \Pi^{\alpha \gamma} (\omega,\boldsymbol{q}) \widetilde{V}_{\gamma \nu} (\omega, \boldsymbol{q}) \end{align} \tag{ A12 }$$

where $\Pi^{\alpha \gamma}$ is the polarization operator,

$$\begin{align} \Pi^{\mu\nu}(\omega, \boldsymbol{q})& = -i \textrm{Tr} \int J^{\,\mu} G(E+\omega, \boldsymbol{k}+\boldsymbol{q}) J^\nu G(E, \boldsymbol{k})\nonumber\\ & \quad \times \frac{dE d^{\,2}\boldsymbol{k}}{(2\pi)^3}, \nonumber \\ G(E, \boldsymbol{k})& = \frac{1}{E +\mu - \upsilon \tau_z \boldsymbol{k} \cdot \boldsymbol{\alpha} +i 0~\textrm{sgn}(E)} \end{align} \tag{ A13 }$$

where the vertices $J^{\,\mu} = \{J_1,\dots, J_{10}\} = \{\tau_0~\alpha_0, \tau_z\alpha_z,$ $\tau_0\alpha_x, \tau_0\alpha_y, \tau_x\alpha_0, \tau_y\alpha_0, \tau_x \alpha_x, \tau_x\alpha_y, \tau_y \alpha_x,\tau_y \alpha_y\}$, as per section 4 of the main text. We choose units in which υ = 1, and upon performing the frequency integral by residues, manipulations presented in [44] result in

$$\begin{align} & \Pi^{\mu\nu}(\omega,q)\nonumber\\ & \;\; = \textrm{Tr}\int \frac{d^{\,2}k}{(2\pi)^2} \sum_{s = \pm}\nonumber \\ & \quad \times \frac{J^{\,\mu} (s\tilde{\omega}+k - \tau_z(\boldsymbol{k} - s\boldsymbol{q} )\cdot \boldsymbol{\alpha})J^\nu(-k+ \tau_z\boldsymbol{k} \cdot \boldsymbol{\alpha})}{2k((s\tilde{\omega}+k)^2-(\boldsymbol{k} - s\boldsymbol{q})^2)} \nonumber \\ & \quad +\frac{J^{\,\mu} (s\tilde{\omega}+k +\tau_z(-\boldsymbol{k} + s\boldsymbol{q} )\cdot \boldsymbol{\alpha})J^\nu(k-\tau_z\boldsymbol{k} \cdot \boldsymbol{\alpha})}{2k((s\tilde{\omega}+k)^2-(\boldsymbol{k} - s\boldsymbol{q})^2)}\nonumber\\ & \quad \times \Theta(\mu-k). \end{align} \tag{ A14 }$$

The first term—the interband polarization operator, denoted $\Pi_+$—contributes when µ = 0, while the second term—the intraband polarization operator, denoted $\Pi_-$—only contributes when µ ≠ 0. Explicitly, we write

$$\begin{align} &\Pi^{\mu\nu}_+(\omega,q)\nonumber\\ & \;\; = \textrm{Tr}\int \frac{d^{\,2}k}{(2\pi)^2} \sum_{s = \pm}\nonumber\\ & \quad \times \frac{J^{\,\mu} (s\tilde{\omega}+k -\tau_z(\boldsymbol{k} - s\boldsymbol{q} )\cdot \boldsymbol{\alpha})J^\nu(k-\tau_z\boldsymbol{k} \cdot \boldsymbol{\alpha})}{2k((s\tilde{\omega}+k)^2-(\boldsymbol{k} - s\boldsymbol{q})^2)}\nonumber\\ & \quad \times \Theta(\mu-k) \end{align} \tag{ A15 }$$

and

$$\begin{align} & \Pi^{\mu\nu}_-(\omega,q)\nonumber\\ & \;\; = -i\textrm{Tr} \int \frac{d^{\,2}k}{(2\pi)^2}\frac{dE}{2\pi}\nonumber\\ & \quad \times \frac{J^{\,\mu}((\omega+E)e^{i0}+\tau_z(\boldsymbol{k}+\boldsymbol{q})\cdot \boldsymbol{\alpha})J^{\nu}(E e^{i0}+\tau_z\boldsymbol{k}\cdot \boldsymbol{\alpha})}{(((\omega+E)e^{i0})^2-(\boldsymbol{k}+\boldsymbol{q})^2)((E e^{i0})^2-k^2)}. \end{align} \tag{ A16 }$$

The polarization operator Π¹¹ is the standard charge polarization operator for graphene [60, 61, 99], which describes the screening of the density–density interaction $\alpha_0\tau_0\otimes \alpha_0\tau_0$. The function Π²² is the pseudospin polarization operator which describes the screening behavior of the pseudospin–pseudospin interaction $\alpha_z\tau_z\otimes \alpha_z\tau_z$ which was first calculated in [44] (denoted there as $\Pi^{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 [100–102], 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^{\,\mu} = (\tau_0\alpha_0,\tau_0\alpha_i)$ where $i = x,y$, and denote Π²² as $\Pi^{zz}$. Considering $J^{\,\mu} = \alpha_i \tau_0$ and $J^\nu = \alpha_j \tau_0$, i.e. the chiral current screening $\Pi^{33},\Pi^{34},\Pi^{44}$, we then write this function as $\Pi^{ij}$ and decompose into longitudinal and transverse parts:

$$\begin{align} \Pi^{ij}(\omega, \boldsymbol{q}) = \Pi_\perp(\omega, \boldsymbol{q})(q^2\delta^{ij} - q^i q^j) +\Pi_\parallel(\omega, \boldsymbol{q}) q^i q^j. \end{align} \tag{ A17 }$$

The chiral current polarization operator obeys an important identity due to electromagnetic gauge invariance. In the Dirac theory, the electromagnetic current operator $\boldsymbol{J} = \boldsymbol{\alpha}$ so that the vector potential of electromagnetism A appears in the Hamiltonian as the perturbation $e\boldsymbol{J} \cdot \boldsymbol{A}$. Gauge invariance, or equivalently the conservation of charge, can be shown to imply the Ward identity,

$$\begin{align} q_\nu\Pi^{\mu\nu}(\boldsymbol{q}, \omega) = 0. \end{align} \tag{ A18 }$$

From this we find

$$\begin{align} \omega \Pi^{0i}(\omega, \boldsymbol{q}) - q_j\Pi^{ji}(\omega, \boldsymbol{q}) = 0~\nonumber \\ \omega \Pi^{00}(\omega, \boldsymbol{q}) - q_i\Pi^{0i}(\omega, \boldsymbol{q}) = 0.\nonumber \end{align}$$

Combining these,

$$\begin{align} \omega^2~\Pi^{00}(\omega, \boldsymbol{q}) = q_i q_j \Pi^{ij}(\omega, \boldsymbol{q}) = q^4~\Pi_\perp(\omega, \boldsymbol{q}) \end{align} \tag{ A19 }$$

which gives us an expression for the longitudinal part of the current-current susceptibility in terms of the density-density response

$$\begin{align} \Pi_\parallel(\omega, \boldsymbol{q}) = \frac{\omega^2}{q^4} \Pi^{00}(\omega, \boldsymbol{q}). \end{align} \tag{ A20 }$$

Comparing our below results with the expression for Π⁰⁰ cite [44] shows our expressions satisfy this relation. We shall also derive an expression for the transverse part of the susceptibility in the dimensional regularization scheme,

$$\begin{align} q^2\Pi_{\perp}(\omega, \boldsymbol{q}) = \frac{q^2-\omega^2}{q^2}\Pi^{00}(\omega, \boldsymbol{q}) + \Pi^{zz}(\omega, \boldsymbol{q}) \end{align} \tag{ A21 }$$

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 $\Pi^{00}(\omega,\boldsymbol{q})$ obeys the exact compressibility sum rule, $\Pi^{00}(\omega = 0,\boldsymbol{q}\rightarrow 0) = -\nu_0$ where ν₀ 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 $\Pi^{ij}(\omega,\boldsymbol{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 $\Pi^{ij}(\omega,\boldsymbol{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 $\Pi^{ij}(\omega,\boldsymbol{q})$, but these UV-contributions originate from physics away from the Dirac point where the Ward identity applies. Similar to $\Pi^{00}(\omega,\boldsymbol{q})$, for small frequencies and momenta, the regularized $\Pi^{ij}(\omega,\boldsymbol{q})$ must be the same in lattice and dimensional regularization.

The pseudospin susceptibility $\Pi^{zz}(\omega,\boldsymbol{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 dependence—contributions which originate from physics near the Dirac point, and contributions which originate from physics away from the Dirac point. While only the latter affects $\Pi^{ij}(\omega,\boldsymbol{q})$, as we explained above, both types of UV-contributions appear in $\Pi^{zz}(\omega,\boldsymbol{q})$; one way to see this is to use a hard cutoff Λ in the Dirac theory, where (restoring the Dirac velocity v) one finds $\Pi^{00}(\omega = 0,\boldsymbol{q}\rightarrow 0) = -\Lambda/(2\pi v)+\nu_0$. Thus by solely modifying physics near the Dirac point—for 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 (equation 4)—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 $\Pi^{ij}(\omega,\boldsymbol{q})$ and $\Pi^{zz}(\omega,\boldsymbol{q})$ receive contributions from physics away from the Dirac point, and $\Pi^{zz}(\omega,\boldsymbol{q})$ receives an additional constant contribution from physics near the Dirac point.

Intravalley current-current polarization operator

We now turn to a calculation of $\Pi^{ij}(\omega,\boldsymbol{q})$. We begin by calculating the interband part through dimensional regularization, using the formulae

$$\begin{align} \int \frac{d^dk}{(2\pi)^d} \frac{(k^2)^a}{(k^2+\Delta)^b}& = \frac{1}{(4\pi)^{d/2}}\nonumber\\ & \quad \times \frac{\Gamma(\frac{d}{2}+a)\Gamma(b-a-\frac{d}{2})}{\Gamma(\frac{d}{2})\Gamma(b)}\nonumber\\ & \quad \times \Delta^{d/2+a-b} \end{align} \tag{ A22 }$$

$$\begin{align} \int \frac{d^dk}{(2\pi)^d} \frac{(k^2)^a}{(k^2-\Delta)^b}& = i\ \frac{(-1)^{a-b}}{(4\pi)^{d/2}}\nonumber\\ & \quad \times \frac{\Gamma(\frac{d}{2}+a)\Gamma(b-a-\frac{d}{2})}{\Gamma(\frac{d}{2})\Gamma(b)}\nonumber\\ & \quad \times \Delta^{d/2+a-b}. \end{align} \tag{ A23 }$$

Focus first on the denominator in equation (A16). Using the Schwinger–Feynman parametrization,

$$\begin{align} \frac{1}{AB} = \int_0^1\frac{dx}{(xA+(1-x)B)^2} \end{align} \tag{ A24 }$$

we write

$$\begin{align} &\frac{1}{((\omega+E)^2-(\boldsymbol{k}+\boldsymbol{q})^2)(E^2-k^2)}\nonumber\\ & \;\; = \int_0^1 dx \ \frac{1}{((k+xq)^2-x(x-1)q^2)^2} \end{align} \tag{ A25 }$$

where we now use relativistic notation $l^\mu = (E,\boldsymbol{k})$, $p^\mu = (\omega,\boldsymbol{q})$, and $l^2 = E^2-\boldsymbol{k}^2$, $p^2 = \omega^2-\boldsymbol{q}^2$. Shifting $l\rightarrow l-xp$ (bear in mind that this will affect the numerator as well), and Wick rotating $E \rightarrow iE$, the expression becomes

$$\begin{align} \int_0^1 dx \ \frac{1}{(l^2+x(x-1)p^2)^2}. \end{align} \tag{ A26 }$$

The corresponding numerator (before Wick rotation) is

$$\begin{align} & \textrm{Tr} \ J^{\,\mu} \left((1-x)\omega+E + [(1-x)\boldsymbol{q}+ \boldsymbol{k}]\cdot \boldsymbol{\alpha} \right)\nonumber\\ & \quad \times J^\nu \left( -x\omega + E +[-x\boldsymbol{q}+\boldsymbol{k}]\cdot\boldsymbol{\alpha}\right). \end{align} \tag{ A27 }$$

Now we substitute $J^{\,\mu},J^\nu = \alpha_i,\alpha_j$ with $i,j = x,y$ and perform the pseudospin trace using the identity

$$\begin{align} & \textrm{Tr} \left[ \alpha_i (A + B^\lambda\alpha_\lambda) \alpha_j (C + D^\rho \alpha_\rho) \right]\nonumber\\ & \quad = 2 AC\delta^{ij} +2\left( B^i D^j + B^j D^i - B\cdot D \delta^{ij} \right) \end{align} \tag{ A28 }$$

where $\Delta = x(x-1)p^2$. Using rotational symmetry to replace $(l^i)^2\rightarrow\frac{1}{3}l^2$, we arrive at

$$\begin{align} \Pi^{ij}_-& = 2N\int_0^1 dx \int \frac{d^dl}{(2\pi)^d}\nonumber\\ & \quad \times \frac{[x(x-1)p^2-\frac{1}{3}l^2]\delta^{ij}-2x(1-x) q^i q^j}{(k^2+\Delta)^2} \end{align} \tag{ A29 }$$

where we account for the valley and spin trace through a factor N = 4. Integrating using equation (A23), we find the imaginary part

$$\begin{align} & \textrm{Im} \ \Pi^{ij}_-\nonumber\\ & \quad = 2N\int_0^1~\frac{dx}{8\pi} \nonumber\\ & \qquad \times \left\{\left\{\frac{x(x-1)p^2}{\sqrt{x(1-x)p^2}}-\sqrt{x(1-x)p^2} \right\}\delta^{ij} \right.\nonumber\\ & \qquad\qquad \left.+\frac{2x(x-1)p^ip^j}{\sqrt{x(1-x)p^2}}\right\}\Theta(p^2). \end{align} \tag{ A30 }$$

Using $\int_0^1 dx \sqrt{x(1-x)} = \pi/8$, we arrive at

$$\begin{align} \textrm{Im} \ \Pi^{ij}_-& = -\frac{N}{16\sqrt{\omega^2-q^2}}\left\{\delta^{ij}(\omega^2-q^2) + q^i q^j \right\}\nonumber\\ & \quad \times \Theta(\omega-q). \end{align} \tag{ A31 }$$

Similarly using equation (A22), the real part is calculated from equation (A29) to be

$$\begin{align} \textrm{Re} \ \Pi^{ij}_-& = -\frac{N}{16\sqrt{q^2-\omega^2}}\left\{\delta^{ij}(\omega^2-q^2) + q^i q^j \right\} \nonumber\\ & \quad \times \Theta(q-\omega). \end{align} \tag{ A32 }$$

The functions $\Pi^{ij}_\perp(\Omega, \boldsymbol{q})$ and $\Pi^{ij}_\parallel(\Omega, \boldsymbol{q})$ are respectively the transverse and longitudinal components of the polarization operator. From above, we see the their interband parts are

$$\begin{align} \Pi^{ij}_{-,\perp}& = -\frac{N}{16 q^2} \sqrt{q^2-\omega^2} \ \Theta(q-\omega)\nonumber\\[5pt] & \quad + i\frac{N}{16 q^2}\sqrt{\omega^2-q^2}\Theta(\omega-q) \end{align} \tag{ A33 }$$

and

$$\begin{align} \ \Pi^{ij}_{-,\parallel}& = -\frac{N}{16q^2}\frac{\omega^2}{\sqrt{q^2-\omega^2}} \ \Theta(q-\omega)\nonumber\\[5pt] & \quad -i\frac{N}{16q^2}\frac{\omega^2}{\sqrt{\omega^2-q^2}} \ \Theta(\omega-q). \end{align} \tag{ A34 }$$

Calculating the density response $J^{\,\mu} = J^\nu = \alpha_0$ through the formalism above, one may see that the above expression for $\Pi^{ij}_\parallel$ satisfies the Ward identity (A20).

Let us now perform the intraband calculation. The numerator in equation (A15) is evaluated using equation (A28), so that

$$\begin{align} \Pi^{ij}_+(\omega,q) = 2N\int \frac{d^{\,2}k}{(2\pi)^2} \sum_{s = \pm} \frac{\{k(s\widetilde{\omega}+k) - \boldsymbol{k} \cdot (\boldsymbol{k}- s \boldsymbol{q})\}\delta^{ij} + 2k^ik^j - s(k^iq^j + k^jq^i) }{2k((s\tilde{\omega}+k)^2-(\boldsymbol{k} - s\boldsymbol{q})^2)} \Theta(\mu-k). \end{align} \tag{ A35 }$$

The denominator equals $\widetilde{\omega}^2 + 2s\widetilde{\omega}k -q^2+2~s \boldsymbol{k} \cdot \boldsymbol{q}$. Shifting the angle of integration $\theta_{\boldsymbol{k} } \rightarrow \theta_{\boldsymbol{k} } + \theta_{\boldsymbol{q} }$ where $\theta_{\boldsymbol{q} }$ is the angle of the vector q , and denoting the angle of integration hereafter as simply θ, we have $\boldsymbol{k} \cdot \boldsymbol{q} \rightarrow kq \cos \theta$. Using the evenness of the denominator in θ, the numerator may be simplified to arrive at

$$\begin{align} \Pi^{ij}_+(\omega,q) = 2N\int \frac{dk d\theta}{(2\pi)^2} \sum_{s = \pm}\frac{\{s\widetilde{\omega}k + skq \cos \theta\}\delta^{ij} + 2k^2\cos\theta \delta^{ij} + 2k^2\cos 2\theta \hat{q}^{\,i}\hat{q}^{\,j} - 2skq\cos\theta\hat{q}^{\,i}\hat{q}^{\,j} }{\widetilde{\omega}^2 + 2s\widetilde{\omega}k -q^2+2~s kq \cos\theta} \Theta(\mu-k) \end{align} \tag{ A36 }$$

where $\hat{q}^{\,i} = q^i/q$. To proceed, we make use of the integral

$$\begin{align} I(a + ib0)& = \int_0^{2\pi} \frac{1}{a+ib0+\cos\theta} \frac{d\theta}{2\pi } \nonumber\\[16pt] & = \frac{\textrm{sgn}(a)}{\sqrt{a^2-1}} \Theta(|a|-1)-i \frac{\textrm{sgn}(b)}{\sqrt{1-a^2}}\nonumber\\[15pt] & \quad \times \Theta(1-|a|) \end{align} \tag{ A37 }$$

and the identities

$$\begin{align} \int_0^{2\pi} \frac{\cos\theta}{a+\cos\theta} \frac{d\theta}{2\pi} & = 1- aI(a ) \end{align} \tag{ A38 }$$

$$\begin{align} \int_0^{2\pi} \frac{\cos 2\theta}{a+\cos\theta} \frac{d\theta}{2\pi} & = -2a + (2a^2-1)I(a ). \end{align} \tag{ A39 }$$

Evaluating (A36), we define $\alpha = \frac{\widetilde{\omega}^2 + 2s\widetilde{\omega}k -q^2}{2kq} = \frac{{\omega}^2 + 2s{\omega}k -q^2}{2kq} + i 0\frac{{\omega}^2 + s{\omega}k}{kq}$ and denote $\tilde{I}(\alpha) = 2kqI(\alpha)$. Then

$$\begin{align} \Pi^{ij}_+(\omega,q) = 2N\int \frac{dk }{2\pi} \sum_{s = \pm} & \Bigg( \left(\frac{1}{2}+ \frac{k\alpha}{q}+\{s\widetilde\omega k-kq \alpha+2k^2(1-\alpha_2)\}\tilde I(\alpha)\right)\delta^{ij}\;\nonumber\\ & \quad + \Bigg(-1+ \frac{2k\alpha}{q} + \{2k^2+2kq\alpha-4k^2(1-\alpha_2)\}\tilde I(\alpha) \Bigg)\hat{q}^{\,i}\hat{q}^{\,j} \Bigg) \Theta(\mu-k). \end{align} \tag{ A40 }$$

Making the substitution $2p = 2k+s\widetilde\omega$, we get

$$\begin{align} \Pi^{ij}_+(\omega,q) = \int \frac{dk }{2\pi} \sum_{s = \pm}\left( \left(\frac{\widetilde{\omega}^2}{2q^2}+ \frac{1}{q^2}(q^2-\widetilde{\omega}^2)\tilde I(\alpha)\right)\left[\delta^{ij}-\hat{q}^{\,i}\hat{q}^{\,j}\right] - \left(\frac{\widetilde{\omega}^2}{2q^2}-\frac{\widetilde{\omega}^2}{2q^2}(4p^2-q^2) \tilde I(\alpha) \right)\hat{q}^{\,i}\hat{q}^{\,j} \right) \Theta(\mu-k). \end{align} \tag{ A41 }$$

We arrive at

$$\begin{align} \Pi^{ij}_+(\omega, \boldsymbol{q})& = \Pi_{+,\perp}(\omega, \boldsymbol{q})(q^2\delta^{ij} - q^i q^j)\nonumber\\ & \quad +\Pi_{+,\parallel}(\omega, \boldsymbol{q}) q^i q^j \end{align} \tag{ A42 }$$

with

$$\begin{align} \Pi_{+,\perp}& = 2N\int \frac{dk }{2\pi} \sum_{s = \pm}\left( \frac{\widetilde{\omega}^2}{2q^4}+ \frac{1}{q^4}(q^2-\widetilde{\omega}^2)\tilde I(\alpha) \right)\nonumber\\ & \quad \times \Theta(\mu-k) \end{align} \tag{ A43 }$$

and

$$\begin{align} \Pi_{+,\parallel}& = N \frac{\widetilde{\omega}^2}{q^4} \int \frac{dk }{2\pi} \sum_{s = \pm}\left(-1 +(4p^2-q^2)\tilde I(\alpha) \right)\nonumber\\ & \quad \times \Theta(\mu-k). \end{align} \tag{ A44 }$$

As stated earlier, the latter function is related to the density response through the Ward identity. Direct comparison with the appendix of [44] shows that the Ward identity is satisfied by the above expression. One also sees that the inter and intraband polarization operators separately satisfy equation (A21). In summary, one finds

$$\begin{align*} \Pi^{ij}(\omega, \boldsymbol{q}) = \Pi_\perp(\omega, \boldsymbol{q})(q^2\delta^{ij} - q^i q^j) +\Pi_\parallel(\omega, \boldsymbol{q}) q^i q^j \end{align*}$$

with

$$\begin{align} q^2\Pi_{\perp}(\omega, \boldsymbol{q}) = \frac{q^2-\omega^2}{q^2}\Pi^{00}(\omega, \boldsymbol{q}) + \Pi^{zz}(\omega, \boldsymbol{q}) \end{align} \tag{ A45 }$$

$$\begin{align} \Pi_\parallel(\omega, \boldsymbol{q}) = \frac{\omega^2}{q^4} \Pi^{00}(\omega, \boldsymbol{q}) \end{align} \tag{ A46 }$$

as claimed in the previous section.

Intervalley polarization operator

Taking $J^{\,\mu} = \tau_\pm \alpha_i$ and $J^\nu = \tau_\mp \alpha_j$ with $\tau_\pm = \frac{1}{2}(\tau_x\pm i\tau_y)$ and performing the trace changes the sign of the terms containing factors of τ_z . If we denote this intervalley polarization operator $\Pi_{II}^{ij}$ then one finds the simple result

$$\begin{align} \Pi_{II, \perp}^{ij} & = \Pi_{\parallel}^{ij} \end{align} \tag{ A47 }$$

$$\begin{align} \Pi_{II, \parallel}^{ij} & = \Pi_{\perp}^{ij}. \end{align} \tag{ A48 }$$

In other words, the transverse and longitudinal response are reversed between intra and intervalley. The same trace relations show that $\Pi^{22} = \Pi^{55} = \Pi^{66} = \Pi^{zz}$, i.e. the screening of the operators $\tau_\pm$ is the same as that of $\tau_z\alpha_z$.

Solution to the RPA equations in the onsite limit

We have the interactions

$$\begin{align} & g_1\alpha_0\otimes \alpha_0+ g_2\alpha_z\tau_z\otimes \alpha_z\tau_z + g_4(\tau_+\otimes \tau_- + \tau_-\otimes \tau_+) \nonumber \\[5pt] & g_3\left( \alpha_x\otimes \alpha_x + \alpha_y \otimes \alpha_y \right)+g_5\left( \alpha_x\otimes \alpha_x + \alpha_y \otimes \alpha_y \right)\nonumber\\[5pt] & \quad \times \left(\tau_+ \otimes \tau_-+\tau_- \otimes \tau_+ \right). \end{align} \tag{ A49 }$$

We consider the simple case of only onsite interactions, so only $U_{AAAA}\neq 0$; in this limit one finds $g_3 = g_4 = 0$ for the honeycomb lattice and $g_2 = g_5 = 0$ for kagome lattices. In each case the RPA equations decouple, so that each g_i is screened separately. The result is the screened interactions:

$$\begin{align} & \frac{g_1}{1-g_1\Pi^{00}}\alpha_0\otimes \alpha_0+ \frac{g_2}{1-g_2\Pi^{zz}}\alpha_z\tau_z\otimes \alpha_z\tau_z \nonumber \\ & \quad+ \left( \frac{g_5 + g_5^2 q^2~\Pi_\parallel }{1-g_5(\Pi_\parallel+\Pi_\perp)q^2 + g^2_5\Pi_\parallel\Pi_\perp q^4}\left( \alpha_x\otimes \alpha_x + \alpha_y \otimes \alpha_y \right) \right.\nonumber \\ & \qquad\quad \left. - \frac{g_5^2(\Pi_\perp-\Pi_\parallel)}{1-g_5(\Pi_\parallel+\Pi_\perp)q^2 + g^2_5\Pi_\parallel\Pi_\perp q^4}(\boldsymbol{q} \cdot \boldsymbol{\alpha}) \otimes (\boldsymbol{q} \cdot \boldsymbol{\alpha}) \right) \left(\tau_+ \otimes \tau_-+\tau_- \otimes \tau_+ \right) \end{align} \tag{ A50 }$$

for the honeycomb lattice and

$$\begin{array}{l} \frac{{{g_1}}}{{1 - {g_1}{\Pi ^{00}}}}{\alpha _0} \otimes {\alpha _0} + + \frac{{{g_4}}}{{1 - {g_4}{\Pi ^{zz}}}}({\tau _ + } \otimes {\tau _ - } + {\tau _ - } \otimes {\tau _ + })\\ + \frac{{{g_3} + g_3^2{q^2}{\Pi _\parallel }}}{{1 - {g_3}({\Pi _\parallel } + {\Pi _ \bot }){q^2} + g_3^2{\Pi _\parallel }{\Pi _ \bot }{q^4}}}\left( {{\alpha _x} \otimes {\alpha _x} + {\alpha _y} \otimes {\alpha _y}} \right)\\ + \frac{{g_3^2({\Pi _ \bot } - {\Pi _\parallel })}}{{1 - {g_3}({\Pi _\parallel } + {\Pi _ \bot }){q^2} + g_3^2{\Pi _\parallel }{\Pi _ \bot }{q^4}}}({\boldsymbol{q}} \cdot {\boldsymbol{\alpha }}) \otimes ({\boldsymbol{q}} \cdot {\boldsymbol{\alpha }})\\ + \left( {\frac{{{g_5} + g_5^2{q^2}{\Pi _\parallel }}}{{1 - {g_5}({\Pi _\parallel } + {\Pi _ \bot }){q^2} + g_5^2{\Pi _\parallel }{\Pi _ \bot }{q^4}}}\left( {{\alpha _x} \otimes {\alpha _x} + {\alpha _y} \otimes {\alpha _y}} \right)} \right.\\ \left. { - \frac{{g_5^2({\Pi _ \bot } - {\Pi _\parallel })}}{{1 - {g_5}({\Pi _\parallel } + {\Pi _ \bot }){q^2} + g_5^2{\Pi _\parallel }{\Pi _ \bot }{q^4}}}({\boldsymbol{q}} \cdot {\boldsymbol{\alpha }}) \otimes ({\boldsymbol{q}} \cdot {\boldsymbol{\alpha }})} \right)\left( {{\tau _ + } \otimes {\tau _ - } + {\tau _ - } \otimes {\tau _ + }} \right) \end{array} \tag{ A51 }$$

for the kagome lattice.

Cooper channel scattering amplitudes

On-site limit for honeycomb systems

For honeycomb systems, the $\ell$-wave scattering amplitude is

$$\begin{align} \Gamma^\ell_{\tau_1\tau_2\tau_3\tau_4}(\omega, k,p) & = \frac{1}{4}\int \frac{d\theta }{2\pi } e^{-i\ell\theta}\left( (1+e^{i\tau_1\theta})^2\frac{g_0}{1-g_0\Pi^{00}(\omega, \boldsymbol{q})} + (1-e^{i\tau_1\theta})^2\frac{g_2}{1-g_2\Pi^{zz}(\omega, \boldsymbol{q})} \right)\delta_{\tau_1,\tau_2}\delta_{\tau_3,\tau_4}\delta_{\tau_1,\tau_3} \nonumber \\ & \quad+2e^{-i\ell\theta}\left( (1+\cos\theta)\frac{g_0}{1-g_0\Pi^{00}(\omega, \boldsymbol{q})} + (1-\cos\theta)\frac{g_2}{1-g_2\Pi^{zz}(\omega, \boldsymbol{q})} \right)\delta_{\tau_1,\tau_2}\delta_{\tau_3,\tau_4}\delta_{\tau_1,-\tau_3} \nonumber \\ & \quad + \ 2e^{-i\ell\theta}\left( \frac{g_5 + g_5^2 q^2~\Pi_\parallel }{1-g_5(\Pi_\parallel+\Pi_\perp)q^2 + g^2_5\Pi_\parallel\Pi_\perp q^4}\right)\delta_{\tau_1,-\tau_2}\delta_{\tau_3,-\tau_4}\delta_{\tau_1,-\tau_3} \end{align} \tag{ A52 }$$

where $\boldsymbol{q} = \boldsymbol{k} -\boldsymbol{p}$, $\theta = \theta_{\boldsymbol{k}}-\theta_{\boldsymbol{p}}$.

On-site limit for kagome systems

For the kagome lattice we have

$$\begin{align} \Gamma^\ell_{\tau_1\tau_2\tau_3\tau_4}(\omega, k,p)& = \frac{1}{4}\int \frac{d\theta }{2\pi } e^{-i\ell\theta}\left( (1+e^{i\tau_1\theta})^2\frac{g_0}{1-g_0\Pi^{00}(\omega, \boldsymbol{q})} - 2e^{i\tau\theta}\frac{g_3 + g_3^2 q^2~\Pi_\parallel }{1-g_3(\Pi_\parallel+\Pi_\perp)q^2 + g^2_3\Pi_\parallel\Pi_\perp q^4} \right. \nonumber \\ & \qquad\qquad\qquad\quad \left. +\ (k-p)^2(1+e^{i\tau_1\theta})^2\frac{g_3^2 (\Pi_\perp-\Pi_\parallel) }{1-g_3(\Pi_\parallel+\Pi_\perp)q^2 + g^2_3\Pi_\parallel\Pi_\perp q^4}\right)\delta_{\tau_1,\tau_2}\delta_{\tau_3,\tau_4}\delta_{\tau_1,\tau_3} \nonumber \\ & \quad +\ e^{-i\ell\theta}\left( (1+\cos\theta)\frac{2g_0}{1-g_0\Pi^{00}(\omega, \boldsymbol{q})} +\frac{2(g_3 + g_3^2 q^2~\Pi_\parallel ) }{1-g_3(\Pi_\parallel+\Pi_\perp)q^2 + g^2_3\Pi_\parallel\Pi_\perp q^4} \right. \nonumber \\ & \qquad\qquad\quad \left.-\ (k-p)^2(1+\cos\theta)\frac{2g_3^2 (\Pi_\perp-\Pi_\parallel) }{1-g_3(\Pi_\parallel+\Pi_\perp)q^2 + g^2_3\Pi_\parallel\Pi_\perp q^4} \right)\delta_{\tau_1,\tau_2}\delta_{\tau_3,\tau_4}\delta_{\tau_1,-\tau_3} \nonumber \\ & \quad+ \ e^{-i\ell\theta}\left( (1-\cos\theta)\frac{2g_4}{1-g_4\Pi^{zz}(\omega, \boldsymbol{q})} +\frac{2(g_5 + g_5^2 q^2~\Pi_\parallel) }{1-g_5(\Pi_\parallel+\Pi_\perp)q^2 + g^2_5\Pi_\parallel\Pi_\perp q^4}\right) \delta_{\tau_1,-\tau_2}\delta_{\tau_3,-\tau_4}\delta_{\tau_1,-\tau_3}. \end{align} \tag{ A53 }$$

Eliashberg equations

Neglecting self-energy corrections, the linearized Gor'kov–Eliashberg frequency dependent gap equation is given

$$\begin{align} \Delta^{\ell}(i\omega_n, k) = -T\sum_{m} \sum_p K^\ell(i\omega_n, k ; i\omega_m, p) \Delta^{\ell}(i\omega_m,p) \end{align} \tag{ A54 }$$

with the Eliashberg kernel,

$$\begin{align} & K^\ell(i\omega_n, k ; i\omega_m, p)\nonumber\\ & \;\;= \frac{1}{\omega_m^2 + \varepsilon_p^2} \int \frac{d\theta}{2\pi}\ e^{i\ell \theta} \ \lambda(i\omega_n - i\omega_m, k, p, \theta) \end{align} \tag{ A55 }$$

where $\omega_n = \pi(2n+1)T$ are the fermionic Matsubara frequencies [103–105]. The frequency dependent coupling $\lambda(\omega, q)$ is defined analogously to the static coupling in the main text: $\lambda(\omega,q) = \Gamma_{\tau\tau\tau\tau}(\omega,q)$ for intravalley pairing or $\lambda(\omega,q) = \Gamma_{+-+-}(\omega,q) \pm \Gamma_{+-+}(\omega,q)$ for intervalley pairing. The functions Δ, K and Γ are tensors in spin and valley space—spin and valley indices are left implicit. The critical temperature can found by finding the value of T for which the above linear equation for Δ has an eigenvalue of $-1/T$. The zero temperature gap can also be found from solving

$$\begin{align} \Delta^\ell(\omega^{^{\prime}}) = -\nu_0\int \frac{\Delta^\ell(\omega)\lambda^\ell(\omega-\omega^{^{\prime}})}{\sqrt{\omega^2 - |\Delta^\ell(\omega)|^2}} \frac{d\omega}{4\pi} \end{align} \tag{ A56 }$$

where $\lambda^\ell(\omega-\omega^{^{\prime}})$ is the $\ell$-wave coupling evaluated onshell, $\varepsilon_k = \omega, \varepsilon_p = \omega^{^{\prime}}$. These equations can be derived through general assumptions about the analytic behavior of the scattering amplitude Γ (see e.g. [103]). The requirements of unitarity and causality imply that as a function of frequency the amplitude must be analytic in the upper half plane, and that the amplitude has branch cuts corresponding to the threshold of particle-hole production. While Γ possesses singularities as a function of $\omega^{^{\prime}}$, these occur either for $\omega^{^{\prime}}\gg v_F|\boldsymbol{k}-\boldsymbol{p}|$, which provides a negligible contribution, or close to the edge of the particle hole continuum $\omega^{^{\prime}}\approx v_F|\boldsymbol{k}-\boldsymbol{p}|$, which is significant for a range of scattering angles $\theta_{\boldsymbol{k}}-\theta_{\boldsymbol{p}}\approx 0$ which only becomes significant when $p\approx k_F$.

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—$\lambda^{^{\ell = 1}}_{intra} = \Gamma^{\ell = 1}_{\tau\tau\tau\tau}(|k|-|p|,|\boldsymbol{k}-\boldsymbol{p}|)$ as well as $\lambda^{\ell = 1}_{inter} = \Gamma_{+-+-}(|k|-|p|,|\boldsymbol{k}-\boldsymbol{p}|) - \Gamma_{+-+}(|k|-|p|,|\boldsymbol{k}-\boldsymbol{p}|)$—for the kagome and honeycomb lattices 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 [106], 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 [44].

Zoom In Zoom Out Reset image size

Appendix. Real space form of the $p+i\tau p$ gap function

In this appendix we derive the real space form of the $p+i\tau p$ gap function used for exact diagonalization in section 4. Starting with the momentum space mean field Hamiltonian, we have

$$\begin{align} H_\textrm{BdG}& = \sum_{\boldsymbol{k}}\varepsilon_{\boldsymbol{k}} f^{\,\dagger}_{\boldsymbol{k},\tau,s}\;f_{\boldsymbol{k},\tau,s} + \frac{1}{2}\sum_{k \lt k_F,\tau,s,s^{^{\prime}}}\nonumber\\ & \quad \times {\Delta_k e^{i\tau \phi} e^{-i\tau \theta_{\boldsymbol{k}}} (i s_y d^\mu s_\mu)_{ss^{^{\prime}}}f^{\,\dagger}_{\boldsymbol{k},\tau,s}\; f^{\,\dagger}_{-\boldsymbol{k},\tau,s^{^{\prime}}}} +\textrm{h.c.} \end{align} \tag{ A57 }$$

where $f^{\,\dagger}_{\boldsymbol{k}, \tau, s}$ creates an electron in the upper band, $\Delta_k$ is an overall factor depending only on the magnitude k, φ is the phase of the pair density wave, and $(d^{\,x},d^{\,y},d^{\,z})$ is a constant 3D vector with unit length. The normal dispersioon $\varepsilon_{\boldsymbol{k}}$ is spin independent, which allows us to perform a spin rotation so that $(d^{\,x},d^{\,y},d^{\,z}) = -i(0,1,0)$, and the Hamiltonian decouples into two spin-diagonal terms,

$$\begin{align} H_{BdG} = H_\uparrow + H_\downarrow &= \sum_{\boldsymbol{k},s}\varepsilon_{\boldsymbol{k}} f^{\,\dagger}_{\boldsymbol{k},\tau,s}f_{\boldsymbol{k},\tau,s} + \frac{1}{2}\sum_{k \lt k_F,\tau,s,s^{^{\prime}}}\nonumber\\ & \quad \times {\Delta_k e^{i\tau\phi} e^{-i\tau \theta_{\boldsymbol{k}}} f^{\,\dagger}_{\boldsymbol{k},\tau,s} f^{\,\dagger}_{-\boldsymbol{k},\tau,s}} +\textrm{h.c.} \end{align} \tag{ A58 }$$

Since the Hamiltonian is a sum of independent spin blocks, we discard the spin index henceforth. We now convert this expression to real space, in which the Hamiltonian becomes

$$\begin{align} H = \sum_{\boldsymbol{r},\boldsymbol{r}^{^{\prime}}}\mathcal{H}(\boldsymbol{r},\boldsymbol{r}^{^{\prime}}) c^\dagger_{\boldsymbol{r}} c^\dagger_{\boldsymbol{r}^{^{\prime}}} + \frac{1}{2} \Delta(\boldsymbol{r},\boldsymbol{r}^{^{\prime}}) c^\dagger_{\boldsymbol{r}} c^\dagger_{\boldsymbol{r}^{^{\prime}}} + \textrm{h.c.} \end{align} \tag{ A59 }$$

where $c^\dagger_{\boldsymbol{r},s}$ is the full electron creation operator (as compared to $f^{\,\dagger}_{\boldsymbol{k}}$ which creates an electron with momentum k in the upper band). We assume that the only hopping term present connects nearest neighbors, i.e. that $\mathcal{H}(\boldsymbol{r},\boldsymbol{r}^{^{\prime}}) = -t$ if r and $\boldsymbol{r}^{^{\prime}}$ are nearest neighbors and zero else. It remains then to calculate the function $\Delta(\boldsymbol{r},\boldsymbol{r}^{^{\prime}})$. In momentum space, the pairing term for each spin block is given by

$$\begin{align} H_\Delta = \frac{1}{2}\sum_{k \lt k_F,\tau}{\Delta_k e^{i\tau \phi} e^{-i\tau \theta_{\boldsymbol{k}}} f^{\,\dagger}_{\boldsymbol{k},\tau} f^{\,\dagger}_{-\boldsymbol{k},\tau}}. \end{align} \tag{ A60 }$$

We relate the creation operator for the upper band to the full electron creation operator by

$$\begin{align} f^{\,\dagger}_{\boldsymbol{k},\tau,s} = \sum_{\boldsymbol{r}}{u_{\boldsymbol{k},\tau}(\boldsymbol{r}) c^\dagger_{\boldsymbol{r},s}} \end{align} \tag{ A61 }$$

where $u_{\boldsymbol{k},\tau}(\boldsymbol{r})$ is the wavefunction for the upper Dirac band. Substituting this relation into the pairing term,

$$\begin{align} H_\Delta = \frac{1}{2} \sum_{k,\tau,\boldsymbol{r},\boldsymbol{r}^{^{\prime}}}{\Delta_k e^{-i\tau\phi}e^{-i\tau\theta_{\boldsymbol{k}}} u_{\boldsymbol{k},\tau}(\boldsymbol{r}) u_{-\boldsymbol{k},\tau}(\boldsymbol{r}^{^{\prime}}) c^\dagger_{\boldsymbol{r},s} c^\dagger_{\boldsymbol{r}^{^{\prime}},s}} \end{align} \tag{ A62 }$$

from which it follows immediately that

$$\begin{align} \Delta(\boldsymbol{r},\boldsymbol{r}^{^{\prime}}) = \sum_{k,\tau}{\Delta_k e^{-i\tau\phi} e^{-i\tau\theta_{\boldsymbol{k}}} u_{\boldsymbol{k},\tau}(\boldsymbol{r}) u_{-\boldsymbol{k},\tau}(\boldsymbol{r}^{^{\prime}})}. \end{align} \tag{ A63 }$$

The upper band eigenstates of the Dirac Hamiltonian $\mathcal{H} = v\tau \boldsymbol{k}\cdot\boldsymbol{\alpha}$ are given in the coordinate representation in terms of the wavefunctions $e^{i\tau\boldsymbol{K}_1\cdot\boldsymbol{r}} \varphi_{\tau \alpha}(\boldsymbol{r})$ by

$$\begin{align} u_{\boldsymbol{k},+}(\boldsymbol{r}) & = \frac{1}{\sqrt{2}} \left(\varphi_{++}(\boldsymbol{r}) + e^{i\theta_{\boldsymbol{k}}} \varphi_{+-}(\boldsymbol{r})\right) e^{i(\boldsymbol{k}+\boldsymbol{K})\cdot\boldsymbol{r}}, \nonumber \\ u_{\boldsymbol{k},-}(\boldsymbol{r})& = u^*_{-\boldsymbol{k},+}(\boldsymbol{r}) = \frac{1}{\sqrt{2}} \left(\varphi^*_{++}(\boldsymbol{r}) - e^{-i\theta_{\boldsymbol{k}}} \varphi^*_{+-}(\boldsymbol{r})\right)\nonumber \\ & \qquad\qquad\qquad \times e^{i(\boldsymbol{k}-\boldsymbol{K})\cdot\boldsymbol{r}} \end{align} \tag{ A64 }$$

where $\varphi_{\tau\alpha}(\boldsymbol{r})$ are periodic functions under lattice translations. Thus

$$\begin{align} \Delta(\boldsymbol{r},\boldsymbol{r}^{^{\prime}})& = \sum_{\boldsymbol{k}}{}\Delta_k\Bigg\{ \frac{1}{2}e^{i\{\boldsymbol{K}\cdot(\boldsymbol{r}+\boldsymbol{r}^{^{\prime}}) + \boldsymbol{k}\cdot(\boldsymbol{r}-\boldsymbol{r}^{^{\prime}}) + \phi -\theta_{\boldsymbol{k}}\}} \left[ (\varphi_{++}(\boldsymbol{r}) + e^{i\theta_{\boldsymbol{k}}}\varphi_{+-}(\boldsymbol{r})) (\varphi_{++}(\boldsymbol{r}^{^{\prime}}) - e^{i\theta_{\boldsymbol{k}}} \varphi_{+-}(\boldsymbol{r}^{^{\prime}}))\right] \nonumber \\ & \qquad \qquad \quad \! + \! \frac{1}{2}e^{i\{-\boldsymbol{K}\cdot(\boldsymbol{r}+\boldsymbol{r}^{^{\prime}})+\boldsymbol{k}\cdot(\boldsymbol{r}-\boldsymbol{r}^{^{\prime}}) -\phi +\theta_{\boldsymbol{k}}\}} \left[ (\varphi^*_{++}(\boldsymbol{r}) - e^{-i\theta_{\boldsymbol{k}}} \varphi^*_{+-}(\boldsymbol{r})) (\varphi^*_{++}(\boldsymbol{r}^{^{\prime}}) + e^{-i\theta_{\boldsymbol{k}}} \varphi^*_{+-}(\boldsymbol{r}^{^{\prime}}))\right]\! \Bigg\}. \end{align} \tag{ A65 }$$

Performing the summation over k and introducing the functions

$$\begin{align} \int{\Delta_k e^{i\boldsymbol{k}\cdot\boldsymbol{r}} \frac{d^{\,2}\boldsymbol{k}}{(2\pi)^2}} = \frac{1}{2\pi}\int{\Delta_k J_0(kr) kdk} = \Phi_0(r), \nonumber \\ i^\mp\int{\Delta_k e^{i\left[\boldsymbol{k}\cdot\boldsymbol{r} \pm \theta\right]} \frac{d^{\,2}\boldsymbol{k}}{(2\pi)^2}} = \frac{1}{2\pi}\int{\Delta_k J_1(kr) kdk} = \Phi_1(r) \end{align} \tag{ A66 }$$

we find

$$\begin{align} \Delta(\boldsymbol{r},\boldsymbol{r}^{^{\prime}}) & = \frac{1}{2}e^{i\{\boldsymbol{K}\cdot(\boldsymbol{r}+\boldsymbol{r}^{^{\prime}}) + \phi \}} \left[ (\varphi_{++}(\boldsymbol{r}) \varphi_{++}(\boldsymbol{r}^{^{\prime}}) - \varphi_{+-}(\boldsymbol{r}) \varphi_{+-}(\boldsymbol{r}^{^{\prime}})) i\Phi_1(r) + (\varphi_{+-}(\boldsymbol{r}) \varphi_{++}(\boldsymbol{r}^{^{\prime}}) \right.\nonumber \\ & \qquad \qquad \qquad \qquad \left. - \varphi_{++}(\boldsymbol{r}) \varphi_{+-}(\boldsymbol{r}^{^{\prime}}))\Phi_0(r)\right] \nonumber \\ & \quad + \ \frac{1}{2} e^{-i\{\boldsymbol{K}\cdot(\boldsymbol{r}+\boldsymbol{r}^{^{\prime}})+ \phi\}} \left[-(\varphi_{++}(\boldsymbol{r}) \varphi^*_{++}(\boldsymbol{r}^{^{\prime}}) - \varphi^*_{+-}(\boldsymbol{r}) \varphi^*_{+-}(\boldsymbol{r}^{^{\prime}})) i\Phi_1(r) + (\varphi^*_{++}(\boldsymbol{r}) \varphi^*_{+-}(\boldsymbol{r}^{^{\prime}}) \right.\nonumber \\ & \qquad\qquad\qquad\qquad\qquad \left. - \varphi^*_{+-}(\boldsymbol{r}) \varphi^*_{++}(\boldsymbol{r}^{^{\prime}}) ) \Phi_0(r) \right] \nonumber \\ & = \textrm{Re}\{ e^{i\{\boldsymbol{K}\cdot(\boldsymbol{r}+\boldsymbol{r}^{^{\prime}})+ \phi+\frac{\pi}{2}\}} \left[ (\varphi_{++}(\boldsymbol{r}) \varphi_{+-}(\boldsymbol{r}^{^{\prime}}) - \varphi_{+-}(\boldsymbol{r}) \varphi_{+-}(\boldsymbol{r}^{^{\prime}})) \right.\nonumber \\ & \quad\;\; \quad\left. \times |\Phi_1(r)| + (\varphi_{++}(\boldsymbol{r})\varphi_{+-}(\boldsymbol{r}^{^{\prime}}) - \varphi_{+-}(\boldsymbol{r}) \varphi_{++}(\boldsymbol{r}^{^{\prime}}))|\Phi_0(r)|\right]\}. \end{align} \tag{ A67 }$$

The function $|\Phi_0(r)|$ is maximum for r = 0 and goes to zero over length scales $r \sim \pi/k_F$ while $|\Phi_1(r)|$ is small for small r and increases to a maximum at $r \sim\pi/k_F$. To simplify $\Delta(\boldsymbol{r},\boldsymbol{r}^{^{\prime}})$ for the purposes of exact diagonalization, we will neglect the pairing correlations at large separations, discarding the term containing Φ₁ and set $\Phi_0 \rightarrow \Delta^{^{\prime}}_{}$ with $\Delta^{^{\prime}}_{}$ being some average value. With these simplifications we finally arrive at the real space Hamiltonian

$$\begin{align} H &= \sum_{\langle \boldsymbol{r},\boldsymbol{r}^{^{\prime}}\rangle}{ -t c^\dagger_{\boldsymbol{r}} c_{\boldsymbol{r}^{^{\prime}}} + \frac{1}{2}\Delta(\boldsymbol{r},\boldsymbol{r}^{^{\prime}}) c^\dagger_{\boldsymbol{r}} c^\dagger_{\boldsymbol{r}^{^{\prime}}} + \textrm{h.c.} },\nonumber \\ \Delta(\boldsymbol{r},\boldsymbol{r}^{^{\prime}})& = \Delta^{^{\prime}}_{} \textrm{Re} \{e^{i\{\boldsymbol{K}\cdot(\boldsymbol{r}+\boldsymbol{r}^{^{\prime}}) + \phi+\frac{\pi}{2}\}} \nonumber \\ & \quad \times \left[\varphi_{++}(\boldsymbol{r}) \varphi_{+-}(\boldsymbol{r}^{^{\prime}}) -\varphi_{+-}(\boldsymbol{r})\varphi_{++}(\boldsymbol{r}^{^{\prime}})\right] \}. \end{align} \tag{ A68 }$$

Note that $\Delta(\boldsymbol{r},\boldsymbol{r}^{^{\prime}}) = -\Delta(\boldsymbol{r}^{^{\prime}},\boldsymbol{r}) = \Delta^*(\boldsymbol{r},\boldsymbol{r}^{^{\prime}})$, and with t real, the Hamiltonian is invariant under complex conjugation, conforming with our claim in the main text that this superconducting phase lies in class BDI.

Below we present $\Delta(\boldsymbol{r},\boldsymbol{r}^{^{\prime}})$ explicitly for the honeycomb and Kagome lattices. The procedure is simply to calculate the eigenfunctions $\varphi_{\tau \alpha}(\boldsymbol{r})$ of the normal state Hamiltonians at the Dirac points, and insert into equation (A68).

Honeycomb lattice

For the honeycomb lattice, we have

$$\begin{align} \varphi_{+ +}(\boldsymbol{r}) = \begin{cases} 1 \ \ \ \boldsymbol{r} \in A\\ 0 \ \ \ \boldsymbol{r} \in B \end{cases} \ \ \ \varphi_{+-}(\boldsymbol{r}) = \begin{cases} 0 \ \ \ \boldsymbol{r} \in A\\ 1 \ \ \ \boldsymbol{r} \in B. \end{cases} \end{align} \tag{ A69 }$$

For $\boldsymbol{r}\in A,\boldsymbol{r}^{^{\prime}}\in B$ we therefore get

$$\begin{align} \Delta(\boldsymbol{r},\boldsymbol{r}^{^{\prime}})& = \textrm{Re} \{e^{i\{\boldsymbol{K}\cdot(\boldsymbol{r}+\boldsymbol{r}^{^{\prime}}) + \phi+\frac{\pi}{2}\}}\}\nonumber \\ & = \sin(\boldsymbol{K}\cdot(\boldsymbol{r}+\boldsymbol{r}^{^{\prime}}) + \phi) \end{align} \tag{ A70 }$$

and hence

$$\begin{align} H& = \sum_{\langle \boldsymbol{r},\boldsymbol{r}^{^{\prime}}\rangle}-t c^\dagger_{\boldsymbol{r}} c^\dagger_{\boldsymbol{r}^{^{\prime}}} + \Delta^{^{\prime}}_{}\nonumber \\ & \quad \times \left[\sin(\boldsymbol{K}\cdot(\boldsymbol{r}+\boldsymbol{r}^{^{\prime}}) + \phi) c^\dagger_{\boldsymbol{r}} c^\dagger_{\boldsymbol{r}^{^{\prime}}} + \textrm{H.c.}\right] \end{align} \tag{ A71 }$$

as per equation (19).

Kagome lattice

For the kagome lattice, we have

$$\begin{align} \varphi_{\tau \alpha}(\boldsymbol{r}) = \begin{cases} 1 \ \ &\ \boldsymbol{r} \in A \\ -e^{-\frac{2\pi i \alpha}{3}}\ \ &\ \boldsymbol{r} \in B \\ -e^{\frac{2\pi i \alpha}{3}}\ \ &\ \boldsymbol{r}\in C \end{cases} \end{align} \tag{ A72 }$$

so the real space gap function is

$$\begin{align} \Delta(\boldsymbol{r},\boldsymbol{r}^{^{\prime}}) = \begin{cases} \cos(\boldsymbol{K}\cdot(\boldsymbol{r}+\boldsymbol{r}^{^{\prime}}) + \phi ), \ \boldsymbol{r}\in A, \boldsymbol{r}^{^{\prime}}\in B \\ \cos(\boldsymbol{K}\cdot(\boldsymbol{r}+\boldsymbol{r}^{^{\prime}}) +\phi+\pi), \ \boldsymbol{r}\in B, \boldsymbol{r}^{^{\prime}}\in C \\ \cos(\boldsymbol{K}\cdot(\boldsymbol{r}+\boldsymbol{r}^{^{\prime}}) + \phi ), \ \boldsymbol{r}\in C, \boldsymbol{r}^{^{\prime}}\in A \end{cases} \end{align} \tag{ A73 }$$

and hence

$$\begin{align} H &= \sum_{\langle \boldsymbol{r},\boldsymbol{r}^{^{\prime}}\rangle}-t c^\dagger_{\boldsymbol{r}} c_{\boldsymbol{r}^{^{\prime}}} + \Delta^{^{\prime}}_{} \cos(\boldsymbol{K}\cdot(\boldsymbol{r}+\boldsymbol{r}^{^{\prime}}) + \phi) c^\dagger_{\boldsymbol{r}} c^\dagger_{\boldsymbol{r}^{^{\prime}}}\nonumber \\ & \quad + \textrm{h.c.} \end{align} \tag{ A74 }$$

with the appropriate orderings of r and $\boldsymbol{r}^{^{\prime}}$ to take account of the case with a relative π phase—as per equation (18).