---
abstract: 'The rise and fall of the superconducting transition temperature $T_c$ upon tuning carrier density or external parameters, such as pressure or magnetic field, is ubiquitously observed in a wide range of quantum materials. In order to investigate such domes of $T_c$, we go beyond the prototypical Hubbard model, and consider a lattice model of electrons coupled via instantaneous, spatially extended, attractive interactions. Using mean-field theory and functional renormalization group (FRG) methods, we find that for a characteristic interaction range $\ell$, there exists a dome in $T_c$ around $k_F \ell \! \sim \! {\cal O}(1)$. This result can intuitively be understood from the geometric relation between the Fermi surface and the interaction range. We show that our results generally hold for single-band as well as multi-band superconductors, in both two and three dimensions. Our model may be relevant for domes of $T_c$ in dilute weakly coupled superconductors or in engineered cold atom systems.'
author:
- Nazim Boudjada
- Finn Lasse Buessen
- Arun Paramekanti
bibliography:
- 'dome.bib'
title: 'Domes of $T_c$ in superconductors with finite-range attractive interactions'
---

[*Introduction.–*]{} Domes in the superconducting (SC) transition temperature $T_c$, observed across a broad range of quantum materials, typically reflect some form of underlying dynamical competition in the electronic fluid. In heavy fermion compounds [@Stewart_RMP1984; @QSi_Science2010; @WirthSteglich_Nature2016] and iron pnictide materials [@Wen_ARCMP2011; @Si_Nature2016; @Fernandes_2016], for instance, the SC dome emerges around magnetic or nematic quantum critical points (QCPs). In SrTiO$_3$, SC domes may possibly be driven by proximity to a ferroelectric QCP [@GASTIASORO2020168107; @PhysRevMaterials.3.091401; @gastiasoro2020anisotropic; @PhysRevB.100.094504; @PhysRevB.98.024521; @PhysRevB.100.094504; @PhysRevLett.115.247002; @Ahadieaaw0120; @PhysRevB.98.104505; @Neaton2019; @atkinson2017influence]. In the cuprates, even aside from the physics of QCPs, a decrease in the hole concentration can enhance spin-fluctuation mediated pair formation while simultaneously suppressing the superfluid density as a consequence of Mott physics or competing orders. This interplay yields the highest $T_c$ at an optimal doping [@Anderson_RVB2004; @Keimer_Nature2015]. SC near QCPs has also been found in numerical simulations and field theory studies [@Berg_Science2012; @Chowdhury_PRB2015; @Raghu_PRB2015; @Wang_PRL2016; @Torroba_PRB2017; @Berg_ARCMP2019; @Chowdhury2019]. Finally, for ultracold atomic fermions, the highest $T_c$ appears near unitary scattering which marks the BCS-BEC crossover from weak to strong coupling SC [@Randeria_ARCMP2014; @STRINATI_BCSBEC2018].

In this Letter, we discuss a geometric picture of superconducting domes in systems with (non-retarded) finite-range attractive interactions. Our proposal is motivated by the following observation. In a system where electrons attract each other over a fixed characteristic range $\ell$ in real space, the typical momentum transfer in electron-electron scattering processes is $\Delta k \!\sim\! 1/\ell$. Thus, in dilute systems with $k_F \!\ll\! \Delta k$, such interactions can efficiently scatter electrons across any two points on the Fermi surface (FS). In the opposite limit, however, when $k_F \!\gg \! \Delta k$ the aforementioned interactions lead to small-angle scattering, making it more challenging for electrons to explore the full FS. Therefore, the phase space which is accessible in a single electron-scattering event initially increases with the size of the FS, before dropping at high densities when the locality of the interactions in momentum space suppresses global SC. Consequently, a dome-like dependence of $T_c$ on the electron density emerges around some intermediate Fermi momentum $k_F^\star$ which satisfies $k_F^\star \ell \!\sim \! {\cal O}(1)$. The dome thus marks the crossover from predominantly global interactions to local interactions in momentum space. From a real-space perspective, the highest $T_c$ occurs when the interaction range $\ell$ becomes comparable to the interparticle spacing. This geometric picture may be most clearly appreciated in dilute systems, when the FS is far from van Hove singularities. Our work does not address the microscopic origin of such a pairing interaction or the length scale $\ell$, which are important issues in their own right [@PhysRevB.99.094524; @PhysRevB.94.224515; @PhysRevB.98.104505; @Gorkov201604145; @Gorkov], but it is reminiscent of the geometric Mott-Ioffe-Regel criterion which marks the crossover from coherent to incoherent electronic transport [@Hussey_PhilMag2004] without reference to an underlying mechanism.

Expanding on the above physical arguments, we study both single band and multiband models in two and three dimensions (2D and 3D). Our work may be useful as a toy model for systems with critical modes or soft bosons which may induce such long-range attractive interactions, such as for fermions experiencing fluctuating zero-momentum orders. We thus make some qualitative comparisons with results on dilute electron gases in bulk SrTiO$_3$. Models similar in spirit to our study have previously been explored in the context of cuprates [@KunYang_PRB2000], FeSe on SrTiO$_3$ [@Rademaker_2016; @dhlee_ARCMP2018], and ultracold atomic fermions [@Parish_PRB2005]. Previous work has also discussed how density-dependent screening might lead to domes of $T_c$ in SrTiO$_3$ [@PhysRevB.94.224515; @Gorkov201604145; @Gorkov]. However, the universal geometric picture for $T_c$ domes we discuss does not appear to have been highlighted. Our work may also be relevant to ultracold Bose-Fermi mixtures, where $\ell$ could be set by the correlation length associated with the superfluid to Mott insulator transition [@KunYang_PRB2008]. We emphasize, however, that the $T_c$ dome we uncover is not inherently a strong-coupling phenomenon. We therefore employ mean-field theory and FRG methods below to study this problem.

[*Model Hamiltonian. –*]{} We consider a tight-binding Hamiltonian parametrized as $$\mathcal{H}_0=\sum_{{{{\bf{k}}}}\mu\nu} c^\dagger_\mu({{{\bf{k}}}})H^{\mu\nu}_0({{{\bf{k}}}})c_\nu({{{\bf{k}}}}) \,,$$ where $\mu,\nu$ stand for generic orbital and spin indices which give a matrix structure to the Hamiltonian $H_0$. The electrons are assumed to interact via an instantaneous attractive interaction $$\label{eqn:interaction}
\mathcal{H}_{\mathrm{int}}=\frac{1}{2} \int\!\! {{\mathrm{d}}}^d{{{\bf{r}}}}\!\! \int \!\! {{\mathrm{d}}}^d{{{\bf{r}}}}' ~{\cal V}({{{\bf{r}}}}-{{{\bf{r}}}}') \hat{n}({{{\bf{r}}}}) \hat{n}({{{\bf{r}}}}') \,,$$ with $\hat{n}({{{\bf{r}}}}) = \sum_\mu c^\dagger_\mu({{{\bf{r}}}})c_\mu({{{\bf{r}}}})$ being the density operator at position ${{{\bf{r}}}}$, and ${\cal V}({{{{\bf{r}}}}}-{{{{\bf{r}}}}}')  \! < \!0$ being the interaction potential.

Anticipating a singlet superconducting instability, we Fourier transform the interaction to momentum space and focus on the zero center of mass pairing channel, which leads to the effective Hamiltonian $$\label{eqn:MFHamiltonian}
\!\!\!\mathcal{H}^{\rm BCS}_{\mathrm{int}}\!=\! \frac{1}{2N} \!\sum_{{{{\bf{k}}}}{{{\bf{k}}}}'}c^\dagger_\mu({{{\bf{k}}}})c^\dagger_\nu(-{{{\bf{k}}}}) V({{{\bf{k}}}}\!-\!{{{\bf{k}}}}') c_\nu(-{{{\bf{k}}}}')c_\mu({{{\bf{k}}}}') \,,$$

where $N$ is the total number of lattice sites and summation over repeated indices is implied. The interaction $V({{{\bf{k}}}}-{{{\bf{k}}}}')$ is the Fourier transform of ${{\cal V}({{{\bf{r}}}}-{{{\bf{r}}}}')}$. We decouple the interaction, using a Hubbard-Stratonovich transformation, via complex bosonic fields $\Delta_{\mu\nu}({{{\bf{k}}}})$ and integrate out the fermions (see the Supplemental Material (SM) for details). The resulting self-consistent matrix gap equation is given by $$\small
\!\!\!\Delta({{{\bf{k}}}})\!=\! - \frac{1}{N}\! \sum_{{{{\bf{k}}}}'} V({{{\bf{k}}}}\!-\!{{{\bf{k}}}}') U({{{\bf{k}}}}')\frac{\tanh\!\left[\!\frac{E({{{\bf{k}}}}')}{2T}\!\right]}{2 E({{{\bf{k}}}}')} U^\dagger({{{\bf{k}}}}')\Delta({{{\bf{k}}}}') \,,
\label{eq:gapeqn}$$ where $T$ is the temperature, $E({{{\bf{k}}}})$ is a diagonal matrix comprising the square roots of the eigenvalues of ${H_0({{{\bf{k}}}})H^\dagger_0({{{\bf{k}}}})+\Delta({{{\bf{k}}}})\Delta^\dagger({{{\bf{k}}}})}$, and $U({{{\bf{k}}}})$ is the corresponding eigenvector matrix. For a one-band model this expression reduces to the familiar single-gap equation. We assume a Gaussian interaction ${{\cal V}({{{\bf{r}}}}) \!=\! -g_0 {{\mathrm{e}}}^{-|{{{\bf{r}}}}|^2/2\ell^2}}$, so that ${V({{\bf{q}}})\!=\!- g_0(2\pi\ell^2)^{d/2}{{\mathrm{e}}}^{-|{{\bf{q}}}|^2\ell^2/2}}$ [^1] in $d$ spatial dimensions. Here $g_0\!>\! 0$ is the pairing strength, and $\ell$ sets the range of the potential in real space (in units of the lattice constant). For $\ell \!\to\! 0$, the interaction reduces to a Hubbard model, while a large value of $\ell$ favors small momentum scattering. Our results are qualitatively unchanged if we use alternative potentials, such as a Lorentzian or a hard-sphere, with a similar characteristic range $\ell$. To explore the full density dependence and multiband examples, we numerically solve the gap equation for $T_c$ and the momentum dependence of the gap $\Delta_{\mu\nu}({{{\bf{k}}}})$. For a fixed density, we also simultaneously solve for the chemical potential.

[1Band\_2D.pdf]{} (80,33)[(0,1)[9]{}]{} (85,15)[(0,-1)[7]{}]{}

We begin by discussing the geometric origin of domes of $T_c$ for interactions with finite range $\ell$. As illustrated in Fig. \[fig:geometry\](a), for a given momentum point on the 2D FS, the arc length of the FS which lies within the (momentum-space) interaction range $\sim\!\! 1/\ell$ depends on the electron density. Up to $k_F \ell \!=\! 1/2$, the full FS circumference $2\pi k_F$ is accessible in a single electron-scattering event. Beyond this value, the accessible part of the FS shrinks to $2/\ell \! \ll \! 2\pi k_F$ at large $k_F$. The functional form of the geometric arc length $\mathcal{R}(\bar{n})$ versus the electron density $\bar{n}$, displayed in Fig. \[fig:geometry\](b) for a 2D FS with $\ell=1$, exhibits a sharp peak at a density corresponding to $k_F \ell \!=\! 1/2$. Since $\mathcal{R}$ is a geometric measure of the phase space available for Cooper pairs, we expect the peak in ${\mathcal R}$ to be reflected as a dome in $T_c$.

Although $T_c$ is not a purely geometric quantity, depending also on the form of the potential and the energy-dependent density of states, we shall demonstrate that a smoothed version of this geometric maximum generically persists. To illustrate this in a simple one-band example, we compute the mean-field singlet pairing $T_c$ for a 2D square lattice with dispersion $\xi({{{\bf{k}}}}) \!=\!-2t_1 (\cos k_x \!+\! \cos k_y) \!-\! \mu$. Using $\ell\!=\! 5$ and $g_0\!=\! 1$ (in units of $t_1$), we observe a peak in $T_c$ at an electron density $n^\star \!\approx\!0.04$ in the dilute limit as shown in Fig. \[fig:geometry\](c). At higher densities $T_c$ exhibits additional peaks: near half-filling $\bar{n}\!=\!1$, and at $\bar{n}=2-n^\star$. The former peak stems from the combination of an enhanced density of states near the van Hove singularity and an increased geometric overlap with states in the second Brillouin zone (BZ) (see the SM). The latter is a geometric peak arising from small hole pockets near $(\pi,\pi)$, when the hole density becomes dilute.

Our numerical solution of the gap equation shows that domes of $T_c$ also appear in 3D. To investigate the role of the coupling strength $g_0$ in the occurrence of the dome, we compute $k_F^\star \ell$ as a function of $g_0$, where $k_F^\star$ is the angle-averaged Fermi wave vector associated with the density $n^\star$ [^2]. The results for the 2D square lattice as well as for the 3D cubic lattice show that the dome shifts towards smaller densities as $g_0$ is reduced, see Fig. \[fig:geometry\](d). However, we emphasize that the dome persists in the weak coupling limit. Indeed, extrapolating our results to $g_0 \!\to \! 0$, we find a finite value $k_F^\star \ell \! \approx \! 0.2$ in 2D and $k_F^\star \ell \! \approx \! 0.86$ in 3D. At the same time, the ratio of the critical temperature and the Fermi energy $T^\star_c/\epsilon^\star_F$ at the geometric peak remains moderate (with $T^\star_c/\epsilon^\star_F \!<\! 1$ for $g_0=1$, and decreasing for smaller $g_0$), implying that the dome is not a strong-coupling phenomenon.

[*Multiband case. –*]{} Next, we study the generalization of our mean-field results to multiband cases and demonstrate that the geometric interpretation of domes still holds. Furthermore, we shall see that in multiband systems, it is possible to obtain multiple domes of $T_c$ as new FSs appear with increasing density. To this end, we consider the two-orbital model $$\!\!\!\mathcal{H}_0= \sum_{{{{\bf{k}}}}} \! \begin{pmatrix}
X^\dagger_{{{{\bf{k}}}}{\uparrow}} & \! Y^\dagger_{{{{\bf{k}}}}{\uparrow}}
\end{pmatrix}\!\!
\begin{pmatrix}
\xi_{X}({{{\bf{k}}}}) & \! \delta \\
\delta & \! \xi_{Y}({{{\bf{k}}}})
\end{pmatrix}\!\!
\begin{pmatrix}
X_{{{{\bf{k}}}}{\uparrow}} \\
Y_{{{{\bf{k}}}}{\uparrow}}
\end{pmatrix},
\label{eqn:twoorbitals}$$ where ${\xi_X({{{\bf{k}}}}) \!=\! -2t_1\cos k_x\!-\! 2t_2\cos k_y\!-\! \mu}$ and ${\xi_Y({{{\bf{k}}}})\!=\!-2t_2\cos k_x\!-\!2t_1\cos k_y\!-\!\mu}$, with $\delta$ being the momentum-independent interorbital hybridization. As before, we set $t_1\!=\!1$, $\ell\!=\!5$, and $g_0\!=\!1$, and choose, for illustrative purposes, $t_2 \! = \! \delta \!=\! 0.2$; additional examples are discussed in the SM. In the low-density regime ($\bar{n}\ll 0.1$), only one band crosses the Fermi level and the physics is analogous to the single orbital model, i.e. a geometric dome forms at a density corresponding to the optimal value of $k_F\ell$ (orange star in Fig. \[fig:Tc2Bands\]). The parameters are chosen such that the Lifshitz transition, i.e. the appearance of the second band at the Fermi level, occurs near the maximum of the dome. The second band then gives rise to a second geometric peak at a slightly higher density (blue star), yielding an overall double peak structure. We note, however, that for different parameter choices, the Lifshitz transition does not necessarily coincide with the first peak in $T_c$: for example, increasing $\ell$ pushes the geometric peak to lower densities (thus keeping $k_F^\star\ell\sim\mathcal{O}(1)$), but has no impact on the Lifshitz transition point.

![Superconducting transition temperature $T_c$ for a two-orbital model showing a double geometric peak in the low-density regime. Colored stars correspond to the FSs in the insets. Vertical dashed line marks the Lifshitz transition when the second band appears at the Fermi level.[]{data-label="fig:Tc2Bands"}](Tc_2Band_2D.pdf){width="\linewidth"}

[*Functional RG. –*]{} In deriving the gap equation Eqn. , we have explicitly assumed Cooper pair formation in the singlet channel. While yielding a structurally simple, self-consistent mean-field theory, the decoupling comes at the price of being inherently biased to favor the specific type of superconductivity encoded in the ansatz, potentially neglecting any competing phases. When $k_F \ell \!\gg\! 1$, however, different patches on the FS could effectively decouple and many angular momentum pairing channels become quasi-degenerate as seen from the eigenfunctions of the subleading instabilities in the linearized mean field gap equation (see SM). Additionally, attractive interactions could make the systems unstable towards phase separation. Such effects can lead to a breakdown of coherent superconductivity.

![Effective interaction in the patching approximation. Normalized color code shows the value of the flowing interaction vertex $u_T(n_1,n_2,1)$, where $(n_1,n_2)$ enumerate momentum patches around the FS. (a) Density $\bar{n}\!=\!0.17$ at $T\!=\!T_\mathrm{max}$, (b) $\bar{n}\!=\!0.17$ at $T\!=\!T_\mathrm{min}$, (c) $\bar{n}\!=\!0.94$ at $T\!=\!T_\mathrm{min}$.[]{data-label="fig:vertex"}](vertex.pdf){width="\linewidth"}

To investigate this breakdown – or, conversely, justify the mean-field ansatz – we employ the FRG approach which treats all competing interaction channels on equal footing [@Metzner2012; @Gersch2008; @Eberlein2014]. The resulting RG flow equations, which relate the bare interaction as defined in Eqn.  to an effective low-energy theory by continuously tracing its evolution under infinitesimal reductions of the temperature [@Honerkamp2001], naturally have a more complex structure than the self-consistent mean-field equation, and in general cannot be solved exactly. For weak coupling, however, it is sufficient to include only the one-loop contributions to the flow equations for the two-particle interaction, neglecting higher-order processes [@Halboth2000], and to treat the interaction vertex in a momentum patching approximation which resolves its angular dependence around the FS. In this way, a finite set of differential equations is obtained which can be solved numerically to determine the effective interaction vertex $u_T(n_1,n_2,n_3)$, where the $n_i$ enumerate the momentum patches around the FS. The specific choice of momentum patches, as well as the detailed FRG flow equations, are outlined in the SM.

We perform calculations on the single band model at fixed $\ell\!=\!1$ and $g_0\!=\!\frac{3}{2\pi}$, while varying the density $\bar{n}$ to assess the role of competing interaction channels. The RG flow is initialized at an upper temperature $T_\mathrm{max}=4t_1$, which is comparable to the bandwidth, and stopped at a temperature scale $T_\mathrm{min}$ when the maximum component of the vertex exceeds $18t_1$, which is large compared to the bandwidth. The onset of strong interactions at $T_\mathrm{min}$ can then be related to a putative phase transition [^3]. In the dilute limit, $k_F \ell \ll 1$, the bare interaction at $T=T_\mathrm{max}$ has negligible momentum dependence on the FS. The Gaussian profile becomes visible only at slightly larger $\bar{n}$ as seen in Fig. \[fig:vertex\](a). The effective low-temperature vertex, however, for a wide range of $\bar{n}$ is dominated by a distinct attractive interaction between momentum patches which lie on opposite sides of the FS, indicating impending zero-momentum Cooper pair formation, see Fig. \[fig:vertex\](b). However, at large densities $\bar{n}>0.88$, the initial Gaussian profile of the bare interaction sharpens throughout the RG flow as shown in Fig. \[fig:vertex\](c), so that the renormalized $\ell \!\to\! \infty$, and forward scattering gets enhanced. This indicates the breakdown of Cooper pairing and coherent superconductivity.

![Characteristic temperature $T_\mathrm{min}$ for the Gaussian interaction potential with $\ell = 1$, as determined from FRG calculations which include only particle-particle scattering (blue curve) or all interaction channels (orange curve). Curves plotted in opaque colors are computed using $N_p=96$ momentum patches, curves in lighter colors are for $N_p=72$ and $N_p=48$. In regimes I and II, the effective vertex at $T_\mathrm{min}$ captures superconducting pairing, while in regime III the interaction becomes increasingly localized in momentum space leading to a breakdown of SC. Inset shows the geometric dome of $T_c$ at low density.[]{data-label="fig:phasediagramFRG"}](phasediagramFRG.pdf){width="\linewidth"}

To better understand this breakdown, we resolve the role of additional interaction channels by comparing the full FRG calculations with reduced flow equations that only include the particle-particle forward scattering as also captured by the mean-field ansatz. We find that we can divide the FRG phase diagram shown in Fig. \[fig:phasediagramFRG\] into three regimes. In regime I ($\bar{n}<0.58$) the superconducting $T_c$ is suppressed by fluctuations in additional interaction channels. Nevertheless, as shown in the inset to Fig. \[fig:phasediagramFRG\], the full FRG calculation yields a dome of $T_c$, in qualitative agreement with the simplified mean-field approach. In regime II ($0.58<\bar{n}<0.88$) on the other hand, unlike what is seen for the attractive Hubbard model, the finite-range character of the interactions leads to an enhancement of $T_\mathrm{min}$ by the additional interaction channels. Finally, in regime III ($\bar{n}>0.88$), mean-field theory formally yields a finite $T_c$, while the full FRG approach reveals the breakdown of superconductivity. We tentatively identify this regime, where the renormalized $\ell \! \to \! \infty$, with phase separation induced by extended attractive interactions.

[*Conclusion. –*]{} In this Letter, we have provided a geometric phase space argument for the formation of $T_c$ domes in systems with spatially extended interactions. We have shown that for multiband systems a scenario with two or more domes can arise naturally. In order to apply this picture to 3D bulk SrTiO$_3$, we note that the first dome with maximum transition temperature $T_c \!\approx\!0.2$K is centered at a density $\bar{n} \!\approx\! 1.2 \times 10^{18}\;{\rm cm}^{-3}$ with a Fermi energy $\epsilon_F\!\approx\! 2$meV. Demanding $k_F \ell \!\sim\! 1$ at the center of the dome yields $\ell \! \sim\! 30$Å, while requiring $T_c/\epsilon_F \! \sim \! 10^{-2}$ at this point fixes $g_0 \!\sim\! 4.5$meV. The inferred length scale $\ell$ may reflect the range of attractive interactions between polaron quasiparticles which have been reported in bulk SrTiO$_3$ [@Swartz_PNAS2018] and its interfaces [@Baumberger_NMat2016; @Strocov_NComms2016]. The microscopic theory of SC of such dilute polarons remains an open issue. It would be interesting to explore such $T_c$ domes in a wider range of experimental systems including atomic Bose-Fermi mixtures, and to extend the FRG results by incorporating the frequency dependence of the interaction vertex. Such studies may shed light on the interplay of spatially extended interactions with retardation effects in driving SC near QCPs.

We thank M. M. Scherer for discussions. This work was funded by NSERC of Canada and FRQNT of Quebec. The numerical simulations were performed on the JURECA cluster at the Forschungszentrum Juelich, and on the Cedar and Niagara clusters enabled by support provided by Compute Ontario, SciNet, Westgrid, and Compute Canada. SciNet is funded by: the Canada Foundation for Innovation; the Government of Ontario; Ontario Research Fund - Research Excellence; and the University of Toronto.

Supplemental Material {#supplemental-material .unnumbered}
=====================

Derivation of multiband gap equation {#derivationGE}
------------------------------------

We start by writing the imaginary-time action for free fermions: $$S_0=\frac{1}{V^2}\int{{\mathrm{d}}}^d{{{\bf{r}}}}\;\int{{\mathrm{d}}}^d{{{\bf{r}}}}'\int_0^\beta{{\mathrm{d}}}\tau\;\bar{\psi}_\mu({{{\bf{r}}}},\tau)[\partial_\tau\delta_{\mu\nu}-H_0^{\mu\nu}({{{\bf{r}}}},{{{\bf{r}}}}')]\psi_\nu({{{\bf{r}}}}',\tau),$$ where $V=Na^d$ is the volume of the $d$-dimensional cubic system with $N$ sites of lattice constant $a$. Working in units where $a=1$, we can Fourier transform $S_0$ to momentum and Matsubara frequency space for a translationally invariant system: $$\begin{aligned}
S_0&=&\frac{1}{N}\sum_{{{{\bf{k}}}}\omega_n}\bar{\psi}_\mu({{{\bf{k}}}},{{\mathrm{i}}}\omega_n)[{{\mathrm{i}}}\omega_n\delta_{\mu\nu}-H_0^{\mu\nu}({{{\bf{k}}}})]\psi_\nu({{{\bf{k}}}},{{\mathrm{i}}}\omega_n)\nonumber\\
&=&\frac{1}{2N}\sum_{{{{\bf{k}}}}\omega_n}\left(\bar{\psi}_\mu({{{\bf{k}}}},{{\mathrm{i}}}\omega_n)\underbrace{[{{\mathrm{i}}}\omega_n\delta_{\mu\nu}-H_0^{\mu\nu}({{{\bf{k}}}})]}_{G_{0p}^{-1}({{{\bf{k}}}},{{\mathrm{i}}}\omega_n)}\psi_\nu({{{\bf{k}}}},{{\mathrm{i}}}\omega_n)+\psi_\mu(-{{{\bf{k}}}},-{{\mathrm{i}}}\omega_n)\underbrace{[{{\mathrm{i}}}\omega_n\delta_{\mu\nu}+H_0^{\nu\mu}(-{{{\bf{k}}}})]}_{G_{0h}^{-1}(-{{{\bf{k}}}},-{{\mathrm{i}}}\omega_n)}\bar{\psi}_\nu(-{{{\bf{k}}}},-{{\mathrm{i}}}\omega_n)\right)\nonumber\\
&=&\frac{1}{2N}\sum_{{{{\bf{k}}}}\omega_n}
\begin{pmatrix}\bar{\psi}({{{\bf{k}}}},{{\mathrm{i}}}\omega_n) & \psi(-{{{\bf{k}}}},-{{\mathrm{i}}}\omega_n)\end{pmatrix}
\begin{pmatrix}G_{0p}^{-1}({{{\bf{k}}}},{{\mathrm{i}}}\omega_n) & 0 \\ 
0 & G_{0h}^{-1}(-{{{\bf{k}}}},-{{\mathrm{i}}}\omega_n)\end{pmatrix}
\begin{pmatrix}\psi({{{\bf{k}}}},{{\mathrm{i}}}\omega_n) \\ \bar{\psi}(-{{{\bf{k}}}},-{{\mathrm{i}}}\omega_n)\end{pmatrix} \,,\end{aligned}$$ where $G^{-1}_{0p}$ and $G^{-1}_{0h}$ are the matrix non-interacting Green’s functions, neglecting self-energy corrections. The instantaneous interaction is, in real-space: $$S_{\mathrm{int}}=\frac{1}{2V^2}\int{{\mathrm{d}}}^d{{{\bf{r}}}}\;\int{{\mathrm{d}}}^d{{{\bf{r}}}}'\int_0^\beta{{\mathrm{d}}}\tau\;\bar{\psi}_\mu({{{\bf{r}}}},\tau)\bar{\psi}_\nu({{{\bf{r}}}}',\tau)\mathcal{V}({{{\bf{r}}}}-{{{\bf{r}}}}')\psi_\nu({{{\bf{r}}}}',\tau)\psi_\mu({{{\bf{r}}}},\tau) \,,$$ and in ${{{\bf{k}}}}$-space, if we only keep zero centre of mass momentum terms: $$S_{\mathrm{int}}=\frac{1}{2\beta N^2}\sum_{{{{\bf{k}}}}{{{\bf{k}}}}'}\sum_{\omega_n\omega_m}\bar{\psi}_\mu({{{\bf{k}}}},{{\mathrm{i}}}\omega_n)\bar{\psi}_\nu(-{{{\bf{k}}}},-{{\mathrm{i}}}\omega_n)V({{{\bf{k}}}}-{{{\bf{k}}}}')\psi_\nu(-{{{\bf{k}}}}',-{{\mathrm{i}}}\omega_m)\psi_\mu({{{\bf{k}}}}',{{\mathrm{i}}}\omega_m) \,,$$ with $V({{{\bf{k}}}}-{{{\bf{k}}}}')<0$ the Fourier transform of $\mathcal{V}({{{\bf{r}}}})$. We now introduce the complex fields $\Delta_{\mu\nu}({{{\bf{k}}}})$: $$\begin{aligned}
{{\mathrm{e}}}^{-S_{\mathrm{int}}}\propto\int\mathcal{D}[\bar{\Delta},\Delta]\;\exp\bigg(&-\frac{1}{2N}\sum_{{{{\bf{k}}}}\omega_n}\Big(\frac{\beta}{N}\sum_{{{{\bf{k}}}}'}\Delta_{\mu\nu}({{{\bf{k}}}}')\Delta_{\nu\mu}^*({{{\bf{k}}}})F({{{\bf{k}}}}-{{{\bf{k}}}}') \nonumber\\
&+\psi_\nu(-{{{\bf{k}}}},-{{\mathrm{i}}}\omega_n)\Delta_{\nu\mu}^*({{{\bf{k}}}})\psi_\mu({{{\bf{k}}}},{{\mathrm{i}}}\omega_n)+\bar{\psi}_\mu({{{\bf{k}}}},{{\mathrm{i}}}\omega_n)\Delta_{\mu\nu}({{{\bf{k}}}})\bar{\psi}_\nu(-{{{\bf{k}}}},-{{\mathrm{i}}}\omega_n)\Big)\bigg)\end{aligned}$$ with $F({{{\bf{k}}}}-{{{\bf{k}}}}')=\frac{1}{V}\int{{\mathrm{d}}}^d{{{\bf{r}}}}\;\frac{{{\mathrm{e}}}^{{{\mathrm{i}}}({{{\bf{k}}}}-{{{\bf{k}}}}')\cdot {{{\bf{r}}}}}}{\mathcal{V}({{{\bf{r}}}})}$. The total partition function, $\mathcal{Z}=\int\mathcal{D}[\bar{\psi},\psi]{{\mathrm{e}}}^{-(S_0+S_{\mathrm{int}})}$ up to normalization constants, becomes quadratic in the fermion fields: $$\begin{aligned}
\hspace{-5mm}
\mathcal{Z}&=&\int\mathcal{D}[\bar{\Delta},\Delta]\exp\left(-\frac{\beta}{2N^2}\sum_{{{{\bf{k}}}}{{{\bf{k}}}}'}\Delta_{\mu\nu}({{{\bf{k}}}}')\Delta_{\nu\mu}^*({{{\bf{k}}}})F({{{\bf{k}}}}-{{{\bf{k}}}}')\right)\nonumber\\
&\times&\int\mathcal{D}[\bar{\psi},\psi]\exp\left(-\frac{1}{2N}\sum_{{{{\bf{k}}}}\omega_n}
\underbrace{\begin{pmatrix}\bar{\psi}({{{\bf{k}}}},{{\mathrm{i}}}\omega_n) & \psi(-{{{\bf{k}}}},-{{\mathrm{i}}}\omega_n)\end{pmatrix}}_{\bar{\Psi}({{{\bf{k}}}},{{\mathrm{i}}}\omega_n)}
\underbrace{\begin{pmatrix}G_{0p}^{-1}({{{\bf{k}}}},{{\mathrm{i}}}\omega_n) & \Delta({{{\bf{k}}}}) \\ 
    \Delta^\dagger({{{\bf{k}}}}) & G_{0h}^{-1}(-{{{\bf{k}}}},-{{\mathrm{i}}}\omega_n)\end{pmatrix}}_{\mathcal{G}^{-1}({{{\bf{k}}}},{{\mathrm{i}}}\omega_n)}
\underbrace{\begin{pmatrix}\psi({{{\bf{k}}}},{{\mathrm{i}}}\omega_n) \\ \bar{\psi}(-{{{\bf{k}}}},-{{\mathrm{i}}}\omega_n)\end{pmatrix}}_{\Psi({{{\bf{k}}}},{{\mathrm{i}}}\omega_n)}\right) \,.\end{aligned}$$ And we can proceed by integrating out the fermions and obtain the effective action $\mathcal{Z}=\int\mathcal{D}[\bar{\Delta},\Delta]{{\mathrm{e}}}^{-S_{\mathrm{eff}}}$: $$S_{\mathrm{eff}}=\underbrace{\frac{\beta}{2N^2}\sum_{{{{\bf{k}}}}{{{\bf{k}}}}'}\Delta_{\mu\nu}({{{\bf{k}}}}')\Delta_{\nu\mu}^*({{{\bf{k}}}})F({{{\bf{k}}}}-{{{\bf{k}}}}')}_{S_1}+\underbrace{\frac{1}{2N}\sum_{{{{\bf{k}}}}\omega_n}\mathrm{tr}\log\begin{pmatrix}{{\mathrm{i}}}\omega_n\mathbb{I}-H_0({{{\bf{k}}}}) & \Delta({{{\bf{k}}}}) \\ 
    \Delta^\dagger({{{\bf{k}}}}) & {{\mathrm{i}}}\omega_n\mathbb{I}+H_0^\mathrm{T}(-{{{\bf{k}}}})\end{pmatrix}}_{S_2},$$ To obtain the equation of motion, we need to set $\frac{\delta S_{\mathrm{eff}}}{\delta \Delta^*_{\sigma\lambda}({{{\bf{p}}}})}=0$. Varying $S_1$ is straighforward: $$\frac{1}{\beta}\frac{\delta S_1}{\delta \Delta^*_{\sigma\lambda}({{{\bf{p}}}})}=\frac{1}{2N^2}\sum_{{{{\bf{k}}}}{{{\bf{k}}}}'}\delta_{\nu\sigma}\delta_{\mu\lambda}\delta^{(d)}({{{\bf{p}}}}-{{{\bf{k}}}})\Delta_{\mu\nu}({{{\bf{k}}}}')F({{{\bf{k}}}}-{{{\bf{k}}}}')=\frac{1}{2N}\sum_{{{{\bf{k}}}}'}\Delta_{\lambda\sigma}({{{\bf{k}}}}')F({{{\bf{p}}}}-{{{\bf{k}}}}') \,.$$ Varying $S_2$ leads to: $$\begin{aligned}
\frac{1}{\beta}\frac{\delta S_2}{\delta \Delta^*_{\sigma\lambda}({{{\bf{p}}}})}&=\frac{1}{2\beta N}\sum_{{{{\bf{k}}}}\omega_n}\mathrm{tr}\left(\frac{\delta\log(\mathcal{G}^{-1}({{{\bf{k}}}},{{\mathrm{i}}}\omega_n))}{\delta\mathcal{G}^{-1}({{{\bf{k}}}},{{\mathrm{i}}}\omega_n)}\frac{\delta\mathcal{G}^{-1}({{{\bf{k}}}},{{\mathrm{i}}}\omega_n)}{\delta \Delta^*_{\sigma\lambda}({{{\bf{p}}}})}\right) \nonumber\\
&=\frac{1}{2\beta N}\sum_{{{{\bf{k}}}}\omega_n}\left[\mathcal{G}({{{\bf{k}}}},{{\mathrm{i}}}\omega_n)\delta^{(d)}({{{\bf{p}}}}-{{{\bf{k}}}})\begin{pmatrix} 0 & 0 \\ \delta_{\sigma\lambda} & 0 \end{pmatrix}\right] \nonumber\\
&=\frac{1}{2\beta}\sum_{\omega_n}\left[\mathcal{G}({{{\bf{p}}}},{{\mathrm{i}}}\omega_n)\right]_{\lambda\sigma}^{12} \,.
\label{eqn:dS2}\end{aligned}$$ $\left[\mathcal{G}({{{\bf{p}}}},{{\mathrm{i}}}\omega_n)\right]_{\lambda\sigma}^{12}$ refers to the $(\lambda,\sigma)$ matrix element of the $(1,2)$ block (i.e. top right block) of the $\mathcal{G}$ matrix. In order to invert a matrix containing square block matrices, we make use of the following identity:

$$\begin{pmatrix} A & B \\ 
C & D\end{pmatrix}^{-1}=\begin{pmatrix} A^{-1}+A^{-1}B(D-CA^{-1}B)^{-1}CA^{-1} & -A^{-1}B(D-CA^{-1}B)^{-1} \\ 
-(D-CA^{-1}B)^{-1}CA^{-1} & (D-CA^{-1}B)^{-1}\end{pmatrix} \,.$$

Focusing on the top right corner, we obtain after some algebra and using the particle-hole symmetry of the Hamiltonian: $$\begin{aligned}
\left[\mathcal{G}({{{\bf{p}}}},{{\mathrm{i}}}\omega_n)\right]^{12}&=&-[\omega_n^2+{{\mathrm{i}}}\omega_n(H_0({{{\bf{p}}}})-\Delta({{{\bf{p}}}}) H_0^\mathrm{T}(-{{{\bf{p}}}})\Delta^{-1}({{{\bf{p}}}}))+\Delta({{{\bf{p}}}})H_0^\mathrm{T}(-{{{\bf{p}}}})\Delta^{-1}({{{\bf{p}}}})H_0({{{\bf{p}}}})+\Delta({{{\bf{p}}}})\Delta^\dagger({{{\bf{p}}}})]^{-1}\Delta({{{\bf{p}}}})\nonumber\\
&=&-({{\mathrm{i}}}\omega_n\mathbb{I}-\mathcal{E}({{{\bf{p}}}}))^{-1}({{\mathrm{i}}}\omega_n\mathbb{I}+\mathcal{E}({{{\bf{p}}}}))^{-1}\Delta({{{\bf{p}}}})\end{aligned}$$ and $\mathcal{E}^2({{{\bf{p}}}})\equiv H_0({{{\bf{p}}}})H^\dagger_0({{{\bf{p}}}})+\Delta({{{\bf{p}}}})\Delta^\dagger({{{\bf{p}}}})$. This matrix can be diagonalized via $\mathcal{E}^2({{{\bf{p}}}})=U({{{\bf{p}}}})E^2({{{\bf{p}}}})U^\dagger({{{\bf{p}}}})$ and the Matsubara frequency summation performed: $$\begin{aligned}
\frac{1}{\beta}\frac{\delta S_2}{\delta \Delta^*_{\sigma\lambda}({{{\bf{p}}}})}&=\frac{1}{2\beta}\sum_{{{\mathrm{i}}}\omega_n}-({{\mathrm{i}}}\omega_n\mathbb{I}-\mathcal{E}({{{\bf{p}}}}))_{\sigma m}^{-1}({{\mathrm{i}}}\omega_n\mathbb{I}+\mathcal{E}({{{\bf{p}}}}))_{m\alpha}^{-1}\Delta_{\alpha\lambda}({{{\bf{p}}}}) \nonumber\\
&=-\frac{1}{2\beta}\sum_{{{\mathrm{i}}}\omega_n}U_{\sigma m}({{{\bf{p}}}})({{\mathrm{i}}}\omega_n-E({{{\bf{p}}}}))_{m}^{-1}({{\mathrm{i}}}\omega_n+E({{{\bf{p}}}}))_m^{-1}U_{m\alpha}^*({{{\bf{p}}}})\Delta_{\alpha\lambda}({{{\bf{p}}}}) \nonumber\\
&=\frac{1}{2}U_{\sigma m}({{{\bf{p}}}})\frac{1}{2E_m({{{\bf{p}}}})}\tanh\left(\frac{\beta E_m({{{\bf{p}}}})}{2}\right)U_{m\alpha}^*({{{\bf{p}}}})\Delta_{\alpha\lambda}({{{\bf{p}}}})\end{aligned}$$ with implied sums over repeated indices. In matrix notation: $$\frac{1}{\beta}\frac{\delta S_\mathrm{eff}}{\delta \Delta^\dagger({{{\bf{p}}}})}=\frac{1}{2N}\sum_{{{{\bf{k}}}}'} F({{{\bf{p}}}}-{{{\bf{k}}}}')\Delta({{{\bf{k}}}}')+\frac{1}{2}U({{{\bf{p}}}})\frac{1}{2E({{{\bf{p}}}})}\tanh\left(\frac{\beta E({{{\bf{p}}}})}{2}\right)U^\dagger({{{\bf{p}}}})\Delta({{{\bf{p}}}})=0 \,.$$ For an inversion symmetric scattering potential $\mathcal{V}({{{\bf{r}}}})=\mathcal{V}(-{{{\bf{r}}}})$, this can be rewritten in terms of the Fourier transform $V({{{\bf{k}}}})$: $$\begin{aligned}
\Delta_{\mu\nu}({{{\bf{k}}}})&=-\frac{1}{N}\sum_{{{{\bf{k}}}}'}\sum_{m\lambda} V({{{\bf{k}}}}-{{{\bf{k}}}}') U_{\mu m}({{{\bf{k}}}}')\frac{1}{2E_m({{{\bf{k}}}}')}\tanh\left(\frac{E_m({{{\bf{k}}}}')}{2T}\right)U^*_{m\lambda}({{{\bf{k}}}}')\Delta_{\lambda\nu}({{{\bf{k}}}}') \nonumber\\
&=\sum_{{{{\bf{k}}}}'}V({{{\bf{k}}}}-{{{\bf{k}}}}')M_{\mu\lambda}({{{\bf{k}}}}')\Delta_{\lambda\nu}({{{\bf{k}}}}') \,,\end{aligned}$$ with $M_{\mu\lambda}({{{\bf{k}}}}')\equiv-\frac{1}{N}\sum_m U_{\mu m}({{{\bf{k}}}}')\frac{1}{2E_m({{{\bf{k}}}}')}\tanh\left(\frac{E_m({{{\bf{k}}}}')}{2T}\right)U^*_{m\lambda}({{{\bf{k}}}}')$. Numerically, it is useful to write this as a matrix equation at temperature $T$: $$\vec{\Delta}(T)=\mathcal{M}(T)\;\vec{\Delta}(T)\iff\Delta_i(T)=\mathcal{M}_{ij}(T)\Delta_j(T) \,.$$ For a Hamiltonian comprising $b$ bands and discretized on a momentum mesh with $N$ points, the indices ${i,j\in[1,...,b^2N]}$ which makes the matrix $\mathcal{M}$ of dimensions $b^2N\times b^2N$. More explicitly: $$\begin{pmatrix}
\Delta_{11}({{{\bf{k}}}}_1) \\ \Delta_{12}({{{\bf{k}}}}_1) \\ \vdots \\ \Delta_{1b}({{{\bf{k}}}}_1) \\ \vdots \\ \Delta_{bb}({{{\bf{k}}}}_1) \\ \vdots \\ \Delta_{bb}({{{\bf{k}}}}_N)
\end{pmatrix}
\!\!=\!\!
\begin{pmatrix}
V({{{\bf{k}}}}_1-{{{\bf{k}}}}_1)[M({{{\bf{k}}}}_1)]\otimes\mathbb{I}_b & V({{{\bf{k}}}}_1-{{{\bf{k}}}}_2)[M({{{\bf{k}}}}_2)]\otimes\mathbb{I}_b & \hdots & V({{{\bf{k}}}}_1-{{{\bf{k}}}}_N)[M({{{\bf{k}}}}_N)]\otimes\mathbb{I}_b\\
V({{{\bf{k}}}}_2-{{{\bf{k}}}}_1)[M({{{\bf{k}}}}_1)]\otimes\mathbb{I}_b & V({{{\bf{k}}}}_2-{{{\bf{k}}}}_2)[M({{{\bf{k}}}}_2)]\otimes\mathbb{I}_b & \hdots & V({{{\bf{k}}}}_2-{{{\bf{k}}}}_N)[M({{{\bf{k}}}}_N)]\otimes\mathbb{I}_b\\
\vdots & \vdots & \ddots & \vdots \\
V({{{\bf{k}}}}_N-{{{\bf{k}}}}_1)[M({{{\bf{k}}}}_1)]\otimes\mathbb{I}_b & V({{{\bf{k}}}}_N-{{{\bf{k}}}}_2)[M({{{\bf{k}}}}_2)]\otimes\mathbb{I}_b & \hdots & V({{{\bf{k}}}}_N-{{{\bf{k}}}}_N)[M({{{\bf{k}}}}_N)]\otimes\mathbb{I}_b
\end{pmatrix}
\!\!
\begin{pmatrix}
\Delta_{11}({{{\bf{k}}}}_1) \\ \Delta_{12}({{{\bf{k}}}}_1) \\ \vdots \\ \Delta_{1b}({{{\bf{k}}}}_1) \\ \vdots \\ \Delta_{bb}({{{\bf{k}}}}_1) \\ \vdots \\ \Delta_{bb}({{{\bf{k}}}}_N)
\end{pmatrix} \,,
\label{eqn:matrixDelta}$$ where $\mathbb{I}_b$ is the $b\times b$ identity matrix.

[GapTc1B.pdf]{}

At $T=T_c$, $E_m({{{\bf{k}}}})\rightarrow\xi_m({{{\bf{k}}}})$ and $U_{\mu m}({{{\bf{k}}}})\rightarrow W_{\mu m}({{{\bf{k}}}})$ where $\xi_m({{{\bf{k}}}})=W^*_{m\lambda}({{{\bf{k}}}})H_0^{\lambda\sigma}({{{\bf{k}}}})W_{\sigma m}({{{\bf{k}}}})$ such that $M$ no longer depends on $\Delta$. Eqn. reduces to an eigenvalue equation and $T_c$ is obtained when the largest eigenvalue of $\mathcal{M}$ reaches 1 (for $T\ll T_c$ all the eigenvalues are larger than 1 while for $T\gg T_c$ all the eigenvalues vanish.)

At $T=0$, the equation is non-linear as the matrix $M$ depends on $\Delta$ and we must solve by (i) guessing an initial $\Delta^{(0)}({{{\bf{k}}}})$, (ii) diagonalizing $H_0({{{\bf{k}}}})H_0^\dagger({{{\bf{k}}}})+\Delta^{(0)}({{{\bf{k}}}})[\Delta^{(0)}({{{\bf{k}}}})]^\dagger=U({{{\bf{k}}}})[E^{(0)}({{{\bf{k}}}})]^2U^{-1}({{{\bf{k}}}}')$, (iii) constructing the $\mathcal{M}$ matrix and (iv) multiplying by the ‘vectorized’ $\vec{\Delta}^{(0)}$ to obtain a new vector $\vec{\Delta}^{(1)}$ which can then be used to repeat the procedure. The gap function $\Delta(T=0)$ is the ‘fixed point’ of this equation.

Finally, since the instability is expected near the Fermi momenta (i.e. where the denominator $E_m({{{\bf{k}}}})\rightarrow 0$) it is useful to only store momenta within a given range of ${{{\bf{k}}}}_F$. The lengthscale $\ell^{-1}$ provides a natural cutoff and we found that only keeping momenta within $\pm 3\ell^{-1}$ of ${{{\bf{k}}}}_F$ is sufficient to reach convergent results.

Momentum structure of the gap functions
---------------------------------------

At $T=T_c$, the spectrum of $\mathcal{M}$ tells us about the modes of the gap. The largest eigenvalue is the lowest energy state into which the system will condense first and the corresponding eigenvector shows the momentum dependence of the gap. When $\ell\rightarrow\infty$, all eigenvalues converge to 1 and the corresponding eigenvectors become localized to single momentum points on the Fermi surface. In that case, any linear combination of the eigenvectors would be a solution to the gap equation and all momenta condense simultaneously. However, for finite values of $\ell$, this doesn’t happen and we generally have one eigenvalue reaching unity before the others. The closest eigenvalues correspond to eigenvectors with higher energy and in general get closer to each other for high densities (as this is equivalent to increasing $\ell$). For a generic density and pairing lengthscale shown in Fig. \[fig:gapTc1B\], we show that the largest eigenvalue is nondegenerate and has $s$-wave symmetry, while the next eigenvalues are doubly degenerate with $p_x$ and $p_y$ symmetry and the next two have $d$-wave symmetry. In Fig. \[fig:gaps\], we show the gaps at various densities at $T=0$, which are obtained by solving the non-linear gap equation.

[Gaps.pdf]{}(2,94)[(a)]{}(53,94)[(b)]{}(2,45)[(c)]{}(53,45)[(d)]{}

The van Hove peak
-----------------

[MFGeometric.pdf]{}(16,95)[(a)]{} (62,95)[(b)]{}(16,44)[(c)]{}

The critical temperature in 2D over all densities is shown in Fig. \[fig:Tc2D\] (a) for the same choice of parameters as Fig. \[fig:geometry\] of the main text. For a one band model with isotropic dispersion there are three peaks. The first peak (dashed line) is the electron geometric peak and was discussed in the main text. We also have the hole version of this geometric peak marked with a dot-dashed line. The corresponding gaps at $T=0$ are shown in Figs. \[fig:gaps\](a) and (d), respectively. Since $1/\ell$ is much larger than the size of the hole pocket, the gap doesn’t peak on the Fermi surface but at the center of the pockets, i.e. the M points. This is similar to the electron pockets where the gap peaks at the $\Gamma$ point in the radial direction. The middle peak (dotted line) is exactly at half-filling when the van Hove singularity occurs in two-dimensions: the density of states (DOS) is largest and we expect an enhancement of $T_c$ from a BCS-like picture. The amplitude of the gap $|\Delta({{{\bf{k}}}},T=0)|$ is shown in Fig. \[fig:gaps\] (c) with the FS overlaid on top. Tangentially to the Fermi surface, it peaks at the van Hove points $(\pi,0)$ and $(0,\pi)$ while in the perpendicular direction it peaks at ${{{\bf{k}}}}_F$ and decays a distance $\sim 1/\ell$ away from the FS. Although the existence of the van Hove peak is expected independently of the lengthscale $\ell$, its appearance overlaps with a different, $\ell$-dependent effect: As the van Hove point is approached, Umklapp processes become allowed and the scattering phase space is enhanced for $|{{{\bf{k}}}}_F-({{{\bf{k}}}}'_F\pm{{{\bf{G}}}})|\ell\sim 1$ where ${{{\bf{G}}}}$ is a reciprocal lattice vector. This is also seen in the geometric weight calculation where $\mathcal{R}$ goes up again as parts of the FS from the neighboring BZ become accessible (Fig. \[fig:Tc2D\] (b)). This situation is shown pictorially in Fig. \[fig:Tc2D\] (c) where a point on the FS of the left Brillouin zone starts to ‘sense’ its neighbor on the other side of the BZ because of the finite potential range $\sim1/\ell$.

Anisotropic dispersions and $\ell$ dependence
---------------------------------------------

[Tc\_1Band\_2D\_L\_Anisotropic.pdf]{}(16,47)[(a)]{} (62,47)[(b)]{}

As we tune the value of $\ell$ we observe that as $\ell$ is increased the location of the geometric peak moves to lower densities as expected (thus keeping $k_F^\star\ell\sim1$). On the other hand, although the van Hove peak remains at half filling (Fig. \[fig:Tc1B2DLAni\] (a)), its tail gets sharper for small $\ell$ because Umklapp scattering only kicks in at increasingly large fillings. For an anisotropic FS with dispersion ${\xi({{{\bf{k}}}})=-2t_1\cos(kx)-2t_2\cos(k_y)-\mu}$ and $t_1\neq t_2$, the location of the van Hove peak shifts to lower densities since the elongated part of the elliptical FS hits the BZ boundary at densities $\bar{n}<1$. For a very elliptical FS ($t_2/t_1\ll 1$ or $t_2/t_1\gg1$), the van Hove point is at smaller densities and $k_F$ acquires a strong angular modulation ${k_F^{\mathrm{long}}\gg k_F^{\mathrm{short}}}$. In turn, this introduces a new condition that as soon as $k_F^{\mathrm{long}}\ell\gg1$, the available phase space for scattering starts to decrease even if ${k_F^{\mathrm{long}}\ell\lesssim1}$. In this scenario, the van Hove peak shifts to lower densities while the geometric peak is pushed to higher densities and the two eventually merge into a single peak (Fig. \[fig:Tc1B2DLAni\] (b)). The peak value of $T_c$ is enhanced when $\ell$ is large and it can be shown analytically that the mean-field $T_c$ reaches an upper limit of $g_0/4$ when $\ell\rightarrow\infty$ ($V({{{\bf{k}}}})\rightarrow-g_0\delta^{(d)}({{{\bf{k}}}})$); this value is however suppressed when vertex corrections are included as discussed in the main text. In all cases, $T_c$ hardly changes by a factor of 2-3 over many orders of magnitude of $\bar{n}$.

Finite potential difference
---------------------------

![Critical temperature for a two-orbital model (black curve) with a finite potential difference $2V_0=2$ and a single orbital model (red curve) showing a perfect agreement at low densities when only one band is populated. The insets show the Fermi surfaces at the correspondingly colored stars. []{data-label="fig:V02D"}](Tc_2Band_2D_V0.pdf){width="\linewidth"}

In this section, we consider a modified version of the Hamiltonian presented in Eqn. of the main text $$\!\!\!\mathcal{H}_0=\sum_{{{{\bf{k}}}}} \! \begin{pmatrix}
X^\dagger_{{{{\bf{k}}}}{\uparrow}} & \! Y^\dagger_{{{{\bf{k}}}}{\uparrow}}
\end{pmatrix}\!\!
\begin{pmatrix}
\xi_{X}({{{\bf{k}}}})+V_0 & \! 0 \\
0 & \! \xi_{Y}({{{\bf{k}}}})-V_0
\end{pmatrix}\!\!
\begin{pmatrix}
X_{{{{\bf{k}}}}{\uparrow}} \\
Y_{{{{\bf{k}}}}{\uparrow}}
\end{pmatrix},$$

and we study the impact of a finite potential difference between the two orbitals, keeping $t_1=t_2=1$. This corresponds to two independent bands separated in energy by $2V_0$. The resulting $T_c$ is illustrated in Fig. \[fig:V02D\] with the Fermi surfaces at the three peak densities in electron-doped regime shown in the insets and the dashed red curve showing the one-orbital model with the same parameters.

In the low-density regime, the two curves agree exactly, showing that only the lower band dictates $T_c$: a first geometric peak (marked by a green star) is reached at $\bar{n}=n^\star\approx0.046$ which corresponds to $k_F^\star\ell\approx2.85$. Near the Lifshitz transition (vertical dashed curve), the one-orbital and two-orbital models start to diverge since the higher band crosses the Fermi level and starts to contribute to $T_c$. A second geometric peak (orange star) is reached precisely when the new band’s Fermi wavevector is such that $k_F^\star\ell\approx2.85$. At higher densities, the lower band reaches the van Hove point and we see the corresponding peak (blue star). As for the one-orbital case, particle-hole symmetry dictates $T_c(\bar{n})=T_c(4-\bar{n})$ and an exact copy of the three peaks is obtained in the hole-doped regime $\bar{n}>2$.

FRG flow equations {#sec:appendix:frg}
------------------

![Partitioning of the Brillouin zone into a set of patches $\mathcal{P}_1,\dots,\mathcal{P}_N$, illustrated for $N=24$. All momentum points within a single patch $\mathcal{P}_m$ are projected onto the representative momentum ${{{\bf{p}}}}_m$ which lies on the Fermi surface. []{data-label="fig:appendix:patching"}](patchingScheme.pdf){width="0.9\linewidth"}

![Diagrammatic representation of contributions to the flow of the effective interaction. Singling out particle-particle scattering (first term) is equivalent to a Bethe-Salpeter like resummation and generates results on the level of self-consistent mean-field theory. The inclusion of additional direct particle-hole (terms two to four) and crossed particle-particle (last term) scattering allows to model the interplay of competing interaction channels.[]{data-label="fig:appendix:FRGchannels"}](flowequation.pdf){width="\linewidth"}

The general form of the effective two-particle interaction vertex is given by $$\label{eq:appendix:frg:vertex}
V(K_1',K_2';K_1,K_2) \,,$$ where the parameters $K_n$ are composite indices denoting tuples $({{{\bf{k}}}}_n,\omega_n,\alpha_n)$ of momentum, Matsubara frequency, and spin, respectively. In this most general form, the fermionic interaction vertex must be antisymmetric under the pairwise exchange of its arguments. For the SU(2)-symmetric model at hand, however, it is more convenient to constrain the effective interaction to a form which is inherently encoded to be SU(2) symmetric; To this end, we parameterize the effective interaction by two terms – spin-conserving and spin-exchange terms – which span a full basis for SU(2) invariant interactions: $$\begin{aligned}
V(K_1',K_2';K_1,K_2)&= U(k_1',k_2';k_1,k_2)\delta_{\alpha_1' \alpha_1}\delta_{\alpha_2'\alpha_2} \nonumber\\
&- U(k_1',k_2';k_2,k_1) \delta_{\alpha_1' \alpha_2}\delta_{\alpha_2'\alpha_1} \,.\end{aligned}$$ Here, the composite indices $k_n$ denote pairs $({{{\bf{k}}}}_n,\omega_n)$ of momentum and Matsubara frequency, while the spin index $\alpha$ is written out explicitly. The basis function $U(k_1',k_2';k_1,k_2)$ is symmetric under simultaneous exchange of ingoing and outgoing indices.

For further simplification of the vertex parametrization we resort to the momentum space patching approximation outlined in Ref. [@Halboth2000], which is suitable in the weak coupling limit. In this approximation, the frequency dependence of the vertex is neglected, while the momentum dependence is parametrized such that it resolves the angular dependence around the Fermi surface, but it neglects any dependence in the radial direction. This is achieved by partitioning the Brillouin zone into a set of *patches* $\{\mathcal{P}_1 \dots \mathcal{P}_N\}$ as shown in Fig. \[fig:appendix:patching\], and projecting all momentum points within a patch $\mathcal{P}_m$ onto a single representative point ${{{\bf{p}}}}_m$ on the Fermi surface, i.e. the vertex function is assumed to be constant within the entire patch. The parametrization of the vertex function can thus be written as $$\begin{aligned}
&U(k_1',k_2';k_1,k_2) = \sum\limits_{i_1, i_2, i_3} u(n_{i_1}, n_{i_2}, n_{i_3}) \nonumber\\
&\quad\times \delta({{{\bf{k}}}}_1'+{{{\bf{k}}}}_2'-{{{\bf{k}}}}_1-{{{\bf{k}}}}_2) \delta_{{{{\bf{k}}}}_1' \in \mathcal{P}_{i_1}} \delta_{{{{\bf{k}}}}_2' \in \mathcal{P}_{i_2}} \delta_{{{{\bf{k}}}}_1 \in \mathcal{P}_{i_3}} \,,\end{aligned}$$ where the indices $n_m$ enumerate momentum patches and the symbol $\delta_{q\in \mathcal{P}_n}=1$ if momentum $q$ lies within patch $\mathcal{P}_n$ and zero otherwise.

The FRG flow equations are obtained by introducing an additional dependence of the interaction vertex on some RG cutoff. We follow the temperature flow RG scheme outlined in Ref. [@Honerkamp2001], where the temperature itself assumes to role of the RG cutoff and the flow equations take the form $$\frac{\mathrm{d}}{\mathrm{d}T}u_T(n_1, n_2, n_3) = \mathcal{T}_{PP,T} + \mathcal{T}_{PH,T}^d + \mathcal{T}_{PH,T}^c \,,$$ where the three interaction channels (particle-particle, direct particle-hole, and crossed particle-hole interaction, respectively) are given by (terms in the same order as shown in Fig. \[fig:appendix:FRGchannels\])

$$\begin{aligned}
    \mathcal{T}_{PP,T}(n_1, n_2, n_3) &= -\sum\limits_{n} u_T(n_1, n_2, n) u_T(n, -n+n_1+n_2, n_3) L^+_T(n, n_1+n_2) \nonumber\\
    \mathcal{T}^d_{PH,T}(n_1, n_2, n_3) &= -\sum\limits_{n} \Big( 
    -2u_T(n_1,n,n_3) u_T(n+n_1-n_3,n_2,n) \nonumber\\ 
    &\quad\quad\quad\quad\quad +u_T(n_1,n,n+n_1-n_3) u_T(n+n_1-n_3,n_2,n) \nonumber\\ 
    &\quad\quad\quad\quad\quad +u_T(n_1,n,n_3) u_T(n_2,n+n_1-n_3,n) 
    \Big) L^-_T(n, n_1-n_3) \nonumber\\
    \mathcal{T}^c_{PH,T}(n_1, n_2, n_3) &= -\sum\limits_{n} u_T(n_1,n+n_2-n_3,n) u_T(n,n_2,n_3) L^-_T(n, n_2-n_3) \,.
    \end{aligned}$$

The internal propagator bubble is defined as $$L^\pm_T(n,m) = \int_{{{{\bf{k}}}}\in\mathcal{P}_n} \mp \frac{\lambda\left(\xi({{{\bf{k}}}})\right)\pm \lambda\left(\xi(\mp {{{\bf{k}}}}+ {{{\bf{p}}}}_m)\right)}{\xi({{{\bf{k}}}})\pm\xi(\mp {{{\bf{k}}}}+ {{{\bf{p}}}}_m)} \,,$$ where $\xi({{{\bf{k}}}})$ is the dispersion of the noninteracting system and $\lambda(\xi)$ is the temperature derivative of the Fermi distribution function $\lambda(\xi) = \xi {{\mathrm{e}}}^{\xi/T}[T^2 \left( {{\mathrm{e}}}^{\xi/T}+1 \right)^2]^{-1}$.

In this form, the flow equations can be solved numerically to connect the high-temperature limit, in which the effective interaction vertex Eqn.  equals the bare interaction as defined by the Hamiltonian $\mathcal{H}_\mathrm{int}$, to the effective low-energy theory.

Hubbard model
-------------

![Benchmark of the characteristic temperature $T_\mathrm{min}$, as determined in FRG calculations (including only particle-particle scattering channel), with the mean-field critical temperature of the Hubbard model. The Hubbard limit $\ell=0$ shows excellent agreement for fillings $n\gtrsim 10^{-3}$. The generalized Gaussian potential ($\ell=1$) approaches the Hubbard limit in the very dilute regime. []{data-label="fig:appendix:hubbard:hubbardFRG"}](hubbardFRG.pdf){width="\linewidth"}

The FRG flow equations derived in Sec. \[sec:appendix:frg\] are suited for the weak-coupling limit [@Honerkamp2001]. In the dilute limit, however, when the Fermi energy scale becomes small compared to the interaction potential, the weak-coupling scenario may be violated. In order to convince ourselves that the approach produces meaningful results nevertheless, we benchmark the implementation against the mean-field solution of the Hubbard model at small densities.

To this end, we consider again the general Gaussian potential introduced in the main article and set $\ell = 0$, while fixing the prefactor to $g_0=\frac{3}{2\pi}$. As displayed in Fig. \[fig:appendix:hubbard:hubbardFRG\], the characteristic temperature scale $T_\mathrm{min}$ obtained from the FRG solution is in excellent agreement with the critical temperature as determined by the mean-field approach. Only at extremely low densities, below fillings relevant for our studies of geometric domes of $T_c$, deviations manifest. The FRG results remain consistent when a finite $\ell=1$ is considered in the sense that the result smoothly connects to the Hubbard limit in the dilute limit where $k_F \ell \ll 1$. This is to be expected since the width of the Gaussian profile becomes large compared to the size of the Fermi surface and the interaction potential effectively appears almost constant.

[^1]: Properly speaking, we set $V({{\bf{q}}})$ to be a periodic Gaussian, given by $V({{\bf{q}}})=  -g_0(2\pi\ell^2)^{d/2} \sum_{{{\bf{G}}}}{{\mathrm{e}}}^{-|{{\bf{q}}}+{{{\bf{G}}}}|^2\ell^2/2}$, where ${{{\bf{G}}}}$ are reciprocal lattice vectors.

[^2]: We define $k_F^\star = (2\pi n^\star)^{1/2}$ in 2D, and $k_F^\star = (3\pi^2 n^\star)^{1/3}$ in 3D

[^3]: For the single band Hubbard model, this choice reproduces the superconducting $T_c$ as shown in the SM
