# The Effect of Repulsion on Superconductivity at Low Density

Dan Phan  Andrey V. Chubukov  School of Physics and Astronomy, University of
Minnesota, Minneapolis, Minnesota 55455, USA

###### Abstract

We examine the effect of repulsion on superconductivity in a three-dimensional
system with a Bardeen-Pines-like interaction in the low-density limit, where
the chemical potential $\mu$ is much smaller than the phonon frequency
$\omega_{L}$. We parameterize the strength of the repulsion by a dimensionless
parameter $f$, and find that the superconducting transition temperature
$T_{c}$ approaches a nonzero value in the $\mu=0$ limit as long as $f$ is
below a certain threshold $f^{*}$. In this limit, we find that $T_{c}$ goes to
zero as a power of $f^{*}-f$, in contrast to the high density limit, where
$T_{c}$ goes to zero exponentially quickly as $f$ approaches $f^{*}$. For all
nonzero $f$, the gap function $\Delta(\omega_{m})$ changes sign along the
Matsubara axis, which allows the system to partially overcome the repulsion at
high frequencies. We trace the position of the gap node with $f$ and show that
it approaches zero frequency as $f$ approaches $f^{*}$. To investigate the
robustness of our conclusions, we then go beyond the Bardeen-Pines model and
include full dynamical screening of the interaction, finding that $T_{c}$
still saturates to a non-zero value at $\mu=0$ when $f<f^{*}$.

##  I Introduction

In recent years, there has been a renewed interest in superconductivity at low
carrier density, arising from experimental advances in an assortment of
materials such as SrTiO3 [1, 2, 3], single-crystal Bi [4, 5], Pb1-xTlxTe [6],
Bi2Se3 [7], and half-Heusler compounds [8] (for recent review see Ref. [9]).
As of present, the origin of superconductivity in many of these low-density
materials is not well-understood. Due to the dilute nature of these compounds,
the repulsive electron-electron interaction is weakly screened, and one
naively expects this repulsion to dominate any attraction arising from the
electron-phonon interaction.

Non-$s$-wave superconductivity arising from electron-electron repulsion (e,g.,
attraction in the $d_{x^{2}-y^{2}}$ channel in the cuprates) is well
understood as the ultimate result of static screening of the Coulomb repulsion
in the particle-hole channel, which leads to Friedel oscillations at large
distances and in most cases generates an attraction in one or more
non-$s$-wave channels (see, e.g., Ref. 10 and references therein). However,
superconducting order parameters in SrTiO3 and other low-density materials are
likely $s$-wave, in which case a static interaction remains repulsive.

A standard formalism, which describes how $s$-wave superconductivity can exist
even for a (sufficiently weak) repulsive interaction, involves dynamical
screening [11, 12, 13, 9]. The reasoning is as follows [14, 15, 16, 17]: when
$\mu$ is much larger than the frequency of the pairing boson, $\omega_{L}$,
the bare Coulomb repulsion is renormalized down by dynamical screening in the
particle-particle channel, (the same channel which accounts for
superconductivity in the case of attraction), and at frequencies of order
$\omega_{L}$ is reduced by a factor $\ln(\mu/\omega_{L})$. If this reduced
Coulomb repulsion is smaller than the electron-phonon attraction, the
effective interaction is attractive at smaller frequencies, and
superconductivity develops. In a more accurate treatment [18, 19] the full
interaction (Coulomb + electron-phonon) remains repulsive, but is reduced at
small frequencies. Superconductivity then develops with a frequency-dependent
gap, which changes sign between small and large frequencies. An effective
description of a “conventional” sign-preserving $s$-wave superconductor with
an attractive interaction (electron-phonon minus Coulomb (reduced by
$\ln(\mu/\omega_{L})$)) emerges once one integrates out higher-energy
fermions, for which the sign of the gap is opposite to that at smaller
frequencies.

This reasoning holds in the high-density limit, where $\mu\gg\omega_{L}$, but
is not applicable to the case when the fermionic density is small and
$\mu\leq\omega_{L}$. There are two reasons for this: first, even for an
attractive interaction, Cooper pairing is thought to arise from fermions near
the Fermi surface, where the density of states can be approximated by its
value on the Fermi surface and the pairing kernel is logarithmically singular.
As $\mu$ decreases, the range where this description is applicable, shrinks.
Keeping contributions from the Cooper logarithm, one then finds that for
$\mu\leq\omega_{L}$, the prefactor of $T_{c}$ scales with $\mu$, leading to
the vanishing of $T_{c}$ as $\mu\to 0$. Second, by the same smallness of
$\mu/\omega_{L}$, the actual pairing interaction is almost independent of
frequency in the range where the density of states is approximated by its
value on the Fermi surface, and is repulsive. For such an interaction, there
is no solution of the gap equation as there is no way to obtain a sign change
in the gap function.

This first argument was re-analyzed by Gastiasoro et al. [20], who solved the
full momentum and frequency-dependent Eliashberg equations for $T_{c}$ in a
model of a three-dimensional electron gas with only the attractive, phonon-
mediated component of the pairing interaction,
$V(\Omega,\mathbf{q})=V(\mathbf{q})\omega^{2}_{L}/(\Omega^{2}+\omega^{2}_{L})$.
They argued that at small $\mu$, typical $q$ for the pairing are much larger
than $k_{F}$. As a result, the prefactor of $T_{c}$ does not vanish at
$\mu=0$. In fact, they found that $T_{c}$ actually increases in the $\mu=0$
limit and argued that this increase reflects the essentially unscreened nature
of the interaction. This result is consistent with the observations made
earlier by Takada and others [21, 15, 22, 5, 23] that even at moderately-large
density, there are features in the gap function away from $k=k_{F}$.

In this communication we analyze $T_{c}$ and the gap structure in the low-
density limit within the effective Bardeen-Pines model [24, 14, 20] with both
electron-phonon attraction and electron-electron repulsion. Specifically, we
investigate the pairing of electrons in three dimensions with a spherical
Fermi surface, interacting via

$V(\Omega,q)=\frac{4\pi
e^{2}}{q^{2}+\kappa^{2}}\big{(}f-\frac{\omega_{L}^{2}}{\Omega^{2}+\omega_{L}^{2}}\big{)}.$
(1)

Here $f$ is a measure of the strength of the repulsive interaction, and
$\kappa$ is the Thomas-Fermi wavevector. Gastiasoro et al. considered the case
$f=0$, while a physically-motivated interaction corresponds to $f\geq 1$. We
treat $f$ as a parameter and obtain results for values of $f$ both below and
above one.

The model of Eq. 1 excluding the momentum-dependence of the interaction has
been studied by a number of authors (see, e.g., Ref. 25 and references
therein). The result of these studies is that for a generic $\mu$,
superconductivity survives up to some critical $f^{*}(\mu)$. However, as $\mu$
approaches zero, $T_{c}$ vanishes for all $f>1$. Our goal is to study how
these results are modified if one takes the full momentum and frequency-
dependent $V(\Omega,q)$ in Eq. 1. To this end, we solve the momentum and
frequency-dependent integral equations for $T_{c}$ for $s$-wave pairing. Our
calculations show that $T_{c}$ does approach a finite value at $\mu=0$, even
in the presence of static repulsion. We argue that a non-zero $T_{c}$ results
from including both the sign change of the gap function between small and
large frequencies and the fact that at small $\mu$ the pairing involves
fermions away from the Fermi surface. As at high density, the pairing holds as
long as $f$ is below a certain $f^{*}(\mu)$. We show that at vanishing $\mu$,
$f^{*}(\mu)$ approaches a finite value $f^{*}(0)>1$. Unlike the high-density
case, where $T_{c}$ goes to zero exponentially quickly as $f$ approaches the
threshold value, we find that in the low-density limit $T_{c}$ goes to zero as
a power law: $T_{c}\propto(f^{*}-f)^{2}$.

The gap function $\Delta(\omega_{m})$ at $f<f^{*}$ changes sign at some
$\omega_{m}=\omega_{0}$. We show that $\omega_{0}$ exists for arbitrarily
small $f>0$, appearing at infinity when $f=0^{+}$. As $f$ increases,
$\omega_{0}$ decreases and ultimately vanishes at $f=f^{*}$. This result is
interesting from a topological perspective, since nodal points of
$\Delta(\omega_{m})$ correspond to centers of dynamical vortices [26, 27]. The
nodeless state at $f=0$ and a state with a nodal gap at $f>0$ have different
numbers of vortices and are therefore topologically distinct. In this respect,
the vanishing of superconductivity at $f=f^{*}$ can be viewed as a topological
transition, when the gap function can no longer hold vortices on the Matsubara
axis.

We also obtain $T_{c}$ as a function of $\mu$ for the more realistic model
which includes the full momentum and frequency dependence of the polarization
bubble, which accounts for the screening of the Coulomb interaction. We find
that $T_{c}$ remains finite in the limit $\mu\to 0$, as long as $f<f^{*}$.

The paper is structured as follows: In Section II, we introduce and motivate
our model. In Sec.II.1 we present the linearized Eliashberg equations, which
we use to calculate $T_{c}$ and the gap function at $T_{c}$. In Sec. II.2 we
briefly review the effect of repulsion on superconductivity at high-density.
In Sec. III we present analytical results in the low-density limit. In Section
IV we discuss our numerical analysis. We review the numerical methods we use
to solve the Eliashberg equations in Section IV.1, and present the results of
our calculations in Sec. IV.2. In particular, we show (i) how $T_{c}$ varies
with the chemical potential $\mu$, (ii) how $T_{c}$ is affected by the
strength of the repulsive component of the interaction $f$, (iii) how the gap
function depends on Matsubara frequency, (iv) how the location of the nodal
point of $\Delta(\omega_{m})$ depends on the strength of the repulsive
interaction, and (v) how the gap function depends on momenta away from
$k=k_{F}$. In Sec. IV.3 we analyze $T_{c}$ and the gap function in a model
with dynamically-screened Coulomb interaction. We present our conclusions in
Sec. V.

##  II Model

We consider an electron gas in 3 dimensions, with dispersion
$\xi(\mathbf{k})=k^{2}/2m-\mu$. Electrons interact via the Coulomb potential
and through exchange of lattice vibrations, which effectively screen the
electron charge. We follow Ref. 20 and approximate the total (direct and
phonon-mediated) interaction between electrons by

$V(\Omega,q)=\frac{4\pi e^{2}}{\varepsilon(\Omega)q^{2}-4\pi
e^{2}\Pi(\Omega,q)},$ (2)

where $\Omega$ is a bosonic Matsubara frequency, $\varepsilon(\Omega)$ is the
dielectric function, which incorporates the screening by phonons, and
$\Pi(\Omega,q)$ is the electron polarization bubble. We take the dressed
dielectric function to be

$\varepsilon(\Omega)=\varepsilon_{\infty}\frac{\Omega^{2}+\omega_{L}^{2}}{\Omega^{2}+\omega_{T}^{2}},$
(3)

where $\omega_{L}$ and $\omega_{T}$ are the frequencies of longitudinal and
transverse optical phonons, respectively, $\omega_{L}>\omega_{T}$, and
$\varepsilon_{\infty}$ is the dielectric constant in the absence of phonons.
In the zero frequency limit,
$\varepsilon(0)=\varepsilon_{\infty}\omega_{L}^{2}/\omega_{T}^{2}$; this is
known as the Lyddane-Sachs-Teller relation [28]. Since the polarization bubble
$\Pi(\Omega,q)$ is negative for all $\Omega$ and $q$, the interaction
$V(\Omega,q)$ is positive (repulsive) at all frequencies. The phonons,
however, make this interaction frequency-dependent, even if we approximate the
polarization bubble by its static, long-wavelength limit
$\Pi(\Omega,q)\approx-2N(\mu)$. In this approximation,

$V(\Omega,q)=\frac{4\pi e^{2}}{\varepsilon(\Omega)q^{2}+\kappa^{2}}$ (4)

where $\kappa=(8\pi e^{2}N(\mu))^{1/2}$ is the Thomas-Fermi wavevector.

The interaction in the form of Eq. 3 can be rigorously justified for polar
insulators, where $\Pi(\Omega,q)=0$ and

$\displaystyle V(\Omega,q)$ $\displaystyle=$ $\displaystyle\frac{4\pi
e^{2}}{\varepsilon_{\infty}q^{2}}~{}\left(1-\frac{\omega^{2}_{L}-\omega^{2}_{T}}{\Omega^{2}+\omega^{2}_{L}}\right)$
(5) $\displaystyle=\frac{4\pi
e^{2}}{{\tilde{\varepsilon}}_{\infty}q^{2}}\left(f-\frac{\omega^{2}_{L}}{\Omega^{2}+\omega^{2}_{L}}\right)$

where $f=1/(1-\omega^{2}_{T}/\omega^{2}_{L})$ and
${\tilde{\varepsilon}}_{\infty}=\varepsilon_{\infty}f$. In Ref. [22] the
authors applied a similar model with $f=1$ to a polar crystal with a finite
density of conduction electrons, appropriate for SrTiO3. Eq. 3 can also be
justified for a non-polar crystal with a monoatomic basis, where
$\varepsilon_{\infty}=1$ and $\omega_{T}=0$ as there are no transverse
phonons. In this case, at finite electron density, we have

$V^{\text{BP}}(\Omega,q)=\frac{4\pi
e^{2}}{q^{2}+\kappa^{2}}\bigg{(}1-\frac{\omega^{2}_{q}}{\Omega^{2}+\omega^{2}_{q}}\bigg{)},$
(6)

where $\omega_{q}=q\omega_{L}/\sqrt{q^{2}+\kappa^{2}}$ is the phonon
frequency. This is known as the Bardeen-Pines model. The frequency
$\omega_{q}\approx\omega_{L}$ at $q\gg\kappa$ and becomes linear in $q$ for
$q\ll\kappa$ due to electronic screening.

For most of the paper, we follow Ref. [20] and use the semi-phenomenological
form of $V(\Omega,q)$:

$V(\Omega,q)=\frac{4\pi
e^{2}}{q^{2}+\kappa^{2}}\big{(}f-\frac{\omega_{L}^{2}}{\Omega^{2}+\omega_{L}^{2}}\big{)}.$
(7)

We set $\omega_{L}$ to be a constant ($\omega_{L}=0.1$eV) and treat $f$ as a
parameter, which we vary. $V(\Omega,q)$ interpolates between Eqs. (5) and (6)
and can be thought of as an extended Bardeen-Pines model. We keep $\kappa$
finite, but will chiefly focus on the low-density limit, where typical $q$ are
much larger then $\kappa$. In this situation, $f=1$ corresponds to the
Bardeen-Pines interaction for a non-polar crystal, while for $f>1$ the
interaction closely mirrors the electron-electron interaction in a polar
crystal.

Keeping $f$ as a parameter will also allow us to connect to previous work
[20], which considered Eq. (1) in the purely-attractive $f=0$ limit. We later
extend the model by replacing $\kappa^{2}$ with $-4\pi e^{2}\Pi(\Omega,q)$ and
show that the key results, obtained with Eq. 7, survive.

###  II.1 Equations for the fermionic self-energy and the pairing vertex

The interaction $V(\Omega,q)$ gives rise to corrections to the fermionic
dispersion and the fermionic residue, while also mediating pairing between
fermions. We assume that the fermionic self-energy can be evaluated in the
one-loop approximation and the pairing vertex can be evaluated in the ladder
approximation, both using dressed Green’s functions for the intermediate
fermions. These approximations amount to neglecting vertex corrections to the
interaction. At large density ($\mu/\omega_{L}\gg 1$), such approximations can
be justified by invoking Migdal’s theorem [29]. However, for
$\mu\leq\omega_{L}$, there is no rigorous justification for neglecting vertex
corrections. The authors of Ref. 20 argued that for $f=0$, vertex corrections
are of order one and do not affect the results qualitatively. We assume that
this holds also for finite $f$.

Neglecting vertex corrections, we obtain a set of three coupled equations for
the inverse quasiparticle residue $Z_{n}(\varepsilon)$, the pairing vertex
$\phi_{n}(\varepsilon)$, and the correction to fermionic dispersion
$\chi_{n}(\varepsilon)$. Here the index $n$ refers to Matsubara frequency
$\omega_{n}=(2n+1)\pi T$ and $\varepsilon$ is the quasiparticle dispersion
$\varepsilon_{k}=k^{2}/2m$. The two variables $n$ and $\varepsilon$
parameterize the frequency and momentum dependence of the residue and the
correction to the dispersion and of the pairing vertex. We consider only
$s$-wave pairing, where the pairing vertex has no angular dependence, and
focus on $T=T_{c}$, where the pairing vertex is infinitesimally small. The
three equations are

$Z_{n}(\varepsilon)-1=-T\frac{1}{\omega_{n}}\sum_{m}\int_{0}^{\infty}d\varepsilon^{\prime}N(\varepsilon^{\prime})V_{n-m}^{\mathrm{s-wave}}(\varepsilon,\varepsilon^{\prime})\frac{\omega_{m}Z_{m}(\varepsilon^{\prime})}{[\omega_{m}Z_{m}(\varepsilon^{\prime})]^{2}+[\varepsilon^{\prime}-\mu+\chi_{m}(\varepsilon^{\prime})]^{2}}$
(8)
$\phi_{n}(\varepsilon)=-T\sum_{m}\int_{0}^{\infty}d\varepsilon^{\prime}N(\varepsilon^{\prime})V_{n-m}^{\mathrm{s-wave}}(\varepsilon,\varepsilon^{\prime})\frac{\phi_{m}(\varepsilon^{\prime})}{[\omega_{m}Z_{m}(\varepsilon^{\prime})]^{2}+[\varepsilon^{\prime}-\mu+\chi_{m}(\varepsilon^{\prime})]^{2}}$
(9)
$\chi_{n}(\varepsilon)=T\sum_{m}\int_{0}^{\infty}d\varepsilon^{\prime}N(\varepsilon^{\prime})V_{n-m}^{\mathrm{s-wave}}(\varepsilon,\varepsilon^{\prime})\frac{\chi_{m}(\varepsilon^{\prime})+\varepsilon^{\prime}-\mu}{[\omega_{m}Z_{m}(\varepsilon^{\prime})]^{2}+[\varepsilon^{\prime}-\mu+\chi_{m}(\varepsilon^{\prime})]^{2}},$
(10)

and the gap function is given by
$\Delta_{n}(\varepsilon)=\phi_{n}(\varepsilon)/Z_{n}(\varepsilon)$.

In our modified Bardeen-Pines model, we obtain

$\displaystyle V_{n-m}^{\mathrm{s-wave}}(\varepsilon_{k},\varepsilon_{q})$
$\displaystyle=\int_{-1}^{1}\frac{d\cos\theta}{2}V_{n-m}(\sqrt{k^{2}+q^{2}-2kq\cos\theta})$
(11) $\displaystyle=\frac{\pi
e^{2}}{kq}\ln\bigg(\frac{(k+q)^{2}+\kappa^{2}}{(k-q)^{2}+\kappa^{2}}\bigg{missing})u_{n-m},$
(12)

where we have defined
$u_{n}=f-\omega_{L}^{2}/(\omega_{n}^{2}+\omega_{L}^{2})$. For notational
convenience, we will drop the s-wave super-script on the interaction.

####  II.1.1 Migdal-Eliashberg Approximation

In the Migdal-Eliashberg approximation, the integrals over the dispersion
(i.e., over $\varepsilon^{\prime}$) are evaluated by linearizing the integrand
about the chemical potential $\mu$. In this approximation, the variation of
the density of states $N(\varepsilon)$ and the interaction
$V_{n-m}(\varepsilon,\varepsilon^{\prime})$ with $\varepsilon$,
$\varepsilon^{\prime}$ is ignored, and both quantities are evaluated at
$\varepsilon=\mu$. With this approximation, the integrals over energy can be
evaluated. One then finds that $\chi_{n}(\varepsilon)=0$, while $Z_{n}\equiv
Z_{n}(\mu)$ and $\phi_{n}\equiv\phi_{n}(\mu)$ only depend on Matsubara
frequency:

$Z_{n}=1-\pi
N(\mu)\frac{T}{\omega_{n}}\sum_{m}V_{n-m}\operatorname{sgn}(\omega_{m})$ (13)
$\phi_{n}=-\pi
N(\mu)T\sum_{m}V_{n-m}\frac{\phi_{m}}{\absolutevalue{\omega_{m}}Z_{m}}.$ (14)

These two equations are known as the Eliashberg equations. The approximation
$\varepsilon,\varepsilon^{\prime}\approx\mu$ used above to obtain these
equations, is valid if pairing only involves fermions near the Fermi surface.
This holds in the adiabatic limit where $\mu\gg\omega_{L}$, but not when
$\mu\leq\omega_{L}$. Since we are concerned with the low-density limit, we
instead work with the full set of equations, Eqs. (8-10).

###  II.2 The High-Density Limit

Before we proceed, we briefly review what is known in the high-density limit,
$\mu\gg\omega_{L}$. In this limit, typical $q$ are large and of order $k_{F}$.
To first approximation, the momentum dependence of the interaction can then be
approximated by $V_{0}\sim e^{2}/k^{2}_{F}$, and the full interaction can be
approximated as $V_{0}(f-\omega^{2}_{L}/(\Omega^{2}+\omega^{2}_{L}))$. At weak
coupling, when $\lambda=N(\mu)V_{0}\sim r_{s}$ is small, the fermionic self-
energy can be neglected, and the Eliashberg equation for the pairing vertex
can be solved with $Z_{n}(\varepsilon)=1$. In this case, $T_{c}$ is finite for
$f$ below a certain cutoff $f_{c}>1$ which depends on $\lambda$ [25]. In other
words, Cooper pairing continues to exist even if the interaction is purely
repulsive ($f>1$), as long as the repulsion is sufficiently weak ($f<f_{c}$).
However, for any $f>0$, the gap function changes sign as a function of
Matsubara frequency. For $f>1$, this sign change allows the system to
partially neutralize the “average” repulsion and gain from reduction of the
repulsion at small frequencies 111This is analogous to the purely momentum-
dependent case, where the gap changes sign as a function of momentum to
orthogonalize against a repulsive s-wave component of the interaction and take
advantage of an attractive component in a different angular-momentum channel
[10]..

An illustrative toy model highlighting these points was introduced by
Rietschel and Sham [19]. It mirrors the frequency dependence of the modified
Bardeen-Pines model:

$V_{mn}=\begin{cases}0,\quad\absolutevalue{\omega_{m}}>E_{c}\quad\text{or}\quad\absolutevalue{\omega_{n}}>E_{c}\\\
U(f-\Theta(\omega_{L}-\absolutevalue{\omega_{m}})\Theta(\omega_{L}-\absolutevalue{\omega_{n}})),&\text{otherwise.}\end{cases}$
(15)

In this model, electrons experience a repulsion of magnitude $Uf$ for all
frequencies below some cutoff $E_{c}$. However, if both electrons have
frequency smaller than $\omega_{L}$, there is an additional attractive term
$-U$. Using this interaction and assuming weak coupling, $UN(\mu)\ll 1$, one
can show that $T_{c}$ comes from fermions in the vicinity of the Fermi level
and is given by

$T_{c}=1.13\omega_{L}e^{-1/\lambda},$ (16)

where

$\lambda=N(\mu)U\bigg{(}1-\frac{f}{1+fN(\mu)U\ln(1.13E_{c}/\omega_{L})}\bigg{)}.$
(17)

We see that an increase in either $N(\mu)$ or $U$ enhances $T_{c}$. In
particular, by increasing $N(\mu)$ or $U$, we increase the prefactor in Eq. 17
and further reduce the repulsive contribution. Conversely, if either $N(\mu)$
or $U$ is reduced, the effective repulsive contribution is enhanced relative
to the attractive term. Superconductivity exists as long as $f<f_{c}$, where

$f_{c}=\frac{1}{1-N(\mu)U\ln(1.13E_{c}/\omega_{L})}>1$ (18)

Note that for $f\leq 1$, superconductivity exists for arbitrarily small values
of $UN(\mu)$ and $T_{c}\propto\exp(-1/N(\mu)U)$. If $f=1$,
$\lambda\propto(N(\mu)U)^{2}$, so that $T_{c}\propto\exp(-1/(N(\mu)U)^{2})$.
This behavior is also seen in the model with
$V(\Omega)\propto\left(f-\frac{\omega_{L}^{2}}{\Omega^{2}+\omega_{L}^{2}}\right)$
[25].

The gap function in the Rietschel-Sham model has a low-frequency component
$\Delta_{1}$ and a high-frequency component $\Delta_{2}$. The two are of
opposite sign, and are related by

$\Delta_{2}=-\frac{N(\mu)Uf}{1+fN(\mu)U\ln(1.13E_{c}/\omega_{L})}\ln(\frac{1.13\omega_{L}}{T_{c}})\Delta_{1}.$
(19)

As $f$ increases towards $f_{c}$, $T_{c}$ decreases. From the above formula,
this implies that the high-frequency gap becomes more and more negative, with
the ratio $\Delta_{2}/\Delta_{1}$ diverging as $T_{c}$ goes to zero. However,
the position of the sign change of the $\Delta(\omega_{m})$ is fixed at
$\omega_{0}=\omega_{L}$. This is a consequence of the fact that the boundary
between low-frequency and high-frequency regimes in the Rietschel-Sham model
is fixed at $\omega_{L}$. In the model with
$V(\Omega)\propto\left(f-\frac{\omega_{L}^{2}}{\Omega^{2}+\omega_{L}^{2}}\right)$,
the position of the zero of $\Delta(\omega_{m})$ at $\omega_{m}=\omega_{0}$ is
set by the solution of the gap equation and varies with $f$. In the high-
density limit ($\mu\gg\omega_{L}$), $\omega_{0}$ is nearly infinite at
infinitesimally small $f$ and tends to zero as $f$ approaches the critical
$f_{c}(\mu)$ from below [31]. We show below that the same holds in the low-
density limit.

##  III Analytical results in the Low Density Limit

In our analytical study we follow Refs. [20, 25, 31]. Taking $\mu\to 0$, our
goal is to find the critical $f^{*}$ where $T_{c}$ vanishes, the relation
between $T_{c}$ and $f^{*}-f$ (in the limit where $T_{c}\ll\omega_{L}$ and
$f^{*}-f\ll 1$), the relation between the position of the gap node
$\omega_{0}$ and $f^{*}-f$, and the frequency dependence of the gap near
$f^{*}$. Note that we use $f^{*}$ instead of $f_{c}$ to distinguish between
the low and high-density limits.

For nonzero $\omega_{L}$, the system is in the Fermi liquid regime, implying
that the inverse quasiparticle residue $Z_{n}(\varepsilon)$ tends to a
constant at small frequencies and approaches $1$ at large frequencies, while
the correction to the dispersion $\chi_{n}(\varepsilon)$ is non-singular. For
an order-of-magnitude analysis, we can then set $Z_{n}(\varepsilon)=1$ and
neglect $\chi_{n}(\varepsilon)$. Eq. 9 for the pairing vertex
$\phi_{n}(\varepsilon)$ is then essentially the gap equation. Introducing
$p=\sqrt{2m\varepsilon}$ and re-scaling all variables by $\omega_{L}$ as
${\bar{T}}=T/\omega_{L}$, $\bar{p}=p/\sqrt{2m\omega_{L}}$, we re-express Eq. 9
as

$\displaystyle{\phi}_{n}({\bar{p}})$ $\displaystyle=$
$\displaystyle-{\bar{T}}\sum_{m}\frac{2\sqrt{\bar{\rho}}}{\pi}\left(f-\frac{1}{1+({\bar{\omega}}_{n}-{\bar{\omega}}_{m})^{2}}\right)$
(20)
$\displaystyle\times\int_{0}^{\infty}d\bar{p}^{\prime}\frac{\frac{\bar{p}^{\prime}}{\bar{p}}\log\left(\frac{\bar{p}+\bar{p}^{\prime}}{|\bar{p}-\bar{p}^{\prime}|}\right)}{\bar{p}^{{}^{\prime}4}+{\bar{\omega}}^{2}_{m}}{\phi}_{m}(\bar{p}^{\prime}),$

where we have introduced $\bar{\rho}=\mathrm{Ry}/\omega_{L}$ and
$\mathrm{Ry}=me^{4}/2=13.6$ eV is the Rydberg energy. One can show that
$\phi_{m}(\bar{p})$ is independent of $\bar{p}$ for $\bar{p}\ll 1$ and decays
as $1/\bar{p}^{2}$ for $\bar{p}\gg 1$. Setting $\bar{p}\ll 1$ and expanding
the logarithm in ${\bar{p}}$, we find

$\displaystyle{\phi}_{n}$ $\displaystyle=$
$\displaystyle-{\bar{T}}\sum_{m}\frac{4\sqrt{\bar{\rho}}}{\pi}\left(f-\frac{1}{1+({\bar{\omega}}_{n}-{\bar{\omega}}_{m})^{2}}\right)$
(21)
$\displaystyle\times\int_{0}^{\infty}d\bar{p}^{\prime}\frac{1}{\bar{p}^{{}^{\prime}4}+{\bar{\omega}}^{2}_{m}}{\phi}_{m}(\bar{p}^{\prime}),$

where we use as shorthand $\phi_{n}\equiv\phi_{n}(\bar{p}=0)$. Since the
majority of the weight in the $\bar{p}^{\prime}$ integral comes from
$\bar{p}^{\prime}\sim|{\bar{\omega}}_{m}|^{1/2}\leq 1$, we can replace
$\phi_{m}(\bar{p}^{\prime})$ with $\phi_{m}(\bar{p}^{\prime}=0)$ on the right-
hand-side. Integrating then over $\bar{p}^{\prime}$, we obtain

$\displaystyle\phi_{n}$
$\displaystyle=-g\pi{\bar{T}}\sum_{m}\left(f-\frac{1}{1+({\bar{\omega}}_{n}-{\bar{\omega}}_{m})^{2}}\right)\frac{\phi_{m}}{|{\bar{\omega}}_{m}|^{3/2}},$
(22)

where $g=(2\bar{\rho})^{1/2}/\pi$.

To analyze the structure of $\phi_{n}=\phi(\omega_{n})$ it is convenient to
replace the sum over Matsubara frequencies by an integral and set the lower
cutoff of the integral over $\bar{\omega}^{\prime}$ at $O({\bar{T}})$. Doing
so, we obtain

$\displaystyle\phi(\bar{\omega}_{n})=-g\int_{O(\bar{T})}^{\infty}d\bar{\omega}_{m}\times$
(23)
$\displaystyle\left(f-\frac{1}{2}\left(\frac{1}{1+(\bar{\omega}_{n}-\bar{\omega}_{m})^{2}}+\frac{1}{1+(\bar{\omega}_{n}+\bar{\omega}_{m})^{2}}\right)\right)\frac{\phi(\bar{\omega}_{m})}{(\bar{\omega}_{m})^{3/2}}.$

From this equation, we have

$\displaystyle\phi(0)$
$\displaystyle=-g\int_{O(\bar{T})}^{\infty}d\bar{\omega}_{m}\left(f-\frac{1}{1+\bar{\omega}_{m}^{2}}\right)\frac{\phi(\bar{\omega}_{m})}{(\bar{\omega}_{m})^{3/2}}$
(24)

and

$\displaystyle\phi(\bar{\omega}_{n})$
$\displaystyle=\phi(0)\left(1-gQ\frac{\bar{\omega}_{n}^{2}}{1+\bar{\omega}_{n}^{2}}\right)+...,$
(25)

where
$Q=\int_{O(\bar{T})}^{\infty}d\bar{\omega}^{\prime}/(\bar{\omega}^{\prime})^{3/2}\sim
1/{\bar{T}}^{1/2}$ and the unwritten terms account for $O(g)$ corrections,
which are irrelevant for $g\leq 1$. Substituting $\phi(\bar{\omega}_{n})$ from
Eq. 25 into Eq. 24, we obtain the following self-consistent equation for
$\bar{T}_{c}$:

$\frac{1-\beta g}{1-\alpha g}-f=f^{*}-f=\frac{1}{gQ(1-\alpha g)}.$ (26)

Here,
$\alpha=\int_{0}^{\infty}d\bar{\omega}^{\prime}(\bar{\omega}^{\prime})^{1/2}/(1+(\bar{\omega}^{\prime})^{2})=\pi/\sqrt{2}\approx
2.22$ and
$\beta=\int_{0}^{\infty}d\bar{\omega}^{\prime}(\bar{\omega}^{\prime})^{1/2}/(1+(\bar{\omega}^{\prime})^{2})^{2}=\pi/(4\sqrt{2})\approx
0.56$. Note that $f^{*}>1$, since $\alpha>\beta$. Using
$Q\sim{\bar{T}}^{1/2}$, we find the scaling relation

${\bar{T}}_{c}\sim(f^{*}-f)^{2}.$ (27)

Next, from Eq. 25 we see that at large $Q$, i.e small ${\bar{T}}_{c}$, the
frequency $\bar{\omega}_{0}$ at which $\phi(\bar{\omega}_{n})$ changes sign,
is

$\bar{\omega}_{0}\approx\frac{1}{(gQ)^{1/2}}=(f^{*}-f)^{1/2}.$ (28)

For smaller $f$, this expression extends to
$\bar{\omega}_{0}\sim((f^{*}-f)/f)^{1/2}$. Since $\bar{\omega}_{0}$ tends to
zero as $f$ approaches $f^{*}$, the gap $\phi(\bar{\omega}_{n})$ changes sign
at progressively smaller $\bar{\omega}_{0}$. The vanishing of $\omega_{0}$ as
$f$ approaches $f^{*}$ is consistent with the behavior near $f_{c}$ in the
high-density limit [31]. However, the power-law behavior of $T_{c}$ and the
relation $\omega_{0}\sim\left(T_{c}\right)^{1/4}$ are specific to the case of
low-density.

Because
$\phi(\bar{\omega}_{n})=\phi(0)(1-gQ\bar{\omega}_{n}^{2}/(1+\bar{\omega}_{n}^{2}))$
and $Q\propto 1/(f^{*}-f)$, the ratio $\phi(\bar{\omega}_{n})/\phi(0)$ becomes
more negative with increasing $\bar{\omega}_{n}$, going as
$\phi(\bar{\omega}_{n}\gg 1)/\phi(0)\propto 1/\bar{T}_{c}^{1/2}$. This is
shown in the top panel of Fig. 1.

Figure 1:  The qualitative behavior of the gap function $\phi_{m}(\bar{p})$ at
$\mu=0$ as a function of $\bar{\omega}_{m}$ in the top panel and as a function
of $\bar{p}$ in bottom panel). In addition, we highlight how the gap node
$\omega_{0}$, the high-frequency gap ratio
$\phi_{\infty}(\bar{p}=\infty)/\phi_{0}(\bar{p}=0)$, and the peak in
$\phi_{0}(\bar{p})/\phi_{0}(\bar{p}=0)$ scale near $f^{*}-f$.

Thus far, we have considered only the $\bar{p}\ll 1$ limit for
$\phi_{n}(\bar{p})$, since this is sufficient to obtain $\bar{T}_{c}$. We now
examine the behavior of $\phi_{n}(\bar{p})$ at $\bar{p}\geq 1$. At large
${\bar{p}}$, one can easily verify that $\phi_{n}(\bar{p})\propto
1/{\bar{p}}^{2}$. Accordingly, we introduce the function $B_{n}$ via
$\phi_{n}(\bar{p}\gg 1)=B_{n}/\bar{p}^{2}$. Substituting this into Eq. 20,
taking $\bar{p}\gg 1$, and setting $n=0$ for definiteness, we obtain

$\displaystyle B_{0}$
$\displaystyle\approx-\bar{T}\sum_{m}\frac{4\sqrt{\bar{\rho}}}{\pi}\left(f-\frac{1}{1+{\bar{\omega}}^{2}_{m}}\right)$
$\displaystyle\times\bigg{(}\phi_{m}\int_{0}^{1}d\bar{p}^{\prime}\frac{\bar{p}^{{}^{\prime}2}}{\bar{p}^{{}^{\prime}4}+\bar{\omega}^{2}_{m}}+B_{m}\int_{1}^{\bar{p}}d\bar{p}^{\prime}\frac{\bar{p}^{{}^{\prime}2}}{\bar{p}^{{}^{\prime}4}+\bar{\omega}^{2}_{m}}$
$\displaystyle+B_{m}\int_{\bar{p}}^{\infty}d\bar{p}^{\prime}\left(\frac{{\bar{p}}}{{\bar{p}}^{\prime}}\right)^{2}\frac{1}{\bar{p}^{{}^{\prime}4}+\bar{\omega}^{2}_{m}}\bigg{)}.$
(29)

To obtain the first term of the second line, we use our earlier result that
for $\bar{p}\ll 1$, $\phi_{m}(\bar{p})$ is independent of momentum and replace
$\phi_{m}(\bar{p})$ with $\phi_{m}\equiv\phi_{m}(\bar{p}=0)$ for $\bar{p}<1$.
Similarly, we replace $\phi_{m}(\bar{p})$ with its asymptotic limit
$B_{n}/\bar{p}^{2}$ for $\bar{p}>1$ to obtain the latter two terms of the
above equation. To simplify this equation, we note that the last term on the
right-hand side is smaller than the first two terms by a factor of order
$O(1/\bar{p}^{3})$ and is therefore irrelevant. Discarding this term and
taking $\bar{p}\to\infty$ in the upper limit of the second integral, we then
find

$\displaystyle B_{0}$
$\displaystyle\approx-\bar{T}\sum_{m}\frac{4\sqrt{\bar{\rho}}}{\pi}\left(f-\frac{1}{1+{\bar{\omega}}^{2}_{m}}\right)$
$\displaystyle\times\bigg{(}\phi_{m}\int_{0}^{1}d\bar{p}^{\prime}\frac{\bar{p}^{{}^{\prime}2}}{\bar{p}^{{}^{\prime}4}+\bar{\omega}^{2}_{m}}+B_{m}\int_{1}^{\infty}d\bar{p}^{\prime}\frac{\bar{p}^{{}^{\prime}2}}{\bar{p}^{{}^{\prime}4}+\bar{\omega}^{2}_{m}}\bigg{)}$
(30)

One can verify that all momentum and frequency integrals appearing in this
equation are $O(1)$, i.e. nonsingular in the $T=0$ limit. Assuming that the
primary contribution to $B_{0}$ on the right-hand side of the above equation
comes from $\phi_{m}$, not $B_{m}$, we find that $B_{0}$ is determined by
$\phi_{m}$ at frequencies ${\bar{\omega}}_{m}=O(1)$, where $\phi_{m}\sim-
gQ\phi_{0}$.

We then have $B_{0}\sim g^{2}Q\phi_{0}$, which, recalling that $Q\sim
1/T_{c}^{1/2}$, implies that $B_{0}/\phi_{0}$ diverges as $T_{c}\to 0$ 222One
can verify that our assumption (that $\phi_{m}$ matters much more than
$B_{m}$) is valid by inserting our result, $B_{m}\sim g\phi_{m}$ back into Eq.
30.. In other words, the ratio $\phi_{0}(\bar{p})/\phi_{0}(\bar{p}=0)$ must
have a large peak as a function of $\bar{p}$, with its magnitude growing as
$1/(f^{*}-f)$. We show this behavior in
$\phi_{0}(\bar{p})/\phi_{0}(\bar{p}=0)$ in the bottom panel of Fig. 1.

Note also that $T_{c}$ is finite even if one does not impose an ultraviolet
cutoff on the frequency integration. This is due to the $1/q^{2}$ momentum
dependence of the interaction, which leads to the momentum integration in the
particle-particle bubble going as $1/\omega^{3/2}$. In this case, the integral
$\int d\omega/\omega^{3/2}$ converges in the ultraviolet.

##  IV Numerical analysis

###  IV.1 Methods

To solve the linearized Eliashberg equations for $T_{c}$, we note that the
equation for the pairing vertex $\phi_{n}(\varepsilon)$ is essentially an
eigenvalue problem. To solve this eigenvalue problem, we create a linear
operator mapping $\phi_{m}(\varepsilon^{\prime})$ to $\phi_{n}(\varepsilon)$
in Eq. 9 for different values of temperature. To find $T_{c}$, one must then
find at what temperature this linear operator’s largest positive eigenvalue is
equal to 1.

To construct this operator at a given temperature $T$, we first solve for
$Z_{n}(\varepsilon)$ and $\chi_{n}(\varepsilon)$. This is done self-
consistently by iterating Eq. 8) and Eq. 10 starting from
$Z_{n}^{\text{initial}}(\varepsilon)=1$ and
$\chi_{n}^{\text{initial}}(\varepsilon)=0$ until convergence is reached. The
energy integrals are obtained using upper cutoffs from $\Lambda=100\omega_{L}$
to $\Lambda=200\omega_{L}$, and a grid of hundreds of sampling points. We
split the energy range into 3 regions, (i) $\varepsilon<\mu-\delta$, (ii)
$\mu-\delta<\varepsilon<\mu+\delta$, and (iii)
$\mu+\delta<\varepsilon<\Lambda$, where we take $\delta=\mu/100$. In region
(ii) near the chemical potential, we use a high density of quadrature points
to account for the peak in the integrand and apply the trapezoidal rule. In
regions (i) and (iii) where the variation in the integrand is smoother, we use
Gauss-Legendre quadrature to calculate the integrals. In the low-density limit
$\mu=0$, we use a composite Gauss-Legendre grid with around 1000 points, to
ensure that we accurately obtain contributions from all values of
$\varepsilon$. We find that trends in $T_{c}$ are well-converged with respect
to variations in $\delta$ and the number of quadrature points.

To calculate the Matsubara sums, we note that all Matsubara sums appearing in
the Eliashberg equations are convolutions. These convolutions can be
efficiently calculated by first transforming to imaginary time, where the
convolution becomes point-wise multiplication. After point-wise
multiplication, one can then transform back to Matsubara frequency. When using
the Fast Fourier transform, this method scales significantly better,
$O(N\log(N))$, than naively calculating the sums directly in Matsubara space,
$O(N^{2})$.

Though this method works well for larger values of $T_{c}$, we run into memory
issues when trying to extend this method to temperatures smaller than
$T_{c}\sim 10^{-4}\omega_{L}$. We obtained $T_{c}$ in this temperature range
by extrapolating from higher temperature data using the implicit
renormalization method [25], which also allows us to infer the gap function at
$T_{c}$ from less memory-intensive high-temperature calculations. To apply
this method, we divide Eqs. (8-10) into low-energy and high-energy components.
We then use the high-energy components to obtain an effective gap equation for
the low-energy component of the gap. A gap component is considered high-energy
if its respective momentum and frequency satisfy
$\sqrt{(k^{2}/2m-\mu)^{2}+\omega_{m}^{2}}>\Omega_{c}$, where $\Omega_{c}$ is
some cutoff frequency. For consistency, we take $\Omega_{c}=\omega_{L}$ in all
calculations. If the largest eigenvalue obtained from this effective gap
equation scales linearly with $\ln(\omega_{L}/T)$, then we can extrapolate
$T_{c}$ by extrapolating the eigenvalues of less-computationally-intensive,
high-temperature calculations.

We find that this method captures the overall trend of $T_{c}$ relatively
well. In particular, the transition temperatures calculated from both the
traditional eigenvalue method and the implicit renormalization method agree
well for large temperatures, where the traditional eigenvalue method is
practicable, and for sufficiently large $\mu$. In the small-$\mu$ limit, we
find that the trend in $T_{c}$ found via the implicit renormalization method
is very similar to that found using the eigenvalue method. However, due to the
shortcomings of this method in the low-density limit, we calculate $T_{c}$ and
other quantities using the standard eigenvalue method at low density when
possible.

Lastly, we find in our calculations that the function $\chi_{n}(\varepsilon)$
does not vary significantly with momentum or frequency, and can be largely
absorbed into the definition of the chemical potential. Therefore, we expect
all results to be relatively insensitive to whether $\chi_{n}(\varepsilon)$ is
included or set to zero. As such, we take $\chi_{n}(\varepsilon)=0$ in the
following calculations and solve only the two coupled equations for
$Z_{n}(\varepsilon)$ and $\phi_{n}(\varepsilon)$.

###  IV.2 Results

####  IV.2.1 $T_{c}$ vs. $\mu$ and $f$

In Fig. 2, we show how $T_{c}$ varies with the chemical potential $\mu$ for
different values of the repulsive term $f$. For consistency, we use the
implicit-renormalization method to extract $T_{c}$ for all values of $f$ and
$\mu$ presented here. We see that $T_{c}$ is enhanced as one approaches
$\mu=0$, regardless of the strength of the repulsive term. One can understand
this enhancement in the same way as was done in the purely attractive case
[20]. Namely, as $\mu$ decreases, $N(\mu)V_{n-m}(\mu,\mu)$ is enhanced at low
density due to the reduction in screening. Additionally, pairing is no longer
restricted to occur in a narrow window around $\mu$. We note that a similar
trend in $T_{c}$ as a function of density has been found by Takada [21], who
solved the full Eliashberg equations in a multi-valley electron gas to study
plasmon-induced superconductivity.

Figure 2: $T_{c}$ as a function of chemical potential $\mu$ for different
strengths $f$ of the repulsive interaction. Figure 3:  The $\mu=0$ scaling of
$T_{c}$ as $f$ approaches $f^{*}$. Overlaid is the line showing that
$T_{c}\propto(f^{*}-f)^{2}$ near $f=f^{*}$.

Another trend we find in Fig. 2 is that $T_{c}$ drops more precipitously with
increasing $\mu$, passing through a local minimum as a function of
$\mu/\omega_{L}$. The presence of this minimum can understood as follows: as
$\mu$ increases from $0$, the effective interaction decreases and the relevant
values of $\varepsilon$ gradually cluster closer to $\varepsilon=\mu$. Both
factors lead to a reduction of $T_{c}$ as $\mu$ is increased. However, as
$\mu$ increases, the range in $\varepsilon$ about which
$1/(Z^{2}\omega_{m}^{2}+(\varepsilon-\mu)^{2})$ in the integrand is large also
increases, leading to an increase in $T_{c}$. The first two factors, which
desire a decrease in $T_{c}$, dominate at small $\mu$, while the latter
factor, which desires an increase in $T_{c}$, dominates at larger $\mu$.
Together, these competing factors lead to the local minimum seen in Fig. 2.

We also see from Fig. 2 that the value of $T_{c}$ at the minimum rapidly
decreases as $f$ is increased. This is particularly prominent at larger values
of $f$, as one can see from the inset of Fig. 2, where we take $f=1.4$. Here,
the minimum in $T_{c}$ is significantly more pronounced. As such, we find that
our model effectively exhibits re-entrant superconductivity for larger values
of $f$. That is, if one starts at large density and lowers the chemical
potential, $T_{c}$ drops to zero, and then grows from zero to a constant as
the chemical potential is further decreased and approaches zero. However, as
is clear from the inset, the values of $T_{c}$ are so small that it is likely
experimentally infeasible to observe this re-entrant superconductivity.

Finally, our results show that at $\mu=0$ there exists a critical $f^{*}$,
above which superconductivity does not develop. We present our results for
$T_{c}$ near $f^{*}$ in Fig. 3. We clearly see a power-law dependence of
$T_{c}$ on $f^{*}-f$ which agrees with our analytical result,
$T_{c}\propto(f^{*}-f)^{2}$. We find $f^{*}\approx 4.9$ by performing a linear
extrapolation of the calculated $(T_{c}/\omega_{L})^{1/2}$ to 0. We emphasize
that since $f^{*}>1$, there is a range of $f$ where superconductivity survives
in the $\mu=0$ limit even when the pairing interaction is repulsive at all
frequencies.

For larger values of $\bar{\mu}$ ($\bar{\mu}\gg 1$ or equivalently
$\mu\gg\omega_{L}$), we return to the adiabatic regime, where only momenta
near the Fermi level are relevant. In this case, we can track the behavior in
$T_{c}$ by following the behavior in the coupling
$\lambda(\bar{\mu})=N(\bar{\mu})V(\bar{\mu},\bar{\mu})$, which for our model
is

$\lambda(\bar{\mu})=\frac{1}{2\pi}\sqrt{\frac{\bar{\rho}}{\bar{\mu}}}\log(1+\pi\sqrt{\frac{\bar{\mu}}{\bar{\rho}}}).$
(31)

From this, we see that $\lambda(\bar{\mu})$ decays
$\log(\bar{\mu})/\sqrt{\bar{\mu}}$ at $\bar{\mu}\gg\bar{\rho}$. Accordingly,
for all $f>1$, we expect $T_{c}$ to go to zero as $\bar{\mu}$ is increased
past some threshold $\bar{\mu}^{*}$, where the coupling
$\lambda(\bar{\mu}^{*})$ becomes too small to stabilize superconductivity.
However, we expect this destruction of superconductivity to occur at extremely
large $\bar{\mu}$ ($\bar{\mu}\gg\bar{\rho}$, where we have set
$\bar{\rho}=\mathrm{Ry}/\omega_{L}=136$) while the main focus of our work is
on the low-density limit.

####  IV.2.2 Behavior of $\Delta_{n}(\varepsilon)$ with Matsubara Frequency

We now turn to the behavior of the gap as a function of Matsubara frequency.
The results are shown in Fig. 4, where we have set $\mu=10^{-5}\omega_{L}$. We
find that for any $0<f<f^{*}$, $\Delta(\omega_{m})$ undergoes a sign change at
some nonzero $\omega_{0}$. The ratio
$\Delta(\omega_{m}\gg\omega_{L})/\Delta(\omega_{m}=0)$ becomes more negative
as $f$ increases. This fully agrees with the analytical result, presented in
Fig. 1.

Figure 4:  $\Delta(\omega_{m},\varepsilon=0)$ for different values of $f$. We
have set $\mu=10^{-5}\omega_{L}$.

We now examine how $\omega_{0}$ behaves as a function of $f$. In Fig. 5 we
show $\omega_{0}$ for small values of $f$. We find that the position of the
node $\omega_{0}$ scales as $1/f^{1/2}$ at small $f$. This scaling is the same
as at large density [31] and can be easily understood as the frequency
dependence of the gap at large $\omega_{m}$ and small $f$ follows

$\Delta(\omega_{m})\propto\frac{\omega_{L}^{2}}{\omega_{m}^{2}}-f.$ (32)

Figure 5:  The position of the gap node as a function of $1/\sqrt{f}$, for
$\mu=10^{-5}\omega_{L}$.  Figure 6:  The $\mu=0$ scaling of $\omega_{0}$ as
$f$ approaches $f^{*}$. Overlaid is the line showing that
$\omega_{0}\propto(f^{*}-f)^{1/2}$ near $f=f^{*}$.

In Fig. 6 we set $\mu=0$ and show how $\omega_{0}$ varies as $f$ approaches
$f^{*}$. We find that $\omega_{0}$ decreases with $f$ nonlinearly, with the
slope of $\omega_{0}^{2}(f)$ decreasing with increasing $f$. The solid line in
the plot is the fit to $\omega_{0}\propto(f^{*}-f)^{1/2}$ that we obtained
analytically in Eq. 28. The fit is somewhat ambiguous as one needs more points
closer to $f^{*}$. However, the agreement with our analytics is quite
reasonable.

####  IV.2.3 Dependence of $\Delta_{n}(\varepsilon)$ on $\varepsilon$

We now turn to the behavior of the gap as a function of
$\varepsilon=k^{2}/2m$. We show the result in Fig. 7 for
$\mu=10^{-5}\omega_{L}$ for various values of $f$. Overlaid are dashed lines
which delineate the values of $T_{c}/\omega_{L}$ for each value of $f$. From
this plot, we see that there are essentially three regions of interest, (i)
$\varepsilon\ll T_{c}$, (ii) $T_{c}\ll\varepsilon\ll\omega_{L}$, and (iii)
$\varepsilon\gg\omega_{L}$. In region (i), the gap is essentially constant. In
the intermediate region (ii) where $T_{c}\ll\varepsilon\ll\omega_{L}$, we find
a smooth increase in the gap as a function of $\varepsilon$, which gets more
pronounced with increasing $f$. Lastly, in region (iii) where
$\varepsilon\gg\omega_{L}$, the gap decays as $B/\varepsilon$, where the
constant $B$ grows with increasing $f$. We note that the asymptotic behavior
we see here agrees with our analytics, where we argued that the gap should be
constant for small momenta and decay as $B/p^{2}$ (or equivalently
$B/\varepsilon$) for large momenta. Additionally, the behavior in the
intermediate region also agrees with our analytical results, where we argued
that there should be a peak in the momentum dependence of the gap, whose
magnitude scales as $1/(f^{*}-f)$ at $\mu=0$. However, we find numerically
that the magnitude of this peak grows rather slowly as $f$ is increased. This
may be due to a game of numbers.

Figure 7:  The gap function $\Delta(\omega_{m}=\pi T,\varepsilon)$ as a
function of $\varepsilon$ for different values of $f$. The dashed lines show
the values of $T_{c}/\omega_{L}$ for different $f$. We have set
$\mu=10^{-5}\omega_{L}$.

We also note that $\Delta_{n}(\varepsilon)$ is smooth near the chemical
potential. We discuss this in more details in the next section.

###  IV.3 The Effect of the Momentum and Frequency Dependence of
$\Pi(\Omega,q)$

To investigate the robustness of our results, we go beyond the extended
Bardeen-Pines model and recalculate $T_{c}$ as a function of $\mu$, by (i)
replacing $\kappa^{2}(\mu)=(-4\pi e^{2}\Pi(0,0))^{1/2}$ with $(-4\pi
e^{2}\Pi(0,q))^{1/2}$ in Eq. 7, and (ii) working with the full interaction Eq.
2. In Case (i), we take $f=1$, while in Case (ii), we set
$\varepsilon_{\infty}=1$ and $\omega_{T}=0$ in the dielectric function
$\varepsilon(\Omega)$; the latter is analogous to setting $f=1$. The results
for these calculations are presented in Fig. 8. In both cases, $T_{c}$
saturates to a nonzero value with decreasing $\mu$, in line with the results
for the extended Bardeen-Pines model. This behavior is to be expected, since
inclusion of the momentum and frequency dependence of $\Pi(\Omega,q)$ only
weakens the screening.

Figure 8:  $T_{c}$ vs $\mu$ for more general interactions. In the static
bubble case, we include the momentum-dependence of the polarization bubble in
the calculation, and set $f=1$. In the full bubble case, we use the
interaction of Eq. 2, setting $\omega_{T}=0$ in the dielectric function; this
is analogous to setting $f=1$.

In Fig.9a we show the results for $\Delta_{n}(\varepsilon)$, obtained with the
full $\Pi(\Omega,q)$, as a function of energy $\varepsilon$ (not to be
confused with the dielectric function $\varepsilon(\Omega)$) for various
Matsubara frequencies, with $\mu=10^{-5}\omega_{L}$. We see that
$\Delta_{n}(\varepsilon)$ is smooth at $\varepsilon\sim\mu$. This is
consistent with the result that we obtained in the previous section for the
extended Bardeen-Pines model.

Figure 9: The gap $\Delta_{n}(\varepsilon)$ as a function of energy
$\varepsilon$ for various values of $\omega_{m}$ in (a) the ultra-low density
limit, and (b) the high-density limit, as obtained from the interaction with
full dynamical screening. We have set $\omega_{T}=0$ in these calculations,
which is analogous to setting $f=1$. The chemical potentials
($\mu=10^{-5}\omega_{L}$ in Panel (a), $\mu=3.59\omega_{L}$ in Panel (b)), are
denoted with dotted lines.

In Fig. 9b we present results of the same calculation, but at much larger
density $\mu=3.59\omega_{L}$. We see that at small Matsubara frequencies
$\Delta_{n}(\varepsilon)$ is again smooth near the Fermi level, but at larger
$\omega_{n}$ develops a strong dip at $\varepsilon=\mu$. Such a dip has been
originally observed by Takada [21]. Richardson and Ashcroft [15] argued that
it arises from the long-range nature of the Coulomb interaction and holds when
$T_{c}/E_{F}\ll 1$ and one can linearize the dispersion around the Fermi
level.

To understand this dip we set $Z_{n}(\varepsilon)=1$ and analyze the
Eliashberg gap equation

$\phi_{n}(\varepsilon)=-T_{c}\sum_{m}\int_{0}^{\infty}d\varepsilon^{\prime}N(\varepsilon^{\prime})V_{n-m}^{\mathrm{s-wave}}(\varepsilon,\varepsilon^{\prime})\frac{\phi_{m}(\varepsilon^{\prime})}{\omega_{m}^{2}+(\varepsilon^{\prime}-\mu)^{2}}$
(33)

using

$V(\Omega,q)=\frac{4\pi e^{2}}{\varepsilon(\Omega)q^{2}-4\pi
e^{2}\Pi(\Omega,q)}.$ (34)

We assume and then verify that in the limit of large Matsubara frequency,
$\omega_{n}\to\infty$, relevant $\omega_{m}$ in the right-hand side of Eq. 33
are finite. The relevant bosonic $\Omega=\omega_{n}-\omega_{m}$ then approach
$\infty$. Since the dynamical $\Pi(\Omega,q)$ vanishes for large $\Omega$, we
have

$V_{n-m}^{\mathrm{s-wave}}(\varepsilon,\varepsilon^{\prime})=V_{\infty}^{\mathrm{s-wave}}(\varepsilon_{k},\varepsilon_{q})=\frac{2\pi
e^{2}}{kq}\ln\frac{k+q}{|k-q|}.$ (35)

becomes purely static. Solving for the gap we then find that
$|\phi_{\infty}(\varepsilon)|$ is logarithmically enhanced at
$\varepsilon=\mu$:

$\phi_{\infty}(\varepsilon=\mu)\sim\phi_{\infty}(\varepsilon\geq\mu)\ln^{2}\big{(}\frac{4\mu}{T_{c}}\big{)}$
(36)

The logarithmic singularity comes from the range near the Fermi surface, where
$k\approx q$, which in turn arises from the long-range, unscreened behavior of
the interaction at large $\omega_{n}-\omega_{m}$. We emphasize that this
singularity in $\phi_{\infty}(\varepsilon)$ exists only in the high-density
limit, and to obtain it one needs to include the full frequency dependence of
the polarization bubble $\Pi(\Omega,q)$.

We also note that in this section we calculated $\Pi(\Omega,q)$ at $T=0$. This
is valid when $T_{c}\ll\mu(T_{c})$ but is questionable as $\mu\to 0$. However,
in light of the results presented here, we expect that including the
temperature dependence of $\Pi(\Omega,q)$ should not qualitatively change our
conclusions, since its inclusion would only serve to more quickly weaken the
screening of the interaction.

##  V Conclusions

In this work, we have studied the effect of a repulsive Coulomb interaction,
on electron-phonon superconductivity in the low-density limit, the case of
pairing interaction $V(\Omega,q)=4\pi
e^{2}/(q^{2}+\kappa^{2})\times(f-\omega^{2}_{L}/(\Omega^{2}+\omega^{2}_{L}))$.
Our results show that as for the $f=0$ case of pure electron-phonon
attraction, studied in Ref. [20], $T_{c}$ is enhanced as $\mu$ decreases,
approaching a constant in the $\mu=0$ limit. We find that the gap function
changes sign at some Matsubara frequency $\omega_{0}$, reducing the effect of
the repulsion and allowing $T_{c}$ to remain nonzero over some range of
$1<f<f^{*}$, when the interaction is repulsive at all frequencies. As $f$
approaches $f^{*}$, we find that both $T_{c}$ and $\omega_{0}$ approach zero
as powers of $f^{*}-f$. This result, which we obtained both analytically and
numerically, is in contrast to the behavior in the high-density limit, where
$T_{c}$ vanishes exponentially in $f^{*}-f$.

Our results suggest that experimentally tuning the chemical potential should
lead to substantial, observable variation in $\Delta(\omega_{n})$, which can
be observed in, e.g., ARPES experiments.

Lastly, we show that the behavior we find in $T_{c}$, namely that it stays
nonzero when we take $\mu=0$, continues to hold when we include dynamical
screening of the interaction. Also, although in this work we focused on a 3D
Galilean-invariant system, the behavior we find here should be relatively
general and most likely continues to hold in two dimensions and for lattice
systems.

In this work, we have not considered the possibility of other phases. Indeed,
Wigner crystallization is also favored at low density. We leave study of the
competition between superconductivity and other phases to future work, noting
only that a superconductor to Wigner-crystal phase transition has been
previously proposed in the three-dimensional electron gas, where the electron-
electron interaction is plasmon-mediated [33], and in twisted bilayer graphene
[34].

Another item for future study is the role of phase fluctuations. The
transition temperature we calculate from the linearized Eliashberg equations
is not the true superconducting transition temperature, but the onset
temperature for pair formation. The superconducting transition temperature
$T_{c}$ is defined as the onset of phase rigidity [35] and in general should
be smaller than the onset temperature for the pairing. For obvious reasons we
expect the effect of phase fluctuations to become progressively more relevant
for quasi-2D systems.

## Acknowledgements

We thank M. Gastiasoro, R. Fernandes, D. Pimenov, J. Ruhman, Y. Wu and S.
Zhang for useful discussions. The work was supported by the Office of Basic
Energy Sciences US Department 791 of Energy under Award No. DE-SC0014402.

## References

  * Schooley  _et al._ [1964] J. F. Schooley, W. R. Hosler, and M. L. Cohen, Superconductivity in semiconducting srtio3, Phys. Rev. Lett. 12, 474 (1964).
  * Schooley  _et al._ [1965] J. F. Schooley, W. R. Hosler, E. Ambler, J. H. Becker, M. L. Cohen, and C. S. Koonce, Dependence of the superconducting transition temperature on carrier concentration in semiconducting srtio3, Phys. Rev. Lett. 14, 305 (1965).
  * Lin  _et al._ [2014] X. Lin, G. Bridoux, A. Gourgout, G. Seyfarth, S. Krämer, M. Nardone, B. Fauqué, and K. Behnia, Critical doping for the onset of a two-band superconducting ground state in ${\mathrm{srtio}}_{3-\delta}$, Phys. Rev. Lett. 112, 207002 (2014).
  * Prakash  _et al._ [2017] O. Prakash, A. Kumar, A. Thamizhavel, and S. Ramakrishnan, Evidence for bulk superconductivity in pure bismuth single crystals at ambient pressure, Science 355, 52 (2017).
  * Ruhman and Lee [2017] J. Ruhman and P. A. Lee, Pairing from dynamically screened coulomb repulsion in bismuth, Phys. Rev. B 96, 235107 (2017).
  * Chernik and Lykov [1981] I. A. Chernik and S. N. Lykov, Sov. Phys. Solid State 23, 817 (1981).
  * Kong  _et al._ [2013] P. Kong, J. Zhang, S. Zhang, J. Zhu, Q. Liu, R. Yu, Z. Fang, C. Jin, W. Yang, X. Yu, _et al._ , Superconductivity of the topological insulator bi2se3 at high pressure, Journal of Physics: Condensed Matter 25, 362204 (2013).
  * Nakajima  _et al._ [2015] Y. Nakajima, R. Hu, K. Kirshenbaum, A. Hughes, P. Syers, X. Wang, K. Wang, R. Wang, S. R. Saha, D. Pratt, J. W. Lynn, and J. Paglione, Topological rpdbi half-heusler semimetals: A new family of noncentrosymmetric magnetic superconductors, Sci. Adv. 1, e1500242 (2015).
  * Gastiasoro  _et al._ [2020] M. N. Gastiasoro, J. Ruhman, and R. M. Fernandes, Superconductivity in dilute srtio3: a review, Annals of Physics 417, 168107 (2020).
  * Maiti and Chubukov [2013] S. Maiti and A. V. Chubukov, Superconductivity from repulsive interaction, in  _AIP Conference Proceedings_ , Vol. 1550 (American Institute of Physics, 2013) pp. 3–73.
  * Marsiglio [2020] F. Marsiglio, Eliashberg theory: A short review, Annals of Physics 417, 168102 (2020), eliashberg theory at 60: Strong-coupling superconductivity and beyond.
  * Chubukov  _et al._ [2020] A. V. Chubukov, A. Abanov, I. Esterlis, and S. A. Kivelson, Eliashberg theory of phonon-mediated superconductivity — when it is valid and how it breaks down, Annals of Physics 417, 168190 (2020), eliashberg theory at 60: Strong-coupling superconductivity and beyond.
  * Chowdhury and Berg [2020] D. Chowdhury and E. Berg, The unreasonable effectiveness of eliashberg theory for pairing of non-fermi liquids, Annals of Physics 417, 168125 (2020), eliashberg theory at 60: Strong-coupling superconductivity and beyond.
  * Grabowski and Sham [1984] M. Grabowski and L. J. Sham, Superconductivity from nonphonon interactions, Phys. Rev. B 29, 6132 (1984).
  * Richardson and Ashcroft [1997] C. Richardson and N. Ashcroft, Effective electron-electron interactions and the theory of superconductivity, Physical Review B 55, 15130 (1997).
  * McMillan [1968] W. L. McMillan, Transition temperature of strong-coupled superconductors, Phys. Rev. 167, 331 (1968).
  * Tolmachev and Tyablikov [1958] V. Tolmachev and S. Tyablikov, Soviet Physics JETP 7 (1958).
  * Morel and Anderson [1962] P. Morel and P. W. Anderson, Calculation of the superconducting state parameters with retarded electron-phonon interaction, Phys. Rev. 125, 1263 (1962).
  * Rietschel and Sham [1983] H. Rietschel and L. Sham, Role of electron coulomb interaction in superconductivity, Physical Review B 28, 5100 (1983).
  * Gastiasoro  _et al._ [2019] M. N. Gastiasoro, A. V. Chubukov, and R. M. Fernandes, Phonon-mediated superconductivity in low carrier-density systems, Physical Review B 99, 094524 (2019).
  * Takada [1992] Y. Takada, Plasmon mechanism of superconductivity in the multivalley electron gas, Journal of the Physical Society of Japan 61, 238 (1992).
  * Ruhman and Lee [2016] J. Ruhman and P. A. Lee, Superconductivity at very low density: The case of strontium titanate, Phys. Rev. B 94, 224515 (2016).
  * Tchoumakov  _et al._ [2020] S. Tchoumakov, L. J. Godbout, and W. Witczak-Krempa, Superconductivity from coulomb repulsion in three-dimensional quadratic band touching luttinger semimetals, Phys. Rev. Research 2, 013230 (2020).
  * Bardeen and Pines [1955] J. Bardeen and D. Pines, Electron-phonon interaction in metals, Phys. Rev. 99, 1140 (1955).
  * Chubukov  _et al._ [2019] A. Chubukov, N. V. Prokof’ev, and B. V. Svistunov, Implicit renormalization approach to the problem of cooper instability, Phys. Rev. B 100, 064513 (2019).
  * Wu  _et al._ [2021] Y.-M. Wu, S.-S. Zhang, A. Abanov, and A. V. Chubukov, Interplay between superconductivity and non-fermi liquid at a quantum critical point in a metal. iv. the $\gamma$ model and its phase diagram at $1<\gamma<2$, Physical Review B 103, 024522 (2021).
  * Christensen and Chubukov [2021] M. H. Christensen and A. V. Chubukov, Dynamical vortices in electron-phonon superconductors, Phys. Rev. B 104, L140501 (2021).
  * Ashcroft  _et al._ [1976] N. W. Ashcroft, N. D. Mermin, _et al._ ,  _Solid state physics_ , Vol. 2005 (holt, rinehart and winston, new york London, 1976).
  * Migdal [1958] A. Migdal, Interaction between electrons and lattice vibrations in a normal metal, Sov. Phys. JETP 7, 996 (1958).
  * Note [1] This is analogous to the purely momentum-dependent case, where the gap changes sign as a function of momentum to orthogonalize against a repulsive s-wave component of the interaction and take advantage of an attractive component in a different angular-momentum channel [10].
  * [31] D. Pimenov and A. V. Chubukov, Quantum phase transition in a clean superconductor with repulsive dynamical interaction, .
  * Note [2] One can verify that our assumption (that $\phi_{m}$ matters much more than $B_{m}$) is valid by inserting our result, $B_{m}\sim g\phi_{m}$ back into Eq. 30.
  * Takada [1993] Y. Takada, s- and p-wave pairings in the dilute electron gas: Superconductivity mediated by the coulomb hole in the vicinity of the wigner-crystal phase, Phys. Rev. B 47, 5202 (1993).
  * Padhi  _et al._ [2018] B. Padhi, C. Setty, and P. W. Phillips, Doped twisted bilayer graphene near magic angles: proximity to wigner crystallization, not mott insulation, Nano letters 18, 6175 (2018).
  * Pokrovsky [1979] V. Pokrovsky, Properties of ordered, continuously degenerate systems, Advances in Physics 28, 595 (1979).