#  Exact theory of superconductivity from repulsion in narrow band systems

Valentin Crépel, Liang Fu  Massachusetts Institute of Technology, 77
Massachusetts Avenue, Cambridge, MA, USA

###### Abstract

We introduce an analytically-controlled theory of superconductivity from
repulsive interaction in narrow band systems, where pairing is induced by
virtual composite excitations at high energy. A quantitative formula relating
$T_{c}$ to interaction strength and single-particle band gap is obtained, with
maximum $T_{c}$ reaching about $10\%$ of the Fermi temperature. Our theory
demonstrates a reliable and robust mechanism for strong coupling
superconductivity, which enables unconventional pairing symmetry and filling-
dependent gap structure.

Superconductivity in conventional metals results from an effective attraction
between electrons mediated by the exchange of phonons [1, 2]. While this
attraction is much weaker than the bare Coulomb repulsion [3], the latter is
drastically renormalized downward by retardation effects [4]. Thanks to the
vast difference between Fermi and Debye energy, the phonon-mediated attraction
can overscreen the Coulomb repulsion to enable electron pairing and
superconductivity [5]. On the other hand, this crucial retardation condition
fails in systems with narrow bands or low carrier density. Yet
superconductivity has been found in a growing number of materials in such
strong-coupling regime. Two famous examples are (1) strontium titanate, the
most dilute bulk superconductor with Fermi energy as small as $1$meV [6]; (2)
magic-angle graphene with a record-low density $n_{\rm 2D}\sim
10^{11}$\text{\,}{\mathrm{cm}}^{-2}$$ and a very small bandwidth $\sim 10$meV
[7, 8, 9, 10]. Remarkably, the ratio of superconducting transition temperature
$T_{c}$ and Fermi temperature $E_{F}/k_{B}$ far exceeds typical values,
reaching as high as 0.01 in strontium titanate [6] and 0.1 in magic-angle
graphene [8]. Finding electronic mechanisms for strong-coupling
superconductivity in narrow band systems has long been a subject of great
interest and challenge [11, 12, 13, 14].

In this work, we introduce a simple two-band model of interacting spin-
polarized fermions on a two-dimensional lattice to establish a general
mechanism for unconventional superconductivity from strong repulsion. In our
model, an insulating state occurs at the filling of $n=1$ fermion per unit
cell, and superconductivity emerges upon particle or hole doping. Based on a
perturbative expansion around the atomic limit, we rigorously show that,
despite the strong bare repulsion, a non-retarded short-range pairing
interaction between doped fermions in the narrow band is generated by coupling
to high-energy composite excitations. The resulting superconductor is
unconventional by all standards. It has $f$-wave pairing symmetry and changes
from having a full gap to point nodes above a critical doping. $T_{c}$ is
controlled by the bare interaction strength and the band gap, reaching as
large as $T_{c}\sim 0.1E_{F}/k_{B}$ when they are of comparable magnitude. Our
theory is analytically controlled by a small coupling constant that emerges in
the narrow band limit. Our work demonstrates a reliable and robust mechanism
for unconventional superconductivity from repulsion in narrow band systems.

Our model is broadly inspired by a recent work by Slagle and one of us [15],
who proposed a mechanism for pairing from purely classical electrostatic
repulsion in a doped charge transfer insulator. The essential ingredient there
is a charge-$2e$ excitation dubbed “trimer”, which is a composite object
consisting of two doped electrons tightly bound to a dipole. Interestingly,
for certain extended Coulomb repulsion, a trimer is energetically more favored
than two separate electrons. Under appropriate conditions, the presence of
preformed trimers can lead to Wigner crystal or superconducting ground states
at small doping.

While we also start with the problem of doping an insulator, our work differs
fundamentally from Ref. [15]. We find superconductivity without invoking
trimers or any preformed pairs at low energy. Instead, pairing arises from
correlated quantum hopping of doped fermions induced by virtual composite
excitations at high energy. We demonstrate this novel mechanism by introducing
and solving the simplest model of spin-polarized fermions interacting on a
bipartite lattice. Last but not the least, our solution reveals distinct
superconducting states at different ranges of doping and provides a
quantitative formula for $T_{c}$ in terms of microscopic parameters.

We consider spin-polarized fermions on the honeycomb lattice with repulsive
interactions, described by the Hamiltonian

$\displaystyle\mathcal{H}$ $\displaystyle=$
$\displaystyle\mathcal{H}_{0}+\mathcal{H}_{t},$ $\displaystyle\mathcal{H}_{0}$
$\displaystyle=$ $\displaystyle V\sum_{\langle
r,r^{\prime}\rangle}n_{r}n_{r^{\prime}}+\Delta\sum_{r\in B}n_{r},$
$\displaystyle\mathcal{H}_{t}$ $\displaystyle=$ $\displaystyle-t\sum_{\langle
r,r^{\prime}\rangle}(c_{r}^{\dagger}c_{r^{\prime}}+{hc}).$ (1)

where $\mathcal{H}_{0}$ contains the nearest-neighbor interaction – the
dominant interaction for spin-polarized fermions on a lattice – and the
sublattice potential difference between the two inequivalent $A$ and $B$
sites, while $\mathcal{H}_{t}$ describes tunneling between adjacent sites.
Despite the simplicity of our model, in this work we unveil its remarkably
rich phase diagram as a function of filling and interaction strength. We shall
derive the low-energy properties of the system with a fully-controlled
perturbative expansion in the narrow-band limit $t\ll\Delta$, which we further
complement with field-theoretic analysis and extensive exact diagonalization
(ED) studies.

Figure 1:  Low energy fermion added above the $n=1$ insulating background live
on the $B$ lattice. Excitations above this $f$-band are dipoles, polarons and
trimers (solid circles), whose virtual occupation leads to the effective model
Eq. (2) (dashed circles).

In the strong coupling limit, the ground state of $\mathcal{H}$ at $n=1$ is an
insulator with all $A$ sites occupied and all $B$ sites empty, as shown in
Fig. 1. Its insulating property is ensured by the large gap $E_{D}=2V+\Delta$,
which corresponds to the energy necssecary to transfer an electron from an $A$
site to $B$, or equivalently, creating a dipole.The inclusion of tunneling,
small compared to $E_{D}$, slightly decreases the charge transfer gap without
any significant change to the insulating ground state.

Since $\mathcal{H}$ is invariant under particle-hole transformation
$c_{A}\rightarrow c_{A}^{\dagger},c_{B}\rightarrow-c_{B}^{\dagger}$ combined
with spatial inversion that interchanges the two sublattices, it suffices to
consider $n>1$ filling below. At finite doping $n=1+\delta$ ($\delta>0$), low
energy configurations of the system remain with all $A$ sites occupied in
order to avoid the large charge transfer gap. Due to Pauli exclusion
principle, the $\delta$ additional fermions must live on the $B$ lattice, and
in the limit $t=0$, form a highly degenerate manifold with an energy
$E_{f}=\Delta+3V$ per doped charge that we refer to as $f$-band.

Besides these fermions on $B$ sites, there exist various types of composite
excitations at higher energy, which involve holes on $A$ sites, as depicted in
Fig. 1. For example, a $B$-fermion can bind with a neighboring dipole to form
a charge-$e$ Fermi polaron, which has energy $E_{P}=E_{f}+V+\Delta$. More
interesting is the charge-$2e$ trimer, which consists of three neighboring
$B$-fermions surrounding a hole on the center $A$ site. It can also be viewed
as two neighboring $B$-fermions tightly bound to a dipole. A trimer costs
energy $E_{T}=2E_{f}+\Delta$, which is greater than the energy of two separate
$B$-fermions by $\Delta$. These composite excitations – dipoles, polarons and
trimers – are hereafter collectively referred to as charge-transfer complex.

In the presence of small quantum tunneling $t\ll\Delta$, doped carriers in the
$f$-band constitute the only low-energy excitations in our system. They
virtually couple to charge-transfer complex at high energy. This coupling
results in a narrow dispersive $f$-band of doped carriers and induces short-
range interactions between them. Remarkably, we shall show that the induced
interaction leads to pairing within the $f$-band. To that purpose, we
analytically carry out a Schrieffer-Wolff transformation
$\mathcal{H}^{\prime}=e^{iS}\mathcal{H}e^{-iS}$ to decouple the $f$-band from
high energy degrees of freedom [16, 17]. As detailed in App. A, this procedure
accounts for all possible virtual processes (see Fig. 1), and leads to the
following effective Hamiltonian for doped fermions, which is exact to second
order in $t/\Delta$ and at any $f$-band filling:

$\begin{split}\mathcal{H}^{\prime}=&\sum_{\langle
i,j\rangle}t_{f}(f_{i}^{\dagger}f_{j}+hc)+V_{f}n_{i}n_{j}\\\
&+\sum_{(ijk)\in\triangle}\lambda(f_{i}^{\dagger}n_{j}f_{k}+P_{ijk})+U_{3}n_{i}n_{j}n_{k}.\end{split}$
(2)

The $f_{i}$ fermionic operators denote the doped fermions on the triangular
$B$-lattice, and their vacuum is the $n=1$ insulating state described above.
The sums labeled by $\langle i,j\rangle$ and $(i,j,k)\in\triangle$
respectively run over all bonds and all upper triangles of the $B$-lattice,
while $P_{ijk}$ stands for the inclusion of $f_{i}^{\dagger}n_{j}f_{k}$ with
all possible permutations of the indices $i$, $j$ and $k$.

The effective Hamiltonian $\mathcal{H}^{\prime}$ for doped fermions consists
of single-particle tunneling, correlated (density-dependent) tunneling, two-
body and three-body density interactions. Their origins can be understood as
follows. The tunneling from $k$ to $i$ in the upper triangle $(ijk)$ arises
from two consecutive hopping processes. The virtual intermediate state
involved is either a polaron or a trimer, depending on the occupation of site
$j$ (see Fig. 1). The resulting tunneling amplitude is thus
$t^{2}[(1-n_{j})/(E_{P}-E_{f})+n_{j}/(E_{T}-2E_{f})]$, from which the
expression of $t_{f}$ and $\lambda$ are derived:

$\displaystyle t_{f}$ $\displaystyle=\frac{t^{2}}{\Delta+V}\,,$ (3a)
$\displaystyle\lambda$
$\displaystyle=\frac{t^{2}}{\Delta}-\frac{t^{2}}{\Delta+V}\,.$ (3b) The
interaction coefficients $V_{f}$ and $U_{3}$ come from processes where an $A$
fermion hops back and forth between neighboring sites (see Fig. 1). For
example, $V_{f}$ measures the difference of energy between two neighboring
doped charges on $B$ sites and two well-separated ones. The former
configuration can couple to a trimer state whereas the latter cannot, thus
leading to an attraction $-t^{2}/\Delta$. Accounting for all processes, we
find $\displaystyle V_{f}$
$\displaystyle=-\frac{t^{2}}{\Delta}+\frac{4t^{2}}{\Delta+V}-\frac{3t^{2}}{\Delta+2V}\,,$
(3c) $\displaystyle U_{3}$
$\displaystyle=\frac{3t^{2}}{\Delta}-\frac{6t^{2}}{\Delta+V}+\frac{3t^{2}}{\Delta+2V}\,.$
(3d)

In Eq. (3a-d), the denominators $\Delta$, $\Delta+V$ and $\Delta+2V$ are the
energy costs of intermediate states involving trimer, polaron and dipole
respectively. Importantly, higher-order corrections to the effective
Hamiltonian $\mathcal{H}^{\prime}$ are small provided that the narrow band
condition $t\ll\Delta$ is satisfied, regardless of interaction strength $V$.

At small $V\ll\Delta$, the effective interactions in the $f$-band are found to
be $V^{0}_{f}=2(t/\Delta)^{2}V,\lambda^{0}=(t/\Delta)^{2}V$ and $U^{0}_{3}=0$
to first order in $V/\Delta$. These values simply correspond to the projection
of the bare repulsion $V$ into the $f$-band, whose wavefunctions have small
amplitudes $\sim t/\Delta$ on $A$ sites. As $V$ increases, interband mixing
quickly becomes important and our exact results (3a-d) reveals a dramatic
departure of the “dressed” interaction from the projected interaction. As
opposed to the projected interactions, $V_{f}$ starts to decrease at $V\approx
0.29\Delta$ and changes sign from repulsive to attractive at $V=\Delta$.
Similarly, $\lambda$ and $U_{3}$ show manifest deviations from the projected
estimates for $V>0.1\Delta$, and saturate at large $V/\Delta$ (see App. A). In
the rest of this work, we shall mainly consider the case $V<\Delta$, where the
induced $V_{f},\lambda,U_{3}$ turn out to be positive and small compared to
the single-particle bandwidth $W=9t_{f}$.

To reveal the tendency towards pairing, we study the formation of two-particle
bound states, the analog of “Cooper problem” in a doped insulator. Pair
formation is evidenced by the positivity of the pair binding energy

$\displaystyle\varepsilon_{b}$ $\displaystyle=$ $\displaystyle
2[E(1)-E(0)]-[E(2)-E(0)]$ (4) $\displaystyle=$ $\displaystyle
2E(1)-E(2)-E(0),$

with $E(m)$ the ground state energy of a system with $m$ charges added above
$n=1$ filling. As detailed in App. B, we analytically solve the lattice
Hamiltonian Eq. (2) in the case of two doped fermions and obtain
$\varepsilon_{b}$ as a function of $V/\Delta$, shown in Fig. 2. It is found
positive in the entire range $V/\Delta$. This result shows an effective
pairing interaction between low-energy fermions. It is worth noting that the
two-particle bound states cannot be captured by the projected interaction
$V_{f}^{0}=2\lambda^{0}$, which proves that pairing is induced by virtual
interband excitations. We further confirm the formation of pairs in the
original model Eq. ( Exact theory of superconductivity from repulsion in
narrow band systems) with ED, see App. E. Interestingly, pairs are already
present for $V<\Delta$ despite a repulsive induced nearest neighbor
interaction ($V_{f}>0$). This highlights the essential role of correlated
hopping $\lambda$ in the effective model for pairing.

Figure 2:  a) The transition from weakly to tightly bound pairs, depicted for
$V/\Delta=0.5$ (left) and $10$ (right), constitutes a probe of the BCS-BEC
crossover at low doping. b) Pair binding energy $\varepsilon_{b}$ and bound
state size as a function of $V/\Delta$ from the exact solution of the two
fermion problem. The prediction of the continuum model (dashed) perfectly
match the full-fledged lattice calculation, a stringent test of our
derivation.

The two-particle bound statee we found has zero total momentum, is symmetric
under three-fold rotation, and changes sign under reflection that flips one of
the primitive vectors ${\bf a}_{j}$, i.e. it has $f$-wave pairing symmetry.
The size of the bound state shrinks with increasing $V/\Delta$, as shown in
the inset of Fig. 2 for $V/\Delta=0.5$ and $100$. At small $V/\Delta$, the
pair wavefunction is highly extended over many lattice sites, while in the
opposite limit $V/\Delta\rightarrow\infty$, two nearest-neighbor fermions form
the most tightly bound pair “resonating” within a single upper triangle (see
App. B). Our solution of the two-particle problem suggests a crossover between
BCS and BEC states at small particle density, tuned by $V/\Delta$.

We now extend our analysis to finite, but small, doping concentrations
$\delta$ above $n=1$. In this regime, the physics of dilute doped fermions is
governed by long-wavelength properties that transcend the details on the
lattice scale. This motivates us to derive a low-energy theory by taking the
continuum limit of the lattice Hamiltonian. This is achieved by rewriting the
lattice Hamiltonian in terms of fermionic fields in momentum space and
retaining only modes near the bottom of the $f$-band. Importantly, the band
dispersion $\varepsilon({\bf k})=2t_{f}\sum_{j=1}^{3}\cos({\bf k}\cdot{\bf
a}_{j})$, with $t_{f}>0$, has two degenerate minima located at the $\pm\bf K$
points of the Brillouin zone. Therefore, low-energy degrees of freedom are
described by two long-wavelength fermionic fields $\psi_{\tau}({\bf
q})=f(\tau{\bf K}+{\bf q})$ with $qa\ll 1$, distinguished by the valley index
$\tau=\pm$. In App. C, we find the following continuum Hamiltonian for
$\psi_{\tau}$:

$\widetilde{\mathcal{H}}=\int{\rm
d}x\,\sum_{\tau=\pm}\psi_{\tau}^{\dagger}\left[\frac{-\nabla^{2}}{2m}\right]\psi_{\tau}+g\psi_{+}^{\dagger}\psi_{+}\psi^{\dagger}_{-}\psi_{-},$
(5)

where the effective mass and interaction strength are entirely determined from
the lattice parameters:

$\displaystyle m$ $\displaystyle=$ $\displaystyle 2/(3t_{f}a^{2}),$
$\displaystyle g$ $\displaystyle=$ $\displaystyle 6a^{2}(V_{f}-2\lambda)<0.$
(6)

The resulting quantum field theory describes a two-flavor fermion gas in the
continuum with _attractive_ contact interaction. The two flavors correspond to
the valley degree of freedom associated with the underlying lattice, from
which the field theory is derived. This theory is asymptotically exact in the
low doping limit where $s$-wave scattering between fermions of opposite
valleys is the dominant interaction.

This attractive interaction leads to the formation valley-singlet two-particle
bound states, which exactly correspond to the $f$-wave pairs observed on the
lattice (see Fig. 2). Indeed, the pair amplitude $f_{+K}f_{-K}$ is odd under
the reflection that interchanges the two valleys. To verify the validity of
our continuum model, we calculate the two-particle binding energy
$\varepsilon_{b}$ in the field theory, and using the parameters $m$ and $g$
given by Eq. ( Exact theory of superconductivity from repulsion in narrow band
systems) and (3), compare it with the exact solution of the lattice model. The
expression for $\varepsilon_{b}$ is

$\displaystyle\frac{\varepsilon_{b}}{\varepsilon_{\rm uv}}$ $\displaystyle=$
$\displaystyle\left[e^{1/g_{0}}-1\right]^{-1},$ $\displaystyle\quad g_{0}$
$\displaystyle=$ $\displaystyle\frac{9}{\pi}\frac{2\lambda-
V_{f}}{W}=\frac{6}{\pi}\frac{V^{2}}{\Delta(\Delta+2V)}\,,$ (7)

with $\varepsilon_{\rm uv}=\pi W/9$ an energy cutoff that we fix with the
exact binding energy at $V\to\infty$ (see App. C). The exponent $g_{0}$ in Eq.
( Exact theory of superconductivity from repulsion in narrow band systems),
defined by the ratio of effective pairing interaction $g$ and the bandwidth,
only depends on the ratio $V/\Delta$ in the narrow band regime $t\ll\Delta$.
Remarkably, a perfect agreement between continuum theory and lattice model is
found at all values of $V/\Delta$ (see Fig. 2). This proves the accuracy of
our mapping from lattice model to continuum theory.

Figure 3: Critical temperature of the continuum model in the exact BCS (solid,
Eq. (8)) and BEC (dashed) limits, as a function of doping and $V/\Delta$. A
gray horizontal line highlights the bound $k_{B}T_{c}<E_{F}/8$. The self-
consistent lattice mean-field solutions (dots) agree with the continuum theory
for small doping, but differ when the $\pm K$ pockets merge (inset).

The attractive Fermi gas in two dimensions is known to be superconducting at
low temperature and exhibits a BCS-BEC crossover as the ratio between pair
binding and Fermi energies changes from small to large values [18, 19, 20]. In
the region of weakly bound pairs $\varepsilon_{b}\ll E_{F}$, the critical
temperature is given by
$k_{B}T_{c}=e^{\gamma-1}\sqrt{2E_{F}\varepsilon_{b}}/\pi$ with $\gamma\simeq
0.577$ Euler’s constant [21, 22, 23]. From the expression for binding energy
Eq. ( Exact theory of superconductivity from repulsion in narrow band
systems), we derive a formula for $T_{c}$ in our model:

$k_{B}T_{c}=e^{\gamma-1}\sqrt{\frac{2E_{F}W}{9\pi}}e^{-1/(2g_{0})}\,.$ (8)

On the other side of the crossover $\varepsilon_{b}\gg E_{F}$, the physics
depends on the interaction between the bosonic pairs. When these bosons repel,
the system exhibits a BKT transition toward a BEC at low temperature [24],
while it collapses if bosons attract [25, 26]. Between the extreme BCS and BEC
limits, the critical temperature satisfies the very general bound
$k_{B}T_{c}\leq E_{F}/8$ [27], which limits the largest achievable $T_{c}$.

The BCS-BEC crossover of the two-dimension Fermi gas can be achieved by tuning
either carrier density or the interaction strength $g_{0}$, which is
controlled by $V/\Delta$ in our model. We plot in Fig. 3 the critical
temperature $T_{c}$ as a function of doping concentration $\delta$ and
$V/\Delta$. At very low doping where $E_{F}<\varepsilon_{b}$, the system lies
in the BEC regime and $T_{c}$ increases rapidly with $\delta$. At some
critical concentration, the system undergoes the BEC-BCS crossover and finally
follows Eq. (8).

From the expression of $g_{0}$ (Eq. ( Exact theory of superconductivity from
repulsion in narrow band systems)), we infer the critical temperature as a
function of $V/\Delta$. It is worth emphasizing that our exact BCS formula Eq.
(8) applies provided that the dimensionless coupling constant $g_{0}$ is
small, even when the bare repulsion $V$ far exceeds the bandwidth $W$. This is
because doped fermions at the conduction band bottom $\pm K$ reside entirely
on $B$ sublattice, and therefore avoid the direct nearest-neighbor repulsion
$V$. In the weak-coupling regime $V\ll\Delta$, the attraction
$g_{0}\propto(V/\Delta)^{2}$ between low-energy carriers is induced by virtual
interband particle-hole pairs or excitons, and leads to exponentially small
$T_{c}$. Most importantly however, our expression of $g_{0}$ is
nonperturbative in $V$ and remains exact at $V\sim\Delta\gg W$, where strong-
coupling superconductivity and maximum $T_{c}$ are attained. For instance, at
doping $\delta=0.1$, $T_{c}$ reaches $0.1E_{F}\simeq 0.032t^{2}/\Delta$ around
$V=0.43\Delta$, which is about $0.5\%$ of the quasiparticle bandwidth
$W=9t_{f}\simeq 6.3t^{2}/\Delta$. More generally, Eq. (8) can safely be
applied for $V<0.7\Delta$, since $\exp(1/g_{0})>10$ and
$\varepsilon_{b}\lesssim E_{F}/2$ [19].

Figure 4:  a) The ground state energy exhibits the characteristic flux
dependence of a superconductor, as shown here for $V=\Delta$ and 6 particles
on a $8\times 6$ lattice. b )The charge stiffness divided by the Fermi energy
is constant in the BCS limit, reaching 1/8 for small doping (gray). b)
Superconducting gap as a function of $V/\Delta$ for different doping
concentrations. It follows the continuum prediction (dashed) up to
$V/\Delta\simeq 0.7$.

We support the emergence of superconductivity with evidence from ED on finite
size lattices. While the original model Eq. ( Exact theory of
superconductivity from repulsion in narrow band systems) exhibits positive
$\varepsilon_{b}$ (see App. E), we focus on the effective model Eq. (2) which
allows to reach larger system sizes. First, the superfluid behavior of the
system is probed by the charge stiffness [28]

$D=\frac{1}{16\pi^{2}}\left.\frac{L_{1}}{L_{2}}\frac{\partial^{2}E(N,\phi)}{\partial\phi^{2}}\right|_{\phi=0}$
(9)

at doping $\delta=N/(L_{1}\times L_{2})$, with $N$ the number of doped
fermions and $(L_{1},L_{2})$ the number of sites along the two basis vectors
of the triangular lattice. $D$ measures the sensitivity of the ground state
energy $E(N,\phi)$ to twisted boundary conditions
$\psi_{r+L_{1}}=e^{2i\pi\phi}\psi_{r}$. A positive value of $D>0$ in the
thermodynamic limit implies dissipationless charge transport and gives a
direct signature of the Meissner effect [29]. In the range of parameters
considered, our system clearly exhibits (with small finite-size effect) the
$h/2e$ flux periodicity of superconductors [30] and shows positive $D$ (see
Fig. 4a-b), which proves superconductivity in the ground state. The charge
stiffness of a BCS superconductor with a parabolic dispersion relation is
known exactly: $D=E_{F}/4\pi$ [27]. Our results for $V<\Delta$ are correctly
captured by this prediction, especially at low doping concentrations where
fermions live close to the band minima.

To demonstrate the strong-coupling nature of the superconducting state, we
consider the superconducting gap

$\Delta_{\rm gap}=\frac{(-1)^{N}}{2}\left[E(N+1)+E(N-1)-2E(N)\right]\,.$ (10)

Our continuum theory Eq. (5) predicts $\Delta_{\rm
gap}=\sqrt{2E_{F}\varepsilon_{b}}$ up to $V\sim\Delta$. This leads to a ratio
of the gap and critical temperature $\Delta_{\rm gap}/k_{B}T_{c}=\pi
e^{1-\gamma}\simeq 4.796$, which is much larger than the universal value
$1.764$ in BCS theory for weak-coupling superconductors [31]. This is because
the phonon-induced retarded attraction in conventional metals is limited to
electrons within a Debye energy from the Fermi surface, whereas in our theory
the induced pairing interaction is instantaneous on the time scale of inverse
bandwidth ($\hbar/W$), so that all carriers in the narrow band are subject to
the pairing interaction.

Our numerical results for $\delta\simeq 0.1$, shown in Fig. 4b, agree with the
mean field prediction up to $V/\Delta=0.7$, where the gap extracted from ED
reaches as large as $\Delta_{\rm gap}=0.84E_{F}$. This allows the system to
reach critical temperature of 0.1$E_{F}$, as described above. For
$V<0.3\Delta$, the finite lattice size in ED study forbids the bound state to
fully form (see Fig. 2), resulting in an overestimate of the superconducting
gap compared to the continuum theory prediction.

In addition to the BCS-BEC superconductivity at low density, our model shows
very rich physics at higher doping concentrations, where lattice effects
become important. By performing mean-field calculation on the model Eq. (2),
we find $f$-wave pairing for doping $\delta<1/3$ (see App. F). The
corresponding critical temperatures, shown in Fig. 3, are calculated with the
linearized gap equation

$\frac{1}{\alpha}=\frac{1}{N_{s}}\sum_{{\bf q}}\frac{\left[\sum_{j}\sin({\bf
q}\cdot{\bf a}_{j})\right]^{2}}{|\xi_{{\bf k}}|}\tanh\left(\frac{|\xi_{{\bf
k}}|}{2k_{B}T_{c}}\right)$ (11)

with $\xi_{{\bf k}}=\varepsilon({\bf k})-\mu$ and $\alpha=2\lambda-
V_{f}-(2\beta+\delta)U_{3}$, where $\beta=(3N_{s})^{-1}\sum_{{\bf q}}f_{\rm
FD}(\xi_{{\bf k}})\sum_{j}\cos({\bf q}\cdot{\bf a}_{j})$ originates from
three-body interactions. Here, $N_{s}$ denotes the total number of sites, and
the chemical potential is fixed by $\delta=N_{s}^{-1}\sum_{{\bf q}}f_{\rm
FD}(\xi_{{\bf k}})$, with $f_{\rm FD}$ the Fermi-Dirac distribution.

At low doping, the $f$-wave superconducting state has a full pairing gap and
its $T_{c}$ obtained from lattice model calculation agrees well with our
previous result based on continuum theory. Interestingly, at higher doping
concentration $\delta>\delta_{c}\approx 1/4$, the gap vanishes at 6 nodes on
the Fermi surface along $\Gamma M$ direction, where the $f$-wave gap function
vanishes. This change in gap structure is due to the change of Fermi surface
topology across the van Hove singularity, where the two pockets around $\pm K$
merge into a single Fermi surface enclosing the $\Gamma$ point (see Fig. 3).

The above conclusions are confirmed by our ED study, which shows clear
evidence of nodal superconductivity at $\frac{1}{4}<\delta\leq\frac{1}{3}$. We
also find that the superconducting state is remarkably robust against longer
range bare repulsion. These numerical results can respectively be found in
App. G and D. Finally, the ground state at sufficiently large $V/\Delta$ is
found to be non-superconducting (see App. H), which we leave for further study
in a separate work.

Our work opens a new route to unconventional superconductivity in atomic Fermi
gas and electron systems. Encouragingly, optical lattice with honeycomb
geometry and tunable band gap have already been realized [32, 33]. Many recent
advances in dipolar or Rydberg atom systems with longer-range interactions
[34, 35] have enabled the implementation of one-dimensional $t-V$ Hamiltonian
[36], and hold great promise for the realization of our model in near future.

Our mechanism for strong-coupling superconductivity mediated by charge-
transfer complex, or interband excitations, may also shed insight on graphene-
based moiré superlattices, where the small bandwidth and high
$k_{B}T_{c}/E_{F}$ ratio (up to $\sim 0.1$) place important constrains on
viable theories. It will be interesting to develop accurate low-energy models
for these systems and analyze superconductivity in the narrow band limit as
exemplified in our work. In this regard, we note that correlated hopping and
direct repulsion also appear in effective Hamiltonian for narrow bands in
twisted bilayer graphene [37, 38, 39]. Our theory suggests that
renormalization by virtual interband exitations is necessary to obtain strong-
coupling superconductivity. Moreover, our simple model may be relevant to
twisted double bilayer graphene [40] and trilayer graphene-boron nitride
heterostructures [41], where signs of spin-polarized superconductivity has
been reported [42, 43]. We leave to future work the extension of our theory to
spinful systems and its application to various strongly-correlated materials.

Acknowledgments. This work was supported by DOE Office of Basic Energy
Sciences, Division of Materials Sciences and Engineering under Award DE-
SC0019275. LF was supported in part by a Simons Investigator Award from the
Simons Foundation.

## References

  * Bardeen  _et al._ [1957] John Bardeen, Leon N Cooper, and John Robert Schrieffer, “Theory of superconductivity,” Physical review 108, 1175 (1957).
  * Fröhlich [1950] Herbert Fröhlich, “Theory of the superconducting state. i. the ground state at the absolute zero of temperature,” Physical Review 79, 845 (1950).
  * Bardeen and Pines [1955] John Bardeen and David Pines, “Electron-phonon interaction in metals,” Physical Review 99, 1140 (1955).
  * Morel and Anderson [1962] P Morel and PW Anderson, “Calculation of the superconducting state parameters with retarded electron-phonon interaction,” Physical Review 125, 1263 (1962).
  * Bogoljubov  _et al._ [1958] NN Bogoljubov, Vladimir Veniaminovic Tolmachov, and DV Širkov, “A new method in the theory of superconductivity,” Fortschritte der physik 6, 605–682 (1958).
  * Lin  _et al._ [2013] Xiao Lin, Zengwei Zhu, Benoît Fauqué, and Kamran Behnia, “Fermi surface of the most dilute superconductor,” Phys. Rev. X 3, 021002 (2013).
  * Cao  _et al._ [2018] Yuan Cao, Valla Fatemi, Shiang Fang, Kenji Watanabe, Takashi Taniguchi, Efthimios Kaxiras, and Pablo Jarillo-Herrero, “Unconventional superconductivity in magic-angle graphene superlattices,” Nature 556, 43–50 (2018).
  * Lu  _et al._ [2019] Xiaobo Lu, Petr Stepanov, Wei Yang, Ming Xie, Mohammed Ali Aamir, Ipsita Das, Carles Urgell, Kenji Watanabe, Takashi Taniguchi, Guangyu Zhang, _et al._ , “Superconductors, orbital magnets and correlated states in magic-angle bilayer graphene,” Nature 574, 653–657 (2019).
  * Park  _et al._ [2020] Jeong Min Park, Yuan Cao, Kenji Watanabe, Takashi Taniguchi, and Pablo Jarillo-Herrero, “Tunable phase boundaries and ultra-strong coupling superconductivity in mirror symmetric magic-angle trilayer graphene,” arXiv preprint arXiv:2012.01434 (2020).
  * Hao  _et al._ [2020] Zeyu Hao, AM Zimmerman, Patrick Ledwith, Eslam Khalaf, Danial Haie Najafabadi, Kenji Watanabe, Takashi Taniguchi, Ashvin Vishwanath, and Philip Kim, “Electric field tunable unconventional superconductivity in alternating twist magic-angle trilayer graphene,” arXiv preprint arXiv:2012.02773 (2020).
  * Ruhman and Lee [2016] Jonathan Ruhman and Patrick A Lee, “Superconductivity at very low density: The case of strontium titanate,” Physical Review B 94, 224515 (2016).
  * Qin  _et al._ [2020] Mingpu Qin, Chia-Min Chung, Hao Shi, Ettore Vitali, Claudius Hubig, Ulrich Schollwöck, Steven R White, Shiwei Zhang, _et al._ , “Absence of superconductivity in the pure two-dimensional hubbard model,” Physical Review X 10, 031016 (2020).
  * Raghu  _et al._ [2010] S Raghu, SA Kivelson, and DJ Scalapino, “Superconductivity in the repulsive hubbard model: An asymptotically exact weak-coupling solution,” Physical Review B 81, 224505 (2010).
  * Maiti and Chubukov [2013] Saurabh Maiti and Andrey V Chubukov, “Superconductivity from repulsive interaction,” in  _AIP Conference Proceedings_ , Vol. 1550 (American Institute of Physics, 2013) pp. 3–73.
  * Slagle and Fu [2020] Kevin Slagle and Liang Fu, “Charge transfer excitations, pair density waves, and superconductivity in moir$\backslash$’e materials,” arXiv preprint arXiv:2003.13690 (2020).
  * Schrieffer and Wolff [1966] John R Schrieffer and Peter A Wolff, “Relation between the anderson and kondo hamiltonians,” Physical Review 149, 491 (1966).
  * MacDonald  _et al._ [1988] Allan H MacDonald, SM Girvin, and D t Yoshioka, “t/u expansion for the hubbard model,” Physical Review B 37, 9753 (1988).
  * Randeria  _et al._ [1989] Mohit Randeria, Ji-Min Duan, and Lih-Yir Shieh, “Bound states, cooper pairing, and bose condensation in two dimensions,” Physical review letters 62, 981 (1989).
  * Bertaina and Giorgini [2011] Gianluca Bertaina and S Giorgini, “Bcs-bec crossover in a two-dimensional fermi gas,” Physical review letters 106, 110403 (2011).
  * Parish [2015] Meera M Parish, “The bcs–bec crossover,” in  _Quantum Gas Experiments: Exploring Many-Body States_ (World Scientific, 2015) pp. 179–197.
  * Miyake [1983] Kazumasa Miyake, “Fermi liquid theory of dilute submonolayer 3he on thin 4he ii film: dimer bound state and cooper pairs,” Progress of theoretical physics 69, 1794–1797 (1983).
  * Gor’kov and Melik-Barkhudarov [1961] LP Gor’kov and TK Melik-Barkhudarov, “Contribution to the theory of superfluidity in an imperfect fermi gas,” Sov. Phys. JETP 13, 1018 (1961).
  * Petrov  _et al._ [2003] DS Petrov, MA Baranov, and GV Shlyapnikov, “Superfluid transition in quasi-two-dimensional fermi gases,” Physical Review A 67, 031601 (2003).
  * Fisher and Hohenberg [1988] Daniel S Fisher and PC Hohenberg, “Dilute bose gas in two dimensions,” Physical Review B 37, 4936 (1988).
  * Ruprecht  _et al._ [1995] 1PA Ruprecht, MJ Holland, K Burnett, and Mark Edwards, “Time-dependent solution of the nonlinear schrödinger equation for bose-condensed trapped neutral atoms,” Physical Review A 51, 4704 (1995).
  * Roberts  _et al._ [2001] Jacob L Roberts, Neil R Claussen, Simon L Cornish, Elizabeth A Donley, Eric A Cornell, and Carl E Wieman, “Controlled collapse of a bose-einstein condensate,” Physical Review Letters 86, 4211 (2001).
  * Hazra  _et al._ [2019] Tamaghna Hazra, Nishchhal Verma, and Mohit Randeria, “Bounds on the superconducting transition temperature: Applications to twisted bilayer graphene and cold atoms,” Physical Review X 9, 031049 (2019).
  * Kohn [1964] Walter Kohn, “Theory of the insulating state,” Phys. Rev. 133, A171–A181 (1964).
  * Scalapino  _et al._ [1993] Douglas J Scalapino, Steven R White, and Shoucheng Zhang, “Insulator, metal, or superconductor: The criteria,” Physical Review B 47, 7995 (1993).
  * Loder  _et al._ [2008] F. Loder, A. P. Kampf, and T. Kopp, “Crossover from $hc/e$ to $hc/2e$ current oscillations in rings of $s$-wave superconductors,” Phys. Rev. B 78, 174526 (2008).
  * Schrieffer [2018] J.R. Schrieffer,  _Theory Of Superconductivity_ (CRC Press, 2018).
  * Uehlinger  _et al._ [2013] Thomas Uehlinger, Gregor Jotzu, Michael Messer, Daniel Greif, Walter Hofstetter, Ulf Bissbort, and Tilman Esslinger, “Artificial graphene with tunable interactions,” Physical review letters 111, 185307 (2013).
  * Fläschner  _et al._ [2016] Nick Fläschner, BS Rem, M Tarnowski, D Vogel, D-S Lühmann, K Sengstock, and Christof Weitenberg, “Experimental reconstruction of the berry curvature in a floquet bloch band,” Science 352, 1091–1094 (2016).
  * Baier  _et al._ [2018] S Baier, D Petter, JH Becher, A Patscheider, G Natale, L Chomaz, MJ Mark, and F Ferlaino, “Realization of a strongly interacting fermi gas of dipolar atoms,” Physical review letters 121, 093602 (2018).
  * Lu  _et al._ [2012] Mingwu Lu, Nathaniel Q Burdick, and Benjamin L Lev, “Quantum degenerate dipolar fermi gas,” Physical Review Letters 108, 215301 (2012).
  * Guardado-Sanchez  _et al._ [2020] Elmer Guardado-Sanchez, Benjamin Spar, Peter Schauss, Ron Belyansky, Jeremy T Young, Przemyslaw Bienias, Alexey V Gorshkov, Thomas Iadecola, and Waseem S Bakr, “Quench dynamics of a fermi gas with strong long-range interactions,” arXiv preprint arXiv:2010.05871 (2020).
  * Guinea and Walet [2018] Francisco Guinea and Niels R Walet, “Electrostatic effects, band distortions, and superconductivity in twisted graphene bilayers,” Proceedings of the National Academy of Sciences 115, 13174–13179 (2018).
  * Kang and Vafek [2019] Jian Kang and Oskar Vafek, “Strong coupling phases of partially filled twisted bilayer graphene narrow bands,” Phys. Rev. Lett. 122, 246401 (2019).
  * Koshino  _et al._ [2018] Mikito Koshino, Noah F. Q. Yuan, Takashi Koretsune, Masayuki Ochi, Kazuhiko Kuroki, and Liang Fu, “Maximally localized wannier orbitals and the extended hubbard model for twisted bilayer graphene,” Phys. Rev. X 8, 031087 (2018).
  * Liu  _et al._ [2020] Xiaomeng Liu, Zeyu Hao, Eslam Khalaf, Jong Yeon Lee, Yuval Ronen, Hyobin Yoo, Danial Haei Najafabadi, Kenji Watanabe, Takashi Taniguchi, Ashvin Vishwanath, _et al._ , “Tunable spin-polarized correlated states in twisted double bilayer graphene,” Nature 583, 221–225 (2020).
  * Chen  _et al._ [2019] Guorui Chen, Aaron L Sharpe, Patrick Gallagher, Ilan T Rosen, Eli J Fox, Lili Jiang, Bosai Lyu, Hongyuan Li, Kenji Watanabe, Takashi Taniguchi, _et al._ , “Signatures of tunable superconductivity in a trilayer graphene moiré superlattice,” Nature 572, 215–219 (2019).
  * Lee  _et al._ [2019] Jong Yeon Lee, Eslam Khalaf, Shang Liu, Xiaomeng Liu, Zeyu Hao, Philip Kim, and Ashvin Vishwanath, “Theory of correlated insulating behaviour and spin-triplet superconductivity in twisted double bilayer graphene,” Nature communications 10, 1–10 (2019).
  * Cornfeld  _et al._ [2020] Eyal Cornfeld, Mark S Rudner, and Erez Berg, “Spin-polarized superconductivity: order parameter topology, current dissipation, and multiple-period josephson effect,” arXiv e-prints , arXiv–2006 (2020).
  * Vidal  _et al._ [2000] Julien Vidal, Benoit Douçot, Rémy Mosseri, and Patrick Butaud, “Interaction induced delocalization for two particles in a periodic potential,” Physical review letters 85, 3906 (2000).
  * Levinsen and Parish [2015] Jesper Levinsen and Meera M Parish, “Strongly interacting two-dimensional fermi gases,” in  _Annual review of cold atoms and molecules_ (World Scientific, 2015) pp. 1–75.

##  Appendix A Schrieffer-Wolff Transformation

In this appendix, we derive effective model Eq. 2 with a Schrieffer-Wolff
transformation. This effective Hamiltonian is exact up to corrections of order
$(t/\Delta)^{2}$, irrespective of the ratio $t/V$.

###  A.1 Canonical Transformation

Single particle tunneling on the honeycomb lattice couples the $f$-band to
trimers, polarons and dipoles. The Schrieffer-Wolff transformation uses a
canonical transformation $\mathcal{H}^{\prime}=e^{iS}\mathcal{H}e^{-iS}$, with
$S$ Hermitian, to treat these couplings as an effective Hamiltonian that
leaves the $f$-band invariant up to second order correction in $t/\Delta$.
This is achieved if $S$ satisfies

$[\mathcal{H}_{0},iS]=\mathcal{H}_{t}\,,$ (12)

as can be seen with the Baker-Campbell-Haussdorf formula. Under the assumption
Eq. 12, we get

$\mathcal{H}^{\prime}=\mathcal{H}_{0}+\left[iS,\mathcal{H}_{t}\right]+\mathcal{O}(\mathcal{H}_{t}S^{2})\,.$
(13)

To find $S$, we follow Ref. [17] and split the tunneling Hamiltonian into a
collection of operators $T_{p,M}$

$\mathcal{H}_{t}=\sum_{p=\pm 1}\sum_{M=-2}^{2}T_{p,M}\,,$ (14)

where $T_{p,M}$ gathers all tunneling operations that change the number of
occupied $A$ sites by $p$ and the number of occupied nearest neighbor pairs by
$M$. The bounds for $p$ and $M$ are determined by the lattice geometry. With
these definition, we can check that

$S=-i\sum_{p,M}\frac{T_{p,M}}{MV+p\Delta}\,,$ (15)

satisfies Eq. 12. Consider an eigenstate $|n\rangle$ of $\mathcal{H}_{0}$ with
energy $E_{n}$. By definition $T_{p,M}|n\rangle$ is also an eigenstate of
$\mathcal{H}_{0}$ with energy $E_{n}+(p\Delta+MV)$. As a consequence, we have

$\displaystyle[\mathcal{H}_{0},iS]|n\rangle$
$\displaystyle=\sum_{p,M}\frac{(E_{n}+p\Delta+VM)T_{p,m}-T_{p,M}E_{n}}{p\Delta+VM}|n\rangle$
$\displaystyle=\sum_{p,M}T_{p,M}|n\rangle=\mathcal{H}_{t}|n\rangle\,.$ (16)

Finally, the relation $T_{p,M}^{\dagger}=T_{-p,-M}$ ensures the hermiticity of
$S$. Plugging Eq. 15 into Eq. 13, we find the generic expression

$\mathcal{H}^{\prime}=\mathcal{H}_{0}+\sum_{p,p^{\prime},M,M^{\prime}}\frac{[T_{p^{\prime},M^{\prime}},T_{p,M}]}{M^{\prime}V+p^{\prime}\Delta}+\mathcal{O}\left(\frac{t^{3}}{\Delta^{2}}\right)\,.$
(17)

To evaluate the correcting terms of this last equation, we have noticed that
all terms in $S$ are smaller than $t/\Delta$ in magnitude.

###  A.2 Low energy projection

To describe the low-energy properties of our model, we project the obtained
Hamiltonian onto the $f$-band which is well separated from other excitations
by the charge transfer gap $\Delta\gg t$. This projection restricts the sum in
Eq. 17 to cases where $p^{\prime}=-p$ and $M^{\prime}=-M$. Furthermore, the
first operator acting on the states of the $f$-band should move an electron
from an $A$ to a $B$ site, i.e. the rightmost $T_{p,M}$ must have $p,M\geq 0$.
This gives

$\mathcal{H}^{\prime}\simeq\mathcal{H}_{0}-\sum_{M=0,1,2}\frac{T_{-1,-M}T_{1,M}}{MV+\Delta}\,.$
(18)

The three terms of the sum $M=0,1,2$ gather all second order processes that
respectively involve the virtual occupation of a trimer, a polaron and a
dipole excitation (shown in Fig. 1). As expected from second order
perturbation theory, these processes occur with rates

$t_{T}=\frac{t^{2}}{\Delta},\quad t_{P}=\frac{t^{2}}{V+\Delta},\quad\text{and
}\quad t_{D}=\frac{t^{2}}{2V+\Delta},$ (19)

that are inversely proportional to their energy difference with the $f$-band.

###  A.3 Simplification

To further simplify $\mathcal{H}^{\prime}$, we isolate the contributions of
Eq. 18 depending on the number of occupied neighbor of each $A$-site. In the
following, we focus on a $A$-site at position $r_{0}$ and denote its three
neighbors as $r_{1}$, $r_{2}$ and $r_{3}$.

#### No neighboring fermions:

If the $A$-site is surrounded by three empty sites, it can couples to three
different dipole excitations in which the fermion at $r_{0}$ is moved on a $B$
neighboring site. This provides an energy shift

$\mathcal{H}_{0}=-3t_{D}(1-n_{r_{1}})(1-n_{r_{1}})(1-n_{r_{1}})\,.$ (20)

#### One neighboring fermion:

If the $A$-site has a neighboring fermion at position $r_{1}$, it can couples
to two polarons by hopping to either $r_{2}$ or $r_{3}$. If the second
tunneling process in Eq. 18 moves it back to $r_{0}$, the process leads to an
energy shift. On the contrary, if the fermion at $r_{1}$ replace the original
one at $r_{0}$, the process can be viewed as an effective tunneling of the
$r_{1}$-fermion on a neighboring $B$-site. Together, they give

$\begin{split}\mathcal{H}_{1}=-&2t_{P}n_{r_{1}}(1-n_{r_{2}})(1-n_{r_{3}})\\\
&+t_{P}[f_{r_{2}}^{\dagger}(1-n_{r_{3}})+f_{r_{3}}^{\dagger}(1-n_{r_{2}})]f_{r_{1}}\,.\end{split}$
(21)

Note that the tunneling terms come with a positive sign due to the
anticommutation relation between fermions. Circular permutation of the indices
$(r_{1},r_{2},r_{3})$ gives all other possible terms involving a polaron.

#### Two neighboring fermions:

Finally, if $r_{1}$ and $r_{2}$ both host a fermion, the system can only
couple to a trimer excitation by moving the $A$ electron to $r_{3}$. The
second tunneling process can either put this fermion back to $r_{0}$, leading
to en energy reduction, or put the ones at $r_{1}$ or $r_{2}$ at $r_{0}$,
giving a density-assisted tunneling term:

$\mathcal{H}_{1}=-t_{T}n_{r_{1}}n_{r_{2}}(1-n_{r_{3}})+t_{T}(f_{r_{2}}^{\dagger}n_{r_{3}}+f_{r_{3}}^{\dagger}n_{r_{2}})f_{r_{1}}\,.$
(22)

Again, circular permutation of the indices $(r_{1},r_{2},r_{3})$ gives all
other possible terms involving a trimer.

Figure S1:  Coefficients of the effective model Eq. 2 (full lines), and
projection of the bare repulsion into the $f$-band (dashed) as a function of
$V/\Delta$.

Gathering the contributions of a-b-c, we obtain the effective Hamiltonian Eq.
2 where the coefficients originates from the processes outlined in Fig. 1. The
single particle fermion tunneling only appears in Eq. 21 and therefore reads
$t_{f}=t_{P}$. The correlated tunneling appears in Eqs. 22 and 21 with
opposite signs: $\lambda=t_{T}-t_{P}$. Two-body and three-body interaction
coefficients come from the three cases a-b-c. Carefully counting all terms, we
find $V_{f}=-t_{T}+4t_{P}-3t_{D}$ and $U_{3}=3t_{T}-6t_{P}+3t_{D}$. These
expression are given in Eq. 3, and plotted against $V/\Delta$ in Fig. S1. For
comparison (see text), we also show the results from the projection of the
bare repulsion $V$ into the $f$-band, valid when $V\ll\Delta$.

##  Appendix B Solution of the two body problem

In this appendix, we solve the effective lattice model Eq. 2 for two fermions.
To understand the competing roles of two-body interaction and correlated
hopping, we first look at the quartic part of the Hamiltonian (Sec. B.1).
Then, we decouple the center of mass and relative motions of the two fermions
(Sec. B.2) and explain in more detail how the results presented in the main
text were obtained (Sec. B.3).

###  B.1 ’Resonating’ triangle

We start with the quadratic terms of our model

$\mathcal{H}_{\rm int}^{\prime}=V_{f}\sum_{\langle
i,j\rangle}n_{i}n_{j}+\lambda\sum_{(ijk)\in\triangle}(f_{i}^{\dagger}n_{j}f_{k}+P_{ijk})\,.$
(23)

All states where the two electrons are not neighbors are zero energy
eigenstates of this operator. The others have two doped charge on the same
upper triangle, and can be divided into irreducible representations of
$C_{3v}$. More precisely, there is a state in the identity representation
$A_{1}$ and an $E$-doublet

$\displaystyle|A_{1},{\bf r}\rangle=\frac{1}{\sqrt{3}}\left[f_{{\bf
r}}^{\dagger}f_{{\bf r}+{\bf a}_{1}}^{\dagger}+f_{{\bf r}+{\bf
a}_{1}}^{\dagger}f_{{\bf r}-{\bf a}_{3}}^{\dagger}+f_{{\bf r}-{\bf
a}_{3}}^{\dagger}f_{r}^{\dagger}\right]|n=1\rangle$ $\displaystyle|E,{\bf
r}\rangle=\frac{1}{\sqrt{6}}\left[f_{{\bf r}}^{\dagger}f_{{\bf r}+{\bf
a}_{1}}^{\dagger}+f_{{\bf r}+{\bf a}_{1}}^{\dagger}f_{{\bf r}-{\bf
a}_{3}}^{\dagger}-2f_{{\bf r}-{\bf
a}_{3}}^{\dagger}f_{r}^{\dagger}\right]|n=1\rangle$
$\displaystyle|E^{\prime},{\bf r}\rangle=\frac{1}{\sqrt{2}}\left[f_{{\bf
r}}^{\dagger}f_{{\bf r}+{\bf a}_{1}}^{\dagger}-f_{{\bf r}+{\bf
a}_{1}}^{\dagger}f_{{\bf r}-{\bf a}_{3}}^{\dagger}\right]|n=1\rangle\,,$ (24)

with $|n=1\rangle$ the insulating state at unit filling, which also
corresponds to the vacuum for the $f$ operators. These three eigenstates of
$\mathcal{H}_{\rm int}^{\prime}$ have energy

$E_{A_{1}}=V_{f}-2\lambda\,,\quad E_{E}=V_{f}+\lambda\,.$ (25)

Because $V_{f}-2\lambda<0$, the lowest energy manifold is made of all states
$\\{|A_{1},{\bf r}\rangle\\}_{r\in B}$ and has _negative_ energy,
corresponding to bound states. As mentioned in the main text, this binding
energy solely comes from correlated hopping when $V\leq\Delta$, for which
$V_{f}>0$. These dimers have the $f$-wave symmetry described in the main text
(see Fig. 2), and can be seen as dimers ’resonating’ within an upper triangle
as shown in Fig. S2.

Figure S2:  a) The lowest energy subspace of $\mathcal{H}_{\rm int}^{\prime}$
is composed of bound states ($E_{A_{1}}=V_{f}-2\lambda<0$), which can be seen
as dimers ’resonating’ within an upper triangle of the $B$-lattice. b) Weight
of the two-body ground state $\Gamma(q)$ in the limit $V\gg\Delta$. It is
maximal in amplitude and opposite in sign at the $K$ and $K^{\prime}$ points,
underlying its valley-singlet nature.

This solution unveils the crucial role of correlated tunneling for pairing and
gives an intuitive real-space picture for the bound pairs. It also allows to
better understand the low energy continuum model Eq. 5. First, the coupling
constant $g$ is directly proportional to the energy $E_{A_{1}}<0$. Moreover,
it also gives a physical intuition for the $s$-wave scattering between the two
valleys. To see this, let’s reintroduce the tunneling $t_{f}$ in the lowest
energy subspace $\\{|A_{1},{\bf r}\rangle\\}_{r\in B}$ to lift the degeneracy
of the low energy manifold (this perturbative treatment of $t_{f}$ is
justified in the limit $V\gg\Delta$). The dispersion relation of the bound
pairs $\varepsilon_{A_{1}}(k)=E_{A_{1}}-(2t_{f}/3)\sum_{j=1}^{3}\cos\left({\bf
k}\cdot{\bf a}_{j}\right)$ has its minimum at the $\Gamma$ point. The lowest
energy state is

$|\Gamma^{2e}\rangle=\frac{i}{\sqrt{3N_{s}}}\sum_{{\bf q}\in BZ}\Gamma({\bf
q})f_{-{\bf q}}^{\dagger}f_{{\bf q}}^{\dagger}|n=1\rangle\,,$ (26)

with $\Gamma({\bf q})=\sum_{j=1}^{3}\sin({\bf q}\cdot{\bf a}_{j})$. The
amplitude of the paired state $\Gamma(q)$ is maximal in magnitude near the $K$
and $K^{\prime}$ points, as shown in Fig. S2. This confirms our intuition that
the dominant scattering channel couples electrons with different isospin
$\tau=\pm$ in $s$-wave to form valley singlets. While this argument gives an
intuitive understanding of the dominant pairing interaction in our system, it
can only be formally applied in the limit where $V\gg\Delta$ where the
tunneling $t_{f}$ can be treated as a perturbation to $\mathcal{H}_{\rm
int}^{\prime}$.

###  B.2 Center of mass and relative motion

To solve the two-particle problem exactly, we now decouple the center of mass
and relative motion of the two fermions. Taking all terms of Eq. 2
simultaneously into account, we obtain the results presented in Fig. 2, which
also agree with the previous perturbative treatment of the previous section in
the limit $V\gg\Delta$.

The Hilbert space with two fermions is spanned by the states

$|r_{1},r_{2}\rangle=f_{r_{1}}^{\dagger}f_{r_{2}}^{\dagger}|n=1\rangle\,.$
(27)

Taking advantage of the translation invariance of the problem, we can
introduce the center of mass momentum $K$ and reorganize the Hilbert space
with the Bloch-waves [44]

$|\varphi(K,r)\rangle=\frac{1}{\sqrt{N_{s}}}\sum_{R}e^{i(K\cdot
R)}|R,R+r\rangle\,.$ (28)

Because $|r_{2},r_{1}\rangle=-|r_{1},r_{2}\rangle$, the state with opposite
relative positions $r$ describe the same physical state
$|\varphi(K,-r)\rangle=-e^{i(K\cdot r)}|\varphi(K,r)\rangle$. To avoid double
counting the states, we restrict our attentions to states with
$(r\cdot\delta_{1})\geq 0$. The action of the Hamiltonian in this new basis
can be directly computed

$\displaystyle\mathcal{H}^{\prime}|\varphi(K,r)\rangle=$ (29)
$\displaystyle\sum_{\begin{subarray}{c}j=1,2,3\\\
\epsilon=\pm\end{subarray}}t_{f}[1+e^{i\epsilon(K\cdot
a_{j})}]|\varphi(K,r+\epsilon a_{j})\rangle+V_{f}\delta_{r,\epsilon
a_{j}}|\varphi(K,r)\rangle$
$\displaystyle+\lambda\sum_{j}\delta_{r,a_{j}}|\varphi(K,-a_{j-1})\rangle+\delta_{r,-a_{j-1}}|\varphi(K,a_{j})\rangle$
$\displaystyle+\lambda\sum_{j}e^{-i(K\cdot
a_{j+1})}\delta_{r,-a_{j}}|\varphi(K,a_{j-1})\rangle$
$\displaystyle+\lambda\sum_{j}e^{i(K\cdot
a_{j-1})}\delta_{r,a_{j}}|\varphi(K,-a_{j+1})\rangle\,.$

As expected from translation invariance, the state $|\varphi(K,r)\rangle$ only
couple to states with the same center of mass momentum $K$. The tunneling on
the lattice acts as an effective hopping term for the relative position $r$.
The local interaction and correlated tunneling terms only changes this tight-
binding Hamiltonian close to the origin for $r=a_{j}$ or $r=-a_{j}$ with
$j=1,2,3$. This tight binding Hamiltonian can be solved numerically for each
center of mass momentum $K$. Its ground state always is at the $\Gamma$ point,
and reads $|\psi_{\rm bd}\rangle=\sum_{r}\psi_{\rm
bd}(r)|\varphi(0,r)\rangle$. Its energy yields the binding energy and bound
state size

$\xi_{\rm bd}=\sqrt{\sum_{r}r^{2}|\psi_{\rm bd}(r)|^{2}}\,,$ (30)

in Fig. 2. Of course, accurate results can only be obtained when the system
size exceeds the spread of the two-fermion bound state. In particular for
$V/\Delta<0.25$, the bound state size sharply increases and a more careful
extraction of the energy is required.

###  B.3 Finite size scaling

Figure S3:  Bindinng energy as a function of the inverse number of site. The
linear behavior is extrapolated to $1/N_{s}\to 0$ for $V/\Delta<0.25$ where
extremely large system sizes are needed to accommodate the very weakly bound
states (see Fig. 2).

We therefore focus on the regime $V/\Delta<0.25$, and extract the two-body
problem as a function of the inverse number of site in the system $1/N_{s}$ as
shown in Fig. S3. Linear interapolation of these curves at $1/N_{s}=0$ gives
the values of the binding energy in the thermodynamic limit. The interpolated
values are reported in Fig. 2. It is worth pointing out that we infer positive
pair binding energy in the entire range $V/\Delta>0$, in agreement with the
continuum model.

##  Appendix C Continuum model

###  C.1 Derivation

In this appendix, we derive the continuum model Eq. 5 in the limit of small
doping concentration. Focusing on the physics near the $K$ and $K^{\prime}$
pockets, and discarding all high energy fermionic degrees of freedom amounts
to write the Fourier transform of the $f$-operators as

$f_{r}=\frac{1}{\sqrt{N_{s}}}\sum_{\tau=\pm}\sum_{k,\,ka\ll 1}e^{-i[(\tau
K+k)\cdot r]}\psi_{\tau,k}\,,$ (31)

with $N_{s}$ the number of honeycomb unit cells. Total momentum conservation
restricts the interaction terms to valley-preserving ones. They can be further
decomposed as intra-valley interactions

$H_{p}^{\prime}=\frac{1}{2N_{s}}\sum_{k,q,q^{\prime}}\left[\sum_{\tau}V_{k,q+\tau
K}\psi_{\tau,q-k}^{\dagger}\psi_{\tau,q^{\prime}+k}^{\dagger}\psi_{\tau,q^{\prime}}\psi_{\tau,q}\right]$
(32)

and inter-valley interactions

$H_{s}^{\prime}=\frac{1}{2N_{s}}\sum_{k,q,q^{\prime}}\left(V_{k,q+K}+V_{-k,q^{\prime}-K}-V_{q^{\prime}-q+k-2K,q^{\prime}-K}-V_{q-q^{\prime}-k+2K,q+K}\right)\psi_{+,q-k}^{\dagger}\psi_{-,q^{\prime}+k}^{\dagger}\psi_{-,q^{\prime}}\psi_{+,q}\,,$
(33)

which both depend on the momentum space representation of the quartic terms of
Eq. 2:

$V_{k,q}=\frac{V_{f}}{2t_{f}}\varepsilon(k)+\lambda\sum_{j=1}^{3}e^{i(q\cdot
a_{j}+k\cdot a_{j+1})}+e^{-i(q\cdot a_{j}+k\cdot a_{j-1})}.$ (34)

Due to the small momenta considered in our continuum theory, we can Taylor
expand the previous expression with the help of the relations

$V_{k,q+\tau K}=3(V_{f}-\lambda),\,V_{k+2\tau K,q+\tau
K}=-\frac{3}{2}V_{f}+6\lambda,$ (35)

which hold true up to $\mathcal{O}(k,q)$ corrections. This leading order
approximation and the fermion anti-commutation relations show that intra-
valley interactions $H_{p}^{\prime}\simeq 0$ are negligible compared to inter-
valley ones

$H_{s}^{\prime}=\frac{V_{0}}{2N_{s}}\sum_{k,q,q^{\prime}}\psi_{+,q-k}^{\dagger}\psi_{-,q^{\prime}+k}^{\dagger}\psi_{-,q^{\prime}}\psi_{+,q}\,,$
(36)

with $V_{0}=9(V_{f}-2\lambda)=-9U_{3}$. Fourier transforming back to real-
space and accounting for the Brillouin zone area, we obtain the effective
interaction announced in Eq. 5:

$\widetilde{H}_{\rm int}=g\int{\rm
d}x\,\psi_{+}^{\dagger}\psi_{-}^{\dagger}\psi_{-}\psi_{+}\,,\quad
g=2V_{0}a^{2}/3\,.$ (37)

###  C.2 Binding Energy

As a consistency check of our continuum model Eq. 5, we can determine the
binding energy and compare it with the exact solution found in App. B. Using a
$T$-matrix approach, the binding energy $\varepsilon_{b}$ of the continuum
theory is solution of the implicit equation [45]

$\frac{1}{g}=-\int\frac{{\rm
d}^{2}q}{(2\pi)^{2}}\frac{1}{\varepsilon_{b}+|q|^{2}/m}\,.$ (38)

Going to polar coordinates $q=q_{r}e^{i\theta}$ and introducing the momentum
UV cutoff $\Lambda$, we find

$\frac{2\pi}{m|g|}=\int_{0}^{\Lambda}{\rm
d}q_{r}\frac{q_{r}}{q_{r}^{2}+m\varepsilon_{b}}\,,$ (39)

which can be integrated to obtain

$\varepsilon_{b}=\frac{\Lambda^{2}}{m\left[e^{\frac{4\pi}{m|g|}}-1\right]}\,.$
(40)

In order to fix the cutoff $\Lambda$, we can use our solution of the problem
in the limit $V\gg\Delta$ of Sec. B.1, $\varepsilon_{b}=3t^{2}/\Delta$. In
that limit, Eq. 40 reduces to $\varepsilon_{b}\simeq(3\Lambda
t)^{2}/(2\pi\Delta)$. Equating the two limits yields:

$\Lambda=\sqrt{\frac{2\pi}{3a^{2}}}\simeq\frac{1.44720}{a}\,.$ (41)

In Eq.  Exact theory of superconductivity from repulsion in narrow band
systems, we have introduce the corresponding energy cutoff $\varepsilon_{\rm
uv}=\Lambda^{2}/m=\pi t_{f}=\pi W/9$.

##  Appendix D Longer range interactions

In this appendix, we show evidence that longer range interactions neither
destroy the effective pairing between doped charges nor the superconducting
phase described in the main text. We include next-nearest neighbor
interactions $\mathcal{H}_{2}=V_{2}\sum_{\langle\langle
r,r^{\prime}\rangle\rangle}n_{r}n_{r^{\prime}}$ in order to describe the tail
of the fermion-fermion interaction.

The range of next-nearest interaction strength relevant for our work is
$\Delta\sim V>V_{2}$. In that regime, the Schrieffer-Wolff transformation of
App. A almost identically applies, up to two corrections. To first order,
$\mathcal{H}_{2}$ introduces a direct interaction between electron in the
$f$-band. To second order in $t$, it also changes the denominators of Eq. 3,
which are determined by the new dipole, polaron and trimer energies:

$\begin{split}E_{D}&=\Delta+2V-6V_{2}\,,\\\ E_{P}&=E_{f}+\Delta+V-5V_{2}\,,\\\
E_{T}&=2E_{f}+\Delta-3V_{2}\,.\end{split}$ (42)

For weak perturbations $V_{2}\ll\Delta$, the second corrections can be safely
discarded, and we obtain the same effective model as Eq. 2 with a two-body
interaction strength

$V_{f}^{\prime}=V_{2}+V_{f}\,.$ (43)

At small doping, we can go to the continuum and the system is described by Eq.
5 with a modified coupling constant
$g^{\prime}=3a^{2}(V_{f}^{\prime}-2\lambda)$ (compare with Eq.  Exact theory
of superconductivity from repulsion in narrow band systems). When the
additional interaction term is smaller than the largest pair binding energy
$\varepsilon_{b}$ of the original problem, i.e. $V_{2}<3t^{2}/\Delta$, we find
a that $g^{\prime}$ is negative when $V/\Delta\geq x_{+}$, where

$x_{+}=\frac{3v_{2}+\sqrt{v_{2}(v_{2}+24)}}{4(3-v_{2})}\,,\quad
v_{2}=\frac{V_{2}\Delta}{t^{2}}\,.$ (44)

In other words, at small $\delta$, longer range interactions in the range
$V_{2}<3t^{2}/\Delta$ only shifts the range of $V/\Delta$ where
superconductivity is observed. For larger $V_{2}$, the bound pairs breaks to
avoid large next-nearest neighbor repulsion and the binding energy is always
negative.

Figure S4:  Charge stiffness as a function of $V_{2}$ for 4 (orange) and 6
(red) particles in a $7\times 7$ lattice with $V=0.6\Delta$.

Even in that regime, the superconducting phase can reappear at larger doping
concentrations. Indeed, if they are closely packed, the fermions of a pair do
not decrease their energy by splitting apart because this send them closer to
all the other surrounding fermions in the system. We evidence the stability of
the superconducting phase against $V_{2}$ at large doping with ED with 4 and 6
particles in a $7\times 7$ lattice, thus at doping $\delta=4/(7\times 7)\simeq
8\%$ and $6/(7\times 7)\simeq 12\%$, with $V=0.6\Delta$. The charge stiffness
extracted as a function of $V_{2}$ is shown in Fig. S4. Despite a slight
decrease, the stiffness remains positive up to $V_{2}=10t^{2}/\Delta$. This
gives support for the stability of the superconducting phase against longer
range interactions at doping $\delta\sim 0.1-0.15$.

Finally, the regime $\Delta\sim V_{2}\gg t$ has been extensively described in
the in Ref. [15]. In that case, the trimer energy $E_{T}=2E_{f}+\Delta-3V_{2}$
can be resonant with the $f$-band. The presence of these preformed pairs in
the Fermi sea of doped charges leads to superconductivity and pair-density
waves. While the origin of pairs is different in that regime, strong $V_{2}$
interactions can still lead to superconductivity.

##  Appendix E Pairing and superconductivity in the full-fledged model

Figure S5:  The pair binding energy and pair correlation function for finite
honeycomb lattice of 30 sites highlights the presence of two-body bound state
and of a superconducting phase in the original model.

We present ED results on the original model Eq.  Exact theory of
superconductivity from repulsion in narrow band systems in Fig. S5, which give
further evidence for the presence of two-body bound pairs at finite doping
above the $n=1$ insulating state, and for the emergence of superconductivity
in our model. The former is probed by the pair binding energy
$\varepsilon_{b}$ that we find, as in the main text, positive in the limit
$t\ll\Delta$ for the entire range of $V/\Delta$ considered. For $\Delta\sim
t$, we also observe a region with positive pair binding energy for
sufficiently large $V/\Delta$.

To investigate the superconducting correlations in the ground state of the
doped $n=1$ insulator, we calculate the $f$-wave pair correlation function

$P(r\in B)=\frac{1}{N_{s}}\sum_{x\in
B}\left\langle\Delta_{f}^{\dagger}(r+x)\Delta_{f}(x)\right\rangle\,,$ (45)

with $\Delta_{f}(x)=\sum_{j=1}^{3}c_{x+\delta_{j}}c_{x}$. A superconducting
state with $f$-wave symmetry is expected to exhibit a power-law decay of
$P(r)$ with distance, and hence to have sizable correlations at large
distances. For the largest accessible finite size clusters $N_{s}=5\times 3$,
the maximum distance is achieved for $r_{m}=(2,2)$. We denote the pair
correlation function at this distance as $P_{m}=P(r_{m})$. In Fig. S5, we show
$P_{m}$ as a function of the parameters $V/\Delta$ for several $\Delta/t$. We
observe large correlations in the perturbative limit $\Delta>5t$ when
$V/\Delta>0.25$, where we indeed expect a strong superconducting order (see,
for instance, Fig. 3). This direct manifestation of superconductivity gives
strong evidence for the reliability of the predictions made from the effective
model Eq. 2 and its mean-field solution detailed below in App. F.

##  Appendix F Mean-field treatment

Figure S6:  Numerical solution $\Delta_{k}$ of the self-consistent gap
equation for infinitesimal doping, decomposed into partial waves of $p$ and
$f$ symmetry. The inset shows the self-consistent estimate of the critical
temperature for $\delta=0.05$, which almost perfectly match the exact formula
Eq. 8.

###  F.1 Dilute limit

We now carry out a self-consistent mean field treatment of the low-energy
model Eq. 2 in the dilute limit, where three-body interaction terms are
negligible. Going to Fourier space, the Hamiltonian becomes

$\mathcal{H}^{\prime}=\sum_{k}\varepsilon_{k}f_{k}^{\dagger}f_{k}+\frac{1}{N_{s}}\sum_{k,q,q^{\prime}}V_{k,q}f_{q-k}^{\dagger}f_{q^{\prime}+k}^{\dagger}f_{q^{\prime}}f_{q}\,,$
(46)

with $V_{k,q}$ defined in Eq. 34. Then we perform the mean-field substitution
$f_{q^{\prime}}f_{q}\simeq\delta_{q+q^{\prime}}\langle f_{-q}f_{q}\rangle$.
This replacement leads to the following quadratic mean-field Hamiltonian

$\begin{split}\mathcal{H}_{\rm
MF}&=\frac{1}{2}\sum_{k}\begin{pmatrix}f_{k}^{\dagger}&f_{-k}\end{pmatrix}\begin{pmatrix}\xi_{k}&\Delta_{k}\\\
\Delta_{k}^{*}&-\xi_{k}\end{pmatrix}\begin{pmatrix}f_{k}\\\
f_{-k}^{\dagger}\end{pmatrix}\,,\\\ \Delta_{k}&=\frac{2}{N_{s}}\sum_{q}{\rm
Re}\,(V_{q-k,q})\langle f_{-q}f_{q}\rangle\,,\end{split}$ (47)

where we have introduced a chemical potential $\mu$ and the short notation
$\xi_{k}=\varepsilon_{k}-\mu$. We have also used the relation
$V_{-k,-q}=V_{k,q}^{*}$. The quadratic Hamiltonian $\mathcal{H}_{\rm MF}$ is
diagonalized as

$\mathcal{H}_{\rm MF}=\sum_{k}E_{k}\gamma_{k}^{\dagger}\gamma_{k}\,,\quad
E_{k}=\sqrt{\xi_{k}^{2}+|\Delta_{k}|^{2}}\,,$ (48)

through the Bogoliubov transformation

$\displaystyle\begin{pmatrix}f_{k}\\\ f_{-k}^{\dagger}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}u_{k}&v_{k}\\\
-v_{k}^{*}&u_{k}^{*}\end{pmatrix}\begin{pmatrix}\gamma_{k}\\\
\gamma_{-k}^{\dagger}\end{pmatrix}\,,$ (49) $\displaystyle u_{k}$
$\displaystyle=\sqrt{\frac{1}{2}\left(1+\frac{\xi_{k}}{E_{k}}\right)}\,,\,\,v_{k}=-\frac{\Delta_{k}}{|\Delta_{k}|}\sqrt{\frac{1}{2}\left(1-\frac{\xi_{k}}{E_{k}}\right)}\,.$

The ground state of the system is annihilated by all quasi-particle operators
$\gamma_{k}$ and takes the form $|\Psi_{\rm
BCS}\rangle=\prod_{k}\gamma_{k}|n=1\rangle$. The order parameter $\Delta_{k}$
may now be computed self-consistently through the calculation of the
correlators $\langle\Psi_{\rm BCS}|f_{-q}f_{q}|\Psi_{\rm BCS}\rangle$. This
gives the so-called gap equation

$\begin{split}\Delta_{k}&=\frac{1}{N_{s}}\sum_{q}{\rm
Re}\,(V_{q-k,q}-V_{q+k,q})u_{q}^{*}v_{q}\\\ &=\frac{1}{2N_{s}}\sum_{q}{\rm
Re}\,(V_{q+k,q}-V_{q-k,q})\frac{\Delta_{q}}{E_{q}}\,.\end{split}$ (50)

Using the explicit form of $V_{k,q}$, this gap equation can be rewritten as

$\Delta_{k}=\sum_{j=1}^{3}u_{j}\sin(k\cdot
a_{j})\,,\,\,u_{j}=\frac{1}{N_{s}}\sum_{q}\phi_{j}(q)\frac{\Delta_{q}}{E_{q}}\,,$
(51)

with

$\phi_{j}(q)=-V_{f}\sin(qa_{j})+\lambda[\sin(qa_{j+1})+\sin(qa_{j-1})]\,.$
(52)

Similarly, the mean particle number fixes the chemical potential $\mu$ through
the relation

$\delta=\frac{1}{2N_{s}}\sum_{q}\left(1-\frac{\xi_{q}}{E_{q}}\right)\,.$ (53)

We solve the self-consistent gap equation numerically for infinitesimal doping
and decompose the order parameter $\Delta_{k}$ into partial $p$ and $f$-wave.
The results are shown in Fig. S6. We observe that the self-consistent always
has a solution, with a substantial gap for $\Delta/V>0.25$. This provides
another test of the continuum results of Fig. 3, where the critical
temperature drastically rises near $\Delta/V=0.25$. Furthermore, our mean-
field solution exhibits a strong $f$-wave symmetry, which corroborates our
exact two-fermion calculation (Fig. 2).

###  F.2 Critical temperature and and three-body interactions

Including finite temperature effects to the previous analysis provides a way
to determine the critical temperature $T_{c}$ at which the pair are all broken
by thermal fluctuation and $\Delta_{k}=0$ [31]. This definition agrees with
the superconducting temperature in the BCS limit of our model. Note that it
does not capture the physics of the BKT transition, explaining why we obtain
$k_{B}T_{c}\geq E_{F}/8$ in Fig. 3. Focusing on an order parameter with pure
$f$-wave symmetry, the gap and number equations become

$\displaystyle 1$
$\displaystyle=\frac{1}{N_{s}}\sum_{q}\phi_{j}(q)\frac{\sum_{p}\sin(q\cdot
a_{p})}{|\xi_{q}|}\tanh\left(\frac{|\xi_{q}|}{2k_{B}T_{c}}\right)$ (54a)
$\displaystyle\delta$
$\displaystyle=\frac{1}{2N_{s}}\sum_{q}\left[1-\frac{\xi_{q}}{|\xi_{q}|}\tanh\left(\frac{|\xi_{q}|}{2k_{B}T_{c}}\right)\right]\,.$
(54b)

We solve this implicit definition of $T_{c}$ numerically. The results,
presented in the inset of Fig. S6 for $\delta=0.05$, show a sharp increase of
$T_{c}$ near $V=\Delta/4$. Excellent agreement is found with the continuum
results of Fig. 3 in the small doping limit. For $\delta\geq 0.1$, the non-
parabolicity of the dispersion relation leads to corrections to Eq. 8.

The three-body interaction terms of Eq. 2 can be partially reincorporated in
Eq. 54 near $T_{c}$ via the substitution

$n_{i}n_{j}n_{k}\simeq\delta(n_{i}n_{j}+n_{j}n_{k}+n_{i}n_{k})-\beta(f_{i}^{\dagger}n_{j}f_{k}+P_{ijk})\,,$

where we have replace one quadratic operator by its expectation in the normal
state $\langle n\rangle=\delta$ and $\langle
f_{k}^{\dagger}f_{i}\rangle=\beta=N_{s}^{-1}\sum_{k}\cos(k\cdot a_{1})f_{\rm
FD}(\xi_{k})$ with $f_{\rm FD}$ the Fermi-Dirac distribution. This simply
renormalizes the effective interaction $V_{f}$ and correlated hopping strength

$V_{f}^{\prime}=V_{f}+\delta U_{3}\,,\quad\lambda^{\prime}=\lambda-
U_{3}\beta\,,$ (55)

which leads to a similar replacement of $\phi_{j}$ by

$\phi_{j}^{\prime}(q)=-V_{f}^{\prime}\sin(qa_{j})+\lambda^{\prime}[\sin(qa_{j+1})+\sin(qa_{j-1})]$
(56)

in Eq. 54, as declared in the main text. The results shown in Fig. 3 are
obtained with this substitution.

For fillings $\delta>1/4$, we must rely on Eq. 54, which correctly describes
the merging of the $\pm K$ pockets, while the two-flavor continuum model Eq. 5
assume them well-separated. At these doping concentrations, our calculations
show that a superconducting state only appears for $V\geq 0.3\Delta$, in sharp
contrast with the dilute regime where superconductivity extends to
$V\ll\Delta$ (see Fig. 4). Another stark difference with the dilute limit is
the presence of nodes of the superconducting gap located at the Fermi surface,
which now encloses the $\Gamma$ point.

Figure S7:  Phase stiffness and superconducting gap obtained by ED of Eq. 8
for 6 (orange) and 8 (red) particles or more on $7\times 7$ finite lattices.

##  Appendix G ED evidence for nodal superconductor

We present more evidence of this nodal superconductor from ED of the effective
model for $\delta>1/4$. We focus on lattices of size $7\times 7$ with more
than 6 particles, and extract the charge stiffness and the superconducting gap
as in the main text. Our results are presented in Fig. S7. There, we observe
that the system is not superconducting for $V\ll\Delta$, in agreement with the
mean-field calculations of App. F. The transition to $D>0$ occurs around
$V=0.2\Delta$, consistent with the mean-field value up to finite-size jitters.
These observations comfort the presence of a superconducting phase at doping
$\delta>1/4$ and moderate $V/\Delta$.

In the entire superconducting phase evidenced above, the superconducting gap
$\Delta_{\rm gap}$ is almost zero and decrease in magnitude with larger system
size. This hints towards nodes of the superconducting gap at the Fermi energy,
as expected from the Fermi surface topology and $f$-wave symmetry of the order
parameter (see text).

##  Appendix H Bosonic effective interactions and phase separation

In this appendix, we present our numerical results obtained outside of the
superconducting region that we focus on in the main text. We have evaluated
the charges stiffness of the ground state in a wide region of doping
concentration and ratio $V/\Delta$. Our results are presented in Fig. S8 for
lattices with $L_{1}\geq 7$ and $L_{2}\geq 6$, which mitigate finite size
effects. A white dashed line indicates where $\varepsilon_{b}=E_{F}$, and
serves as an indicator to distinguish between the BCS and BEC limits predicted
by the continuum model Eq. 5. As shown in the main text, the stiffness in the
BCS region ($V<\Delta$) slowly increases but remains close to $E_{F}/8$,
especially at small doping concentrations. In the BEC limit, we observe two
different behaviors depending on the value of $\delta$. First, the results at
$\delta<10\%$ cannot reliably describe the thermodynamic behavior of the
system as they were obtain for less than four particles. In that case, the
small BEC is either a two-body or four-body bound state with a large effective
mass. For larger doping, the lattice effects observed in the main text rule
the physics and the charge stiffness drops down with increasing $V/\Delta$. We
observe that this diminution comes with a reduction of the may-body spectrum
gap, until the ground state becomes highly degenerate for the $V>70\Delta$.

To understand these features, we investigate the scattering properties of the
local bosons in the $V\gg\Delta$ limit. To distinguish between attractive and
repulsive interactions of pairs, we evaluate the second binding energy
$\varepsilon_{b}^{\prime}=2E(2)-E(4)-E(0)$, which plays the same role as
$\varepsilon_{b}$ for the tightly bound pairs. The results presented to the
right of Fig. S8, show a positive $\varepsilon_{b}^{\prime}$ for $V>\Delta$.
This signals attractive interaction, which leads to the collapse of the quasi-
BEC and to phase-separation. The collapse of the BEC in favor of phase
separation explains the large degeneracy observed for $V\gg\Delta$. The BEC
collapse at large doping concentrations for $V>\Delta$ is also evidenced by a
decay of the superconducting gap for $V>\Delta$, as shown in Fig. S9. Note
that the points in the range $V<1.2\Delta$ reproduce those of Fig. 4 in the
main text.

Figure S8:  Phase stiffness in a wider region of doping and ratio $V/\Delta$.
The white dashed line distinguishes between the BCS and BEC limits. The right
panel shows the second binding energy. Its positivity indicates effective
attractions between bosons deep in the BEC phase.  Figure S9:  Superconducting
gap in a wider range of $V/\Delta$. It reduces to almost zero at $V=4\Delta$
for large doping concentrations.