# Interplay between superconductivity and non-Fermi liquid at a quantum-
critical point in a metal.
V. The $\gamma$ model and its phase diagram. The case $\gamma=2$.

Yi-Ming Wu  School of Physics and Astronomy and William I. Fine Theoretical
Physics Institute, University of Minnesota, Minneapolis, MN 55455, USA  Shang-
Shun Zhang  School of Physics and Astronomy and William I. Fine Theoretical
Physics Institute, University of Minnesota, Minneapolis, MN 55455, USA  Artem
Abanov  Department of Physics, Texas A&M University, College Station, 77843
USA  Andrey V. Chubukov  School of Physics and Astronomy and William I. Fine
Theoretical Physics Institute, University of Minnesota, Minneapolis, MN 55455,
USA

###### Abstract

This paper is a continuation and a partial summary of our analysis of the
pairing at a quantum-critical point (QCP) in a metal for a set of quantum-
critical systems, whose low-energy physics is described by an effective model
with dynamical electron-electron interaction
$V(\Omega_{m})=({\bar{g}}/|\Omega_{m}|)^{\gamma}$ (the $\gamma$-model).
Examples include pairing at the onset of various spin and charge density-wave
and nematic orders and pairing in SYK-type models. In previous papers, we
analyzed the physics for $\gamma<2$. We have shown that the onset temperature
for the pairing $T_{p}$ is finite, of order ${\bar{g}}$, yet the gap equation
at $T=0$ has an infinite set of solutions within the same spatial symmetry. As
the consequence, the condensation energy $E_{c}$ has an infinite number of
minima. The spectrum of $E_{c}$ is discrete, but becomes more dense as
$\gamma$ increases. Here we consider the case $\gamma=2$. The $\gamma=2$ model
attracted special interest in the past as it describes the pairing by an
Einstein phonon in the limit when the dressed phonon mass $\omega_{D}$
vanishes. We show that for $\gamma=2$, the spectrum of $E_{c}$ becomes
continuous. We argue that the associated gapless ”longitudinal” fluctuations
destroy superconducting phase coherence at a finite $T$, such that at
$0<T<T_{p}$ the system displays pseudogap behavior of preformed pairs. We show
that for each gap function from the continuum spectrum, there is an infinite
array of dynamical vortices in the upper half-plane of frequency. For the
electron-phonon case, our results show that $T_{p}=0.1827{\bar{g}}$, obtained
in earlier studies, marks the onset of the pseudogap behavior, while the
actual superconducting $T_{c}$ vanishes at $\omega_{D}\to 0$.

##  I  Introduction.

This work presents a continuation and a partial summary of our analysis of the
interplay between non-Fermi liquid (NFL) physics and superconductivity for a
set of quantum-critical (QC) itinerant systems, whose low-energy dynamics can
be described by an effective model of spin-full electrons with dynamical
interaction on the Matsubara axis
$V(\Omega_{m})=({\bar{g}}/|\Omega_{m}|)^{\gamma}$ (the $\gamma$-model).
Examples include pairing near spin-density-wave, charge-density wave, and
Ising-nematic instabilities in isotropic and anisotropic 3D and 2D systems,
pairing of fermions at a half-filled Landau level, pairing of dispersion-less
fermions randomly coupled to phonons, and so on. We listed the examples in the
first paper in the series, Ref. [1], and discussed earlier works. In that and
subsequent papers [1, 2, 3, 4], hereafter called Papers I-IV, we analyzed the
$\gamma$-model with exponents $0<\gamma<2$. We rationalized Eliashberg-type
approach, solved generalized Eliashberg equations, and found that the solution
with a non-zero pairing gap develops below a finite $T_{p}$, which for a
generic $\gamma=O(1)$ is of order ${\bar{g}}$. The corresponding gap function
$\Delta(k,\omega_{m})$ can be roughly approximated as
$\Delta(\omega_{m})f(k)$, where a normalized $f(k)$ has a particular spatial
symmetry ($d-$wave, $s^{+-}$, etc), and $\Delta(\omega_{m})$ is a sign-
preserving function of Matsubara frequency, whose amplitude increases with
decreasing $T$. At $T=0$, $\Delta(0)$ is of order $T_{p}$ (the ratio
$\Delta(0)/T_{p}$ is a $\gamma-$dependent number [5, 6]). In this respect, the
pairing at a QCP is similar to pairing away from a QCP, when $V(\Omega_{m})$
saturates at a finite value at $\Omega_{m}=0$. However, on a more careful
look, we found qualitative difference between the two cases. Namely, away from
a QCP, Eliashberg gap equation at $T=0$ has at most a finite number of
solutions with a given spatial symmetry. At a QCP, it has an infinite number
of solutions for $\Delta(\omega)$, i.e., the condensation energy, $E_{c}$, has
an infinite number of local minima. The solutions $\Delta_{n}(\omega_{m})$,
labeled by an integer $n$, are topologically distinct in the sense that
$\Delta_{n}(\omega_{m})$ changes sign $n$ times along the positive Matsubara
axis, and each such point is a vortex in the upper half-plane of frequency.

The ultimate goal of our studies of the $\gamma-$model is to understand how
the existence of an infinite set of solutions affects the interplay between
pairing (i.e., the appearance of a non-zero $\Delta(\omega)$) and a true
superconductivity. In a conventional Eliashberg theory (ET) of SC out of a
non-critical Fermi liquid, phase fluctuations are small in the same parameter
that allows one to neglect vertex corrections. In this situation, the onset
temperature for the pairing, $T_{p}$, and the actual superconducting $T_{c}$
nearly coincide. To address a possible reduction of superconducting $T_{c}$ by
phase fluctuations, one then has to include the effects not incorporated into
ET, e.g., a localization of electrons when the interaction exceeds the
fermionic bandwidth. We analyze whether the emergence of an infinite number of
solutions for the gap at a QCP gives rise to a substantial reduction of
$T_{c}/T_{p}$ ratio already when the interaction is smaller than the fermionic
bandwidth. For $\gamma<2$, which we analyzed in Papers I-IV, the spectrum of
the condensation energies, $E_{c,n}$, is discrete, and $E_{c,0}$ is the
largest. In this situation, the physics at small $T$ is determined by a single
solution $\Delta_{0}(\omega_{m})$, and phase fluctuations are weak. We showed,
however, that the set becomes more dense at $\gamma$ increases: the ratios
$(E_{c,n+1}-E_{c,n})/E_{c,n}$ get progressively smaller.

In this paper, we analyze the case $\gamma=2$. We argue that for this
$\gamma$, the set of the gap functions and the spectrum of the condensation
energies become continuous. Specifically, at $\gamma=2-\delta$, and
$\delta=0+$, $\Delta_{n}(\omega_{m})$ with $n<1/\delta$ become equal to
$\Delta_{0}(\omega_{m})$ for all $\omega_{m}>0$, and $E_{c,n}$ become equal to
$E_{c,0}$, while $\Delta_{n}(\omega_{m})$ and $E_{c,n}$ with infinite
$n>1/\delta$ form a continuous, one-parameter gapless spectra,
$\Delta_{\xi}(\omega_{m})$ $E_{c,\xi}$. Here $\xi$ is a continuous variable,
that runs between zero and infinity and depends on how the double limit
$n\to\infty$ and $\delta\to 0$ is taken (we define $\xi$ such that the minimum
of $E_{c,\xi}$ is at $\xi=0$, and $E_{c,\infty}=0$). We argue that fluctuation
corrections to superconducting order parameter from $E_{c,\xi}$ destroy long-
range superconducting order at any finite $T$. We emphasize that this holds
for itinerant fermions, in the limit when the interaction is smaller than the
bandwidth, and the ET is rigorously justified.

We present a corroborative evidence that the $\gamma=2$ model is critical. It
comes from the analysis of the gap equation on the real frequency axis and in
the upper half-plane of frequency,
$z=\omega^{\prime}+i\omega^{{}^{\prime\prime}}$,
$\omega^{{}^{\prime\prime}}>0$. The gap function $\Delta(z)$ generally cannot
be obtained from $\Delta(\omega_{m})$ by just replacing $i\omega_{m}$ by $z$
as such gap function is not guaranteed to be analytic. To obtain an analytic
function, one has to perform a more sophisticated analysis [7, 8, 9, 10]. As a
consequence, $\Delta(\omega)$ on the real axis can be quite different from
$\Delta(\omega_{m})$. For the $\gamma$-model, some difference is expected on
general grounds, particularly for $\gamma>1$, because while the interaction
$V(\Omega_{m})$ on the Matsubara axis is positive (attractive) for all
$\gamma$, the one on the real axis is complex,
$V(\Omega)=e^{i\pi\gamma/2~{}{\text{sgn}}\Omega}/|\Omega|^{\gamma}$, and its
real part $V^{\prime}(\Omega)=({\bar{g}}/|\Omega|)^{\gamma}\cos{\pi\gamma/2}$
becomes repulsive for $\gamma>1$.

In Paper IV we compared the forms of $\Delta_{n}(\omega_{m})$ and
$\Delta_{n}(\omega)$ for $1<\gamma<2$. We found that at small
$\omega,\omega_{m}<{\bar{g}}$ and at large
$\omega,\omega_{m}>{\bar{g}}(|\log{(2-\gamma)}|/(2-\gamma))^{1/2}$, the two
gap functions transfer into each other under a rotation
$i\omega_{m}\to\omega$. However, at intermediate
${\bar{g}}<\omega,\omega_{m}<(|\log{(2-\gamma)}|/(2-\gamma))^{1/2}$,
$\Delta_{n}(\omega_{m})=a_{n}/|\omega_{m}|^{\gamma}$ is a sign-preserving
function of frequency, while
$\Delta_{n}(\omega)=|\Delta_{n}(\omega)|e^{i\eta_{n}(\omega)}$ oscillates, and
its phase $\eta_{n}(\omega)$ winds up by $2\pi k_{\gamma}$, where $k_{\gamma}$
is an integer, which depends on $\gamma$, but not on $n$. We extended the
analysis to complex $z$ in the upper frequency half-plane and showed that
there exists an array of $k$ dynamical vortices, centered at some complex
$z_{i}$.

Here we show that for $\gamma=2$, the gap functions on the real axis form a
continuous set, each $\Delta_{\xi}(\omega)$ oscillates up to an infinite
frequency, and its phase winds up by an infinite number of $2\pi$.
Accordingly, the number of vortices at $z_{i}$ becomes infinite, and the array
of $z_{i}$ extends up to an infinite frequency, where, we argue, each
$\Delta_{\xi}(\omega)$ develops an essential singularity. We show that for
each gap function from the set, the density of states (DoS)
$N_{\xi}(\omega+i0)$ has an infinite number of maxima and minima, and does not
recover the normal state form up to $\omega=\infty$. For the solution with
$\xi=0$, which was studied before [7, 8, 9, 11, 12], $N(\omega+i0)$ reduces to
a set of $\delta-$functions at some $\omega_{i}$.

We combine the results for $\gamma=2$ and earlier results for $\gamma<2$
(Papers I-IV) and present the phase diagram of the $\gamma-$model for
$\gamma\leq 2$, Fig.16. For all $\gamma$, the ground state is a superconductor
with a finite superfluid stiffness $\rho_{s}$, and the onset temperature for
the pairing, $T_{p}$, is finite. However, superconducting $T_{c}$ decreases
with $\gamma$ and vanishes for $\gamma=2$. In between $T_{p}$ and $T_{c}$, the
system displays a pseudogap (preformed pairs) behavior. One feature of this
phase is ”gap filling” behavior, as $T$ increases towards $T_{p}$. In the next
paper we consider the case $\gamma>2$. We show that the behavior at a finite
$T$ remains largely the same as for $\gamma=2$, however new physics emerges at
$T=0$ and gives rise to a reduction and eventual vanishing of $\rho_{s}$ even
in the ground state.

The model with the pairing interaction
$V(\Omega_{m})=({\bar{g}}/|\Omega_{m}|)^{2}$ attracted a substantial attention
on its own as it describes the pairing, mediated by an Einstein boson, in the
limit where the effective (dressed) Debye frequency $\omega_{D}$ vanishes
111The model also describes strong coupling limit of the interaction between
dispersion-less fermions and phonons (SYK-Yukawa model), Refs. [11, 12, 10,
26, 60]. Electron-phonon model at $\omega_{D}\to 0$ has been studied before by
a large number of authors [14, *Scalapino_69, 16, 17, 18, 7, 8, 9]. We use the
results of these studies, particularly the works by Karakozov, Maksimov, and
Mikhailovsky [7], Marsiglio and Carbotte [8], and Combescot [9] as the input
for some of our calculations. This limit is often termed strong coupling as
the dimensionless coupling constant $\lambda=({\bar{g}}/\omega_{D})^{2}$
diverges at $\omega_{D}\to 0$. However, the interaction ${\bar{g}}$ is still
assumed to be smaller than the Fermi energy, $E_{F}$. Indeed, ET includes
contributions to all orders in $\lambda$ within the ladder approximation, but
neglects vertex corrections to ladder series. The latter hold in powers of
Migdal-Eliashberg parameter
$\lambda_{E}={\bar{g}}^{2}N_{0}/\omega_{D}=\lambda(N_{0}\omega_{D})$, where
$N_{0}\sim 1/E_{F}$ in the DoS per unit volume. For small enough
${\bar{g}}/E_{F}$, $\lambda_{E}$ remains small even when $\lambda$ is large.
From this perspective, the strong coupling limit of the ET is the double limit
in which $\omega_{D}$ and ${\bar{g}}/E_{F}$ tend to zero simultaneously, such
that $\lambda_{E}$ remains small. In physical terms, the smallness of
$\lambda_{E}\ll\lambda$ comes about because in a process that gives rise to a
vertex correction, fermions are forced to vibrate at a phonon frequency, far
away from their own resonance, while in the processes, which form series in
$\lambda$, fermions are vibrating near their resonance frequencies. The
smallness of $\lambda_{E}$ also allows one to neglect the renormalization of
the bosonic propagator by fermions, both in the normal and in the
superconducting state.

Previous studies have found that a non-zero gap function emerges at
$T_{p}\approx 0.25\omega_{D}e^{-1/\lambda}$ at weak coupling (Refs. [19, 20,
21, 22, 23, 24, 25]) and at $T_{p}=0.1827{\bar{g}}$ at strong coupling [16, 8,
9, 26]. 222This formula was originally obtained semi-analytically by Allen and
Dynes [16]. They expressed it as $T_{p}\sim\omega_{D}\sqrt{\lambda}$ to
emphasize that at strong coupling $T_{p}$ becomes larger than $\omega_{D}$.
Given that $\lambda=({\bar{g}}/\omega_{D})^{2}$, their formula reduces to
$T_{p}\sim{\bar{g}}$.. To understand the interplay between the onset of
pairing and $T_{c}$, one has to also compute superfluid stiffness, $\rho_{s}$.
At weak coupling, $\rho_{s}\sim E_{F}\gg T_{p}$ [28, 29, 30, 31, 32, 33]. In
this situation, $T_{c}$ and $T_{p}$ almost coincide. At strong coupling, the
situation is more complex. At $T=0$, the $\rho_{s}\sim T_{p}/\lambda_{E}$ (see
below). Within the validity of ET, this stiffness exceeds $T_{p}$. If we were
to neglect the continuum spectrum of the condensation energy, we would obtain
that $T_{p}$ and $T_{c}$ again also coincide, as thermal corrections to SC
order parameter are of order $T/\rho_{s}$ and hence remain small for all
$T<T_{p}$. Including the additional corrections from the continuum of
$E_{c,\xi}$, we find that thermal corrections actually hold in powers of
$T/(\omega_{D}\lambda_{E})$ and become of order one at
$T\sim\omega_{D}\lambda_{E}$, which we identify with the actual $T_{c}$. At
small $\omega_{D}/{\bar{g}}$, this $T_{c}$ is much smaller than $T_{p}$, even
if we set $\lambda_{E}=O(1)$. In between $T=T_{p}$ and $T_{c}$ the system
displays a preformed pairs behavior. When $\omega_{D}$ increases and becomes
of order ${\bar{g}}$, the pseudogap region shrinks and the system gradually
recovers BCS-like behavior (Fig.15).

The structure of the paper is the following. In Sec.II we present the
Eliashberg gap equations that we use in this paper. In Sec. III we discuss the
solution of the gap equation along the Matsubara axis at $T=0$ and $\gamma\to
2$. We first obtain, in Secs. III.1 and III.2, the exact solution of the
linearized gap equation, $\Delta_{\infty}(\omega_{m})$, which changes sign an
infinite number of times between $\omega_{m}=0$ and $\omega_{m}\sim{\bar{g}}$,
and sign-preserving solution $\Delta_{0}(\omega_{m})$, which tends to a finite
value at $\omega_{m}\to 0$. At larger $\omega_{m}>{\bar{g}}$, both
$\Delta_{\infty}(\omega_{m})$ and $\Delta_{0}(\omega_{m})$ scale as
$1/|\omega_{m}|^{2}$. In Sec. III.3, we obtain the solutions of the non-linear
gap equation in the order-by-order expansion in the gap magnitude and show
that they form a one-parameter continuum set $\Delta_{\xi}(\omega)$, for which
$\Delta_{\infty}(\omega_{m})$ and $\Delta_{0}(\omega_{m})$ are the two
limiting cases. In Sec. IV, we analyze the properties of the gap function
$\Delta(\omega)$ along the real frequency axis. We first obtain, in Sec. IV.1,
the exact solution of the linearized gap equation on the real axis,
$\Delta_{\infty}(\omega)$, and show that it oscillates not only at
$\omega<{\bar{g}}$, but also at $\omega>{\bar{g}}$, with a different period.
In Sec. IV.2 we consider the real-frequency form of $\Delta_{0}(\omega)$,
which does not change sign on the Matsubara axis. We use as an input the
results from earlier works [7, 8, 9], which demonstrated that
$\Delta_{0}(\omega)=|\Delta_{0}(\omega)|e^{i\eta(\omega)}$ oscillates at
$\omega>{\bar{g}}$, and argue that the phase $\eta(\omega)$ winds up by an
infinite number of $2\pi$ between $\omega=O({\bar{g}})$ and $\omega=\infty$.
In Sec. IV.3 we present a one-parameter continuum set of
$\Delta_{\xi}(\omega)$, which in the two limits reduces to
$\Delta_{\infty}(\omega)$ and $\Delta_{0}(\omega)$. In Sec. V we extend
$\Delta_{\xi}(\omega)$ into the upper frequency half-plane ($\omega\to z$) and
show that for each $\xi$, there is an infinite array of vortices in the upper
frequency half-plane and an essential singularity at $|z|=\infty$. In Sec. VI
we consider the gap equation at a finite $\omega_{D}$. We argue that the
number of vortices becomes finite and the high-frequency behavior of the gap
function becomes regular; however, this holds only above a frequency, which
scales inversely with $\omega_{D}$. In Sec. VII we consider fluctuation
corrections to superconducting order parameter $\Delta\langle
e^{i\eta}\rangle$ We argue that the ground state is a superconductor, however
corrections to $\langle e^{i\eta}\rangle$ become $O(1)$ already at
$T\leq\omega_{D}$. We identify this scale with the actual superconducting
$T_{c}$ and discus pseudogap behavior in between $T_{p}\sim{\bar{g}}$ and
$T_{c}$. In Sec. VIII we combine the results for $\gamma=2$ and for $\gamma<2$
from Papers I-IV and obtain the full phase diagram of the $\gamma$ model for
$\gamma\leq 2$. We present our conclusions in Sec. IX. Some technical aspects
are discussed in the Appendices.

##  II Eliashberg equations

The Eliashberg gap equation for the $\gamma$-model is obtained by combining
the equations for the pairing vertex $\Phi$ and the self-energy $\Sigma$. The
two equations are obtained in a standard way, by summing up ladder series and
neglecting vertex corrections (see Paper I and the text below for
justification). On the Matsubara axis we have
($\Phi=\Phi(\omega_{m}),\Sigma=\Sigma(\omega_{m})$):

$\displaystyle\Phi(\omega_{m})={\bar{g}}^{\gamma}\pi
T\sum_{m^{\prime}}\frac{\Phi(\omega_{m^{\prime}})}{\sqrt{(\omega_{m^{\prime}}+\Sigma(\omega_{m^{\prime}}))^{2}+\Phi^{2}(\omega_{m^{\prime}})}}~{}\frac{1}{(|\omega_{m}-\omega_{m^{\prime}}|^{2}+\omega^{2}_{D})^{\gamma/2}}$
$\displaystyle\Sigma(\omega_{m})={\bar{g}}^{\gamma}\pi
T\sum_{m^{\prime}}\frac{\omega_{m^{\prime}}+\Sigma(\omega_{m^{\prime}})}{\sqrt{(\omega_{m^{\prime}}+\Sigma(\omega_{m^{\prime}}))^{2}+\Phi^{2}(\omega_{m^{\prime}})}}~{}\frac{1}{(|\omega_{m}-\omega_{m^{\prime}}|^{2}+\omega^{2}_{D})^{\gamma/2}}$
(1)

Introducing
$\Delta(\omega_{m})=\Phi(\omega_{m})\omega_{m}/(\omega_{m}+\Sigma(\omega_{m}))$
and substituting into (1), we obtain after a simple algebra the equation that
contains only $\Delta(\omega_{m})$:

$\Delta(\omega_{m})={\bar{g}}^{\gamma}\pi
T\sum_{m^{\prime}}\frac{\Delta(\omega_{m^{\prime}})-\Delta(\omega_{m})\frac{\omega_{m^{\prime}}}{\omega_{m}}}{\sqrt{(\omega_{m^{\prime}})^{2}+\Delta^{2}(\omega_{m^{\prime}})}}~{}\frac{1}{(|\omega_{m}-\omega_{m^{\prime}}|^{2}+\omega^{2}_{D})^{\gamma/2}}.$
(2)

For $\gamma=2$ this reduces to

$\Delta(\omega_{m})={\bar{g}}^{2}\pi
T\sum_{m^{\prime}}\frac{\Delta(\omega_{m^{\prime}})-\Delta(\omega_{m})\frac{\omega_{m^{\prime}}}{\omega_{m}}}{\sqrt{(\omega_{m^{\prime}})^{2}+\Delta^{2}(\omega_{m^{\prime}})}}~{}\frac{1}{|\omega_{m}-\omega_{m^{\prime}}|^{2}+\omega^{2}_{D}}.$
(3)

This is the same equation as for the interaction with an Einstein phonon [14,
*Scalapino_69, 34, 35, 16, 17, 18, 36, 7, 8, 9, 11, 12, 26]. As we said, we
consider the limit $\omega_{D}\to 0$. The self-action term with $m^{\prime}=m$
in the r.h.s. of (3) can be safely eliminated because the numerator vanishes
at $m=m^{\prime}$. Setting then $\omega_{D}=0$, we obtain the gap equation at
a QCP:

$\Delta(\omega_{m})={\bar{g}}^{2}\pi T\sum_{m^{\prime}\neq
m}\frac{\Delta(\omega_{m^{\prime}})-\Delta(\omega_{m})\frac{\omega_{m^{\prime}}}{\omega_{m}}}{\sqrt{(\omega_{m^{\prime}})^{2}+\Delta^{2}(\omega_{m^{\prime}})}}~{}\frac{1}{(\omega_{m}-\omega_{m^{\prime}})^{2}}.$
(4)

At $T=0$, $\pi T\sum_{m^{\prime}\neq m}\to(1/2)\int d\omega^{\prime}_{m}$.

The gap equation on the real axis is more conveniently expressed in terms of

$D(\omega)=\Delta(\omega)/\omega.$ (5)

The equation has the form [9]

$D(\omega)\omega B(\omega)=A(\omega)+C(\omega)$ (6)

where

$\displaystyle
A(\omega)=-\frac{{\bar{g}}^{2}}{2}\int_{0}^{\infty}d\omega^{\prime}\Re\frac{D(\omega^{\prime})}{\sqrt{1-D^{2}(\omega^{\prime})}}A_{T}$
$\displaystyle
A_{T}=\frac{\tanh{\frac{\omega^{\prime}}{2T}}+\tanh{\frac{\omega}{2T}}}{(\omega^{\prime}+\omega)^{2}}+\frac{\tanh{\frac{\omega^{\prime}}{2T}}-\tanh{\frac{\omega}{2T}}}{(\omega^{\prime}-\omega)^{2}}-\frac{1}{T\cosh^{2}{\frac{\omega}{2T}}}~{}\frac{\omega^{\prime}}{(\omega^{\prime})^{2}-\omega^{2}}$
$\displaystyle
B(\omega)=1+\frac{{\bar{g}}^{2}}{2\omega}\int_{0}^{\infty}d\omega^{\prime}\left[\Re\frac{1}{\sqrt{1-D^{2}(\omega^{\prime})}}\right]B_{T}$
$\displaystyle
B_{T}=\frac{\tanh{\frac{\omega^{\prime}}{2T}}+\tanh{\frac{\omega}{2T}}}{(\omega^{\prime}+\omega)^{2}}-\frac{\tanh{\frac{\omega^{\prime}}{2T}}-\tanh{\frac{\omega}{2T}}}{(\omega^{\prime}-\omega)^{2}}+\frac{1}{T\cosh^{2}{\frac{\omega}{2T}}}~{}\frac{\omega}{(\omega^{\prime})^{2}-\omega^{2}}$
$\displaystyle
C(\omega)=-i\frac{{\bar{g}}^{2}\pi}{2\sqrt{1-D^{2}(\omega)}}~{}~{}\left[\frac{dD(\omega)}{d\omega}\tanh{\frac{\omega}{2T}}-T~{}\left(\frac{dD^{2}(\omega)}{d\omega^{2}}+\left(\frac{dD(\omega)}{d\omega}\right)^{2}~{}\frac{D(\omega)}{1-D^{2}(\omega)}\right)\right]$
(7)

where the integrals are principal values. At $T=0$, the expressions simplify
to

$\displaystyle A(\omega)$ $\displaystyle=$
$\displaystyle-{\bar{g}}^{2}\int_{0}^{\infty}\frac{d\omega^{\prime}}{(|\omega|+\omega^{\prime})^{2}}\Re\frac{D(\omega^{\prime})}{\sqrt{1-D^{2}(\omega^{\prime})}}$
$\displaystyle B(\omega)$ $\displaystyle=$ $\displaystyle
1+\frac{{\bar{g}}^{2}}{|\omega|}\int_{0}^{\infty}\frac{d\omega^{\prime}}{(|\omega|+\omega^{\prime})^{2}}\Re\frac{1}{\sqrt{1-D^{2}(\omega^{\prime})}}$
$\displaystyle C(\omega)$ $\displaystyle=$
$\displaystyle-i\frac{\pi{\bar{g}}^{2}}{2\sqrt{1-D^{2}(\omega)}}~{}\frac{dD(\omega)}{d\omega}~{}\mbox{sign}\omega,$
(8)

and the gap equation becomes

$\displaystyle-i\frac{\pi{\bar{g}}^{2}}{2}~{}\frac{\frac{dD(\omega)}{d\omega}}{\sqrt{1-D^{2}(\omega)}}~{}\mbox{sign}\omega=$
(9) $\displaystyle
D(\omega)\omega~{}\left(1+\frac{{\bar{g}}^{2}}{|\omega|}\int_{0}^{\infty}\frac{d\omega^{\prime}}{(|\omega|+\omega^{\prime})^{2}}\Re\frac{1}{\sqrt{1-D^{2}(\omega^{\prime})}}\right)+\int_{0}^{\infty}\frac{d\omega^{\prime}}{(|\omega|+\omega^{\prime})^{2}}\Re\frac{D(\omega^{\prime})}{\sqrt{1-D^{2}(\omega^{\prime})}}$

The functions $A(\omega)$ and $B(\omega)$ can be equivalently expressed in
terms of the solution of the gap equation on the Matsubara axis [36, 8, 9]:

$\displaystyle A(\omega)$ $\displaystyle=$ $\displaystyle 2\pi
T\sum_{m=0}^{\infty}\frac{\Delta(\omega_{m})}{\sqrt{\Omega^{2}_{m}+\Delta^{2}(\omega_{m})}}\frac{\omega^{2}_{m}-\omega^{2}}{(\omega^{2}_{m}+\omega^{2})^{2}}$
$\displaystyle B(\omega)$ $\displaystyle=$ $\displaystyle 1+4\pi
T\sum_{m=0}^{\infty}\frac{1}{\sqrt{\Omega^{2}_{m}+\Delta^{2}(\omega_{m})}}\frac{\omega^{2}_{m}}{(\omega^{2}_{m}+\omega^{2})^{2}}.$
(10)

This simplifies numerical calculations: the recipe is to first solve for the
gap at the Matsubata points $\omega_{m}=\pi T(2m+1)$ and then use Eqs. (10) as
an input for the calculation of $D(\omega)$ on the real axis.

##  III Gap equation along the Matsubara axis at $T=0$

In Papers I-IV we analyzed the gap equation for $\gamma<2$ and found that at
$T=0$ it has an infinite, discrete set of solutions at
$\Delta(\omega_{m})=\Delta_{n}(\omega_{m})$. A gap function
$\Delta_{n}(\omega_{m})$ changes sign $n$ times between $\omega_{m}=0$ and
$\omega_{m}=O({\bar{g}})$ and decays as $1/|\omega_{m}|^{\gamma}$ at larger
frequencies. The two end points of the set are the sign-preserving solution
$\Delta_{0}(\omega_{m})$ and the solution of the linearized gap equation
$\Delta_{\infty}(\omega_{m})$, which changes sign an infinite number of times.
The existence of this infinite set is a distinct feature of the pairing at a
QCP. Away from a QCP, the number of solutions becomes finite
$(n=0,1..n_{max}$), and far away from a QCP only the $n=0$ solution remains,
like in a conventional BCS/Eliashberg theory.

Here we extend this analysis to $\gamma=2$. We show that for this $\gamma$,
the set of gap functions becomes $\Delta_{\xi}(\omega_{m})$, where
$0\leq\xi\leq\infty$ is a continuous variable. We first analyze the two end
points, $\Delta_{\infty}(\omega_{m})$ and $\Delta_{0}(\omega_{m})$, and then
obtain the gap function for arbitrary $\xi$.

###  III.1 Linearized gap equation

The linearized gap equation at $T=0$ is obtained from (4) by assuming that the
gap function is infinitesimally small,
$\Delta(\omega_{m})=\Delta_{\infty}(\omega_{m})$. In terms of
$D_{\infty}(\omega_{m})=\Delta_{\infty}(\omega_{m})/\omega_{m}$ we have

$D_{\infty}(\omega_{m})=\frac{{\bar{g}}^{2}}{2\omega_{m}}\int
d\omega^{\prime}_{m}\frac{D_{\infty}(\omega^{\prime}_{m})-D_{\infty}(\omega_{m})}{(\omega_{m}-\omega_{m^{\prime}})^{2}}~{}{\text{sign}}\omega^{\prime}_{m}.$
(11)

One can verify that the leading term in $D_{\infty}(\omega_{m})$ at small
$\omega_{m}\ll{\bar{g}}$ is obtained by neglecting the l.h.s. of (11), i.e.,
by solving

$\int_{0}^{\infty}d\omega^{\prime}_{m}\left[\frac{D_{\infty}(\omega^{\prime}_{m})-D_{\infty}(\omega_{m})}{(\omega_{m}-\omega_{m^{\prime}})^{2}}+\frac{D_{\infty}(\omega^{\prime}_{m})+D_{\infty}(\omega_{m})}{(\omega_{m}+\omega_{m^{\prime}})^{2}}\right]=0$
(12)

This approximation is equivalent to neglecting the bare $\omega$ in the
fermionic propagator in comparison with the NFL fermionic self-energy without
the self-action term, $\Sigma(\omega_{m})=-{\bar{g}}^{2}/\omega_{m}$.

The solution of (12) is

$D_{\infty}(\omega_{m})=2\epsilon~{}{\text{Re}}\left[e^{i(\beta\log{\left(\frac{|\omega_{m}|}{{\bar{g}}}\right)^{2}}+\phi)}\right]~{}\mbox{sign}\omega=2\epsilon\cos{\left(\beta\log{\left(\frac{|\omega_{m}|}{{\bar{g}}}\right)^{2}}+\phi\right)}~{}\mbox{sign}\omega$
(13)

where $\epsilon$ is an infinitesimally small real overall factor, $\phi$ is a
phase factor, which is arbitrary at this stage, and $\beta=0.38187$ satisfies
$\pi\beta\tanh(\pi\beta)=1$ and is the solution of

$\int_{-\infty}^{\infty}dx\frac{|x|^{2i\beta}-{\text{sign}}x}{(x-1)^{2}}=0,$
(14)

The function $D_{\infty}(\omega_{m})$ is scale-invariant (an arbitrary phase
factor $\phi$ can be absorbed into the prefactor for $\omega$ under the
logarithm). This is the consequence of the fact that ${\bar{g}}$ falls off
from the gap equation (11), once we neglect $D_{\infty}(\omega_{m})$ in the
l.h.s.

We now analyze the full gap equation. By power counting, the r.h.s of (11) is
of order $D_{\infty}(\omega_{m})({\bar{g}}/|\omega_{m}|)$. This justifies
neglecting $D_{\infty}(\omega_{m})$ in the l.h.s. for
$|\omega_{m}|<O({\bar{g}})$, but for larger frequencies it must be kept.

We obtained the exact solution of Eq. (11). The derivation parallels the one
for smaller $\gamma$ in Papers I and IV. We skip the details and present the
final result:

$D_{\infty}(\omega_{m})=\epsilon\frac{{\bar{g}}}{\omega_{m}}\int_{-\infty}^{\infty}dkb_{k}e^{-ik\log{(\omega_{m}/{\bar{g}})^{2}}},$
(15)

where

$b_{k}=\frac{e^{-iI_{k}}}{\left[\cosh(\pi(k-\beta))\cosh(\pi(k+\beta))\right]^{1/2}},$
(16)

and

$I_{k}=\frac{1}{2}\int_{-\infty}^{\infty}dk^{\prime}\log{|\epsilon_{k^{\prime}}-1|}\tanh{\pi(k^{\prime}-k+i\delta)},$
(17) $\epsilon_{k^{\prime}}=\pi k^{\prime}\tanh(\pi k^{\prime}),$ (18)

Here $\beta\simeq 0.38187$ is the same as in Eq. (13).

Figure 1: $D_{\infty}(\omega_{m})$ as a function of $\omega_{m}/{\bar{g}}$.
The scale is logarithmic for $\omega_{m}<{\bar{g}}$ and linear at
$\omega_{m}>{\bar{g}}$.

At $\omega_{m}\ll{\bar{g}}$, the exact $D_{\infty}(\omega_{m})$ has the form
of Eq. (13) with some particular $\phi$. At $\omega_{m}\gg{\bar{g}}$,
$D_{\infty}(\omega_{m})$ does not oscillate and decreases as
$1/(\omega_{m})^{3}$ ($\Delta_{\infty}(\omega_{m})$ decreases as
$1/(\omega_{m})^{2}$). We plot the exact $D_{\infty}(\omega_{m})$ in Fig. 1.
The crossover between the two forms occurs at $\omega_{m}\sim{\bar{g}}$, as
expected.

The corrections to Eq. (13) at small $\omega_{m}$ hold in powers of
$|\omega_{m}|/{\bar{g}}$; the leading correction scales as
$(|\omega_{m}|/{\bar{g}})^{2.78}$. The corrections to $1/(\omega_{m})^{3}$ at
large $\omega_{m}$ hold in powers of ${\bar{g}}/|\omega_{m}|$; the leading
correction scales as
$({\bar{g}}/\rvert\omega_{m}\rvert)^{5}\log(\rvert\omega_{m}\rvert/{\bar{g}})$.
We present the details of the analysis in Appendix B. There, we also show that
at $\omega_{m}\gg{\bar{g}}$ there exists an exponentially small, oscillating
component $D_{\infty;u}(\omega_{m})$ in the form

$D_{\infty;u}(\omega_{m})\propto 2\sqrt{2}\epsilon
e^{-(|\omega_{m}|/{\bar{g}})^{2}}\cos{\left[\frac{(\pi^{2}-2)}{2\pi}\left(\frac{|\omega_{m}|}{{\bar{g}}}\right)^{2}+\frac{\pi}{4}\right]}.$
(19)

This term is the contribution to $D_{\infty}$ from large $k$ and $k^{\prime}$
in Eqs. (16) and (17), It is completely irrelevant on the Matsubara axis, but
we will see that it gives the dominant contribution to $D_{\infty}(\omega)$ on
the real axis.

###  III.2 Non-linear gap equation. Sign-preserving solution.

We now analyse the full non-linear gap equation, Eq. (4). We first search for
a “conventional” sign-preserving solution $\Delta_{0}(\omega_{m})$

The analytical analysis uses the same computational steps as in Paper IV and
we will be brief. We use the identity

$\int_{-\infty}^{\infty}d\omega^{\prime}_{m}\frac{1-\frac{\omega^{\prime}_{m}}{\omega_{m}}}{|\omega_{m}-\omega^{\prime}_{m}|^{\gamma}}=0,$
(20)

valid for $\gamma>1$, and re-express Eq. (4) as

$\displaystyle\Delta_{0}(\omega_{m})\left[1-\frac{{\bar{g}}^{2}}{2}\int_{-\infty}^{\infty}d\omega^{\prime}_{m}\frac{1-\frac{\omega^{\prime}_{m}}{\omega_{m}}}{|\omega_{m}-\omega^{\prime}_{m}|^{2}}\left(\frac{1}{\sqrt{\Delta^{2}_{0}(\omega^{\prime}_{m})+(\omega^{\prime}_{m})^{2}}}-\frac{1}{\Delta_{0}(\omega_{m})}\right)\right]$
$\displaystyle=\frac{{\bar{g}}^{2}}{2}\int_{-\infty}^{\infty}d\omega^{\prime}_{m}\frac{\Delta_{0}(\omega^{\prime}_{m})-\Delta_{0}(\omega_{m})}{|\omega_{m}-\omega^{\prime}_{m}|^{2}\sqrt{\Delta^{2}_{0}(\omega^{\prime}_{m})+(\omega^{\prime}_{m})^{2}}}$
(21)

Both integrals in (21) are infra-red convergent and are determined by
$\omega_{m}^{\prime}\leq\Delta(\omega^{\prime}_{m})$. In the limit
$\omega_{m}\to 0$, Eq. (21) reduces to

$\displaystyle\Delta_{0}(0)\left[1-{\bar{g}}^{2}\int_{0}^{\infty}d\omega^{\prime}_{m}\frac{\sqrt{\Delta^{2}_{0}(\omega^{\prime}_{m})+(\omega^{\prime}_{m})^{2}}-\Delta_{0}(0)}{\Delta_{0}(0)\sqrt{\Delta^{2}_{0}(\omega^{\prime}_{m})+(\omega^{\prime}_{m})^{2}}|\omega^{\prime}_{m}|^{2}}\right]$
$\displaystyle={\bar{g}}^{2}\int_{0}^{\infty}\frac{d\omega^{\prime}_{m}~{}\left(\Delta_{0}(\omega^{\prime}_{m})-\Delta_{0}(0)\right)}{\sqrt{\Delta^{2}_{0}(\omega^{\prime}_{m})+(\omega^{\prime}_{m})^{2}}|\omega^{\prime}_{m}|^{2}}$
(22)

We assume and then verify that
$\Delta_{0}(\omega^{\prime}_{m})\approx\Delta_{0}(0)$ for
$\omega^{\prime}_{m}\leq\Delta_{0}(\omega^{\prime}_{m})$, relevant for both
integrals in (22). Substituting into (22), we find

$1\approx{\bar{g}}^{2}\int_{0}^{\infty}d\omega^{\prime}_{m}\frac{\sqrt{\Delta^{2}_{0}(0)+(\omega^{\prime}_{m})^{2}}-\Delta_{0}(0)}{\Delta_{0}(0)\sqrt{\Delta^{2}_{0}(0)+(\omega^{\prime}_{m})^{2}}|\omega^{\prime}_{m}|^{2}}$
(23)

The integral can be evaluated analytically and yields
$\Delta_{0}(0)={\bar{g}}$. Substituting further
$\Delta_{0}(\omega^{\prime}_{m})={\bar{g}}$ into the r.h.s. of (21), we find
that $\Delta_{0}(\omega_{m})$ varies quadratically with $\omega_{m}$ at small
$\omega_{m}$ and for $\omega_{m}\leq{\bar{g}}$ remains comparable to
$\Delta_{0}(0)$.

Figure 2: Sign-preserving solution $\Delta_{0}(\omega_{m})$ of the nonlinear
gap equation along the Matsubara axis. We obtained $\Delta_{0}(\omega_{m})$ by
solving the non-linear gap equation numerically and taking the limit $T\to 0$.
At $\omega_{m}<\bar{g}$, $\Delta_{0}(\omega_{m})$ remains comparable to
$\bar{g}$; at larger frequencies it decays as $1/\omega_{m}^{2}$.

In the opposite limit of large $\omega_{m}$, the prefactor for
$\Delta_{0}(\omega)$ in the l.h.s. of (21) is approximately $1$, and in the
r.h.s. of this equation $1/|\omega_{m}|^{2}$ can be pulled out from the
integral. This yields

$\Delta_{0}(\omega_{m})\approx
Q\left(\frac{{\bar{g}}}{|\omega_{m}|}\right)^{2}$ (24)

where

$Q=\int_{0}^{\infty}\frac{d\omega^{\prime}_{m}\Delta_{0}(\omega^{\prime}_{m})}{\sqrt{\Delta^{2}_{0}(\omega^{\prime}_{m})+(\omega^{\prime}_{m})^{2}}}$
(25)

The integral is determined by $\omega_{m}^{\prime}\sim{\bar{g}}$ and is of
order ${\bar{g}}$. Then
$\Delta_{0}(\omega_{m})\sim{\bar{g}}^{3}/|\omega_{m}|^{2}$ at high
frequencies. The full gap function is sign-preserving. We show the numerical
result for $\Delta_{0}(\omega_{m})$ in Fig.2. At small $\omega_{m}$ we find
$\Delta_{0}(0)\approx 0.75{\bar{g}}$. This fully agrees with the earlier
result, Ref. [8]. We note in passing that the first numerical evidence that
$\Delta_{0}(0)$ scales with ${\bar{g}}$ has been obtained in Ref. [frank_1].

###  III.3 Continuous set of solutions. Expansion in the gap magnitude

The solutions $\Delta_{\infty}(\omega_{m})$ and $\Delta_{0}(\omega_{m})$ (or,
equivalently, $D_{\infty}(\omega_{m})$ and $D_{0}(\omega_{m})$) also exist for
$\gamma<2$. For such $\gamma$, these two solutions are the end points of a
discrete set of topologically distinct solutions $\Delta_{n}(\omega)$. We
argue below that the set becomes continuous for $\gamma=2$. For a continuous
set there is no one-to-one correspondence between a particular member of the
set and integer $n$, and we will show how this correspondence gets lost at
$\gamma=2-0$.

Comparing $D_{\infty}(\omega_{m})$ and $D_{0}(\omega_{m})$, we see that they
have the same form $1/(\omega_{m})^{3}$ for $\omega_{m}>{\bar{g}}$, but are
very different for $\omega_{m}<{\bar{g}}$. We therefore focus on the range
$\omega_{m}<{\bar{g}}$ and use to our advantage the fact that we know the
analytic form of $D_{\infty}(\omega_{m})$ in this range, Eq. (13). We use this
$D_{\infty}(\omega_{m})$ as an input and expand it in powers of
$D^{2}(\omega^{\prime}_{m})$ in the r.h.s. of the gap equation (4).
Specifically, we will be searching for the solution of (4) in the form

$D(\omega_{m})=\sum_{j=0}^{\infty}D^{(2j+1)}(\omega_{m})$ (26)

where

$D^{(1)}(\omega_{m})=D_{\infty}(\omega_{m})=2\epsilon\cos{f(\omega_{m})}~{}{\text{\text{sign}}}\omega_{m}$
(27)

with

$f(\omega_{m})=\beta\log{\left(\frac{|\omega_{m}|}{\bar{g}}\right)^{2}}+\phi$
(28)

We will see that $D^{(2j+1)}\sim\epsilon^{2j+1}$.

Substituting $D(\omega_{m})$ from (26) into (4) and expanding in
$D^{2}(\omega^{\prime}_{m})$ in the r.h.s. of (4), we obtain the set of
equations, which express $D^{(2j+1)}$ for a given $j$ in terms of $D^{(2j+1)}$
with smaller $j$. For $j=1$, we have

$D^{(3)}(\omega_{m})\omega_{m}-\frac{{\bar{g}}^{2}}{2}\int
d\omega^{\prime}_{m}\left(D^{(3)}(\omega^{\prime}_{m})-D^{(3)}(\omega_{m})\right)\frac{\text{sign}{\omega^{\prime}_{m}}}{(\omega_{m}-\omega^{\prime}_{m})^{2}}=K_{3}(\omega_{m})$
(29)

where the source term is

$K_{3}(\omega_{m})=-\frac{{\bar{g}}^{2}}{4}\int
d\omega^{\prime}_{m}\left(D^{(1)}(\omega^{\prime}_{m})-D^{(1)}(\omega_{m})\right)[D^{(1)}(\omega^{\prime}_{m})]^{2}\frac{\text{sign}{\omega^{\prime}_{m}}}{(\omega_{m}-\omega^{\prime}_{m})^{2}}$
(30)

The source term is of order $\epsilon^{3}$, hence $D^{(3)}\propto\epsilon^{3}$
($D^{(5)}\propto\epsilon^{5}$ and so on). Substituting $D^{(1)}(\omega_{m})$
from Eq. (27) and evaluating the integrals, we find the source term for
$D^{(3)}$ as the sum of the two terms, $K_{3}=K_{3a}+K_{3b}$, where

$K_{3a}(\omega_{m})=-\epsilon^{3}\frac{\bar{g}}{\omega_{m}}\cos\left(3f(\omega)\right)\left(2\pi\beta\coth(2\pi\beta)-3\pi\beta\tanh(3\pi\beta)\right)~{}{\text{\text{sign}}}\omega_{m}$
(31)

and

$K_{3b}(\omega_{m})=-\epsilon^{3}\frac{\bar{g}}{\omega_{m}}\cos\left(f(\omega)\right)\frac{1+\sinh^{2}(\pi\beta)}{\sinh^{2}(\pi\beta)}~{}{\text{\text{sign}}}\omega_{m}$
(32)

Solving for $D^{(3)}$ we find that the first term gives rise to
$\epsilon^{3}\cos\left(3f(\omega)\right)$, while the second term accounts for
the renormalization of the prefactor for $\log(\omega^{2}_{m})$ in
$f(\omega_{m})$ in (28) To order $\epsilon^{2}$, the dressed $f(\omega_{m})$,
which we label $f_{\epsilon}(\omega_{m})$, becomes

$f_{\epsilon}(\omega_{m})=\beta_{\epsilon}\log{\left(\frac{|\omega_{m}|}{\bar{g}}\right)^{2}}+\phi_{\epsilon}$
(33)

where

$\beta_{\epsilon}=\beta\left(1-\epsilon^{2}/2\right)\approx\beta(1-\epsilon^{2})^{1/2}$
(34)

The full $D(\omega_{m})$ to order $\epsilon^{3}$ is

$D(\omega_{m})=2\left(\epsilon\cos{f_{\epsilon}(\omega_{m})}+Q_{3}\epsilon^{3}\cos{3f(\omega_{m})}\right)~{}{\text{\text{sign}}}\omega_{m}$
(35)

where

$Q_{3}=\frac{2\pi\beta\coth{(2\pi\beta)}-3\pi\beta\tanh{(3\pi\beta)}}{2(1-3\pi\beta\tanh{(3\pi\beta)})}=\frac{5-(\pi\beta)^{2}}{16}\approx
0.222$ (36)

Expanding to next order, we find (i) $\epsilon^{5}\cos{5f(\omega_{m})}$ term
with the prefactor $Q_{5}=0.043$, (ii) $O(\epsilon^{4})$ corrections to
$\beta_{\epsilon}$ in (34)
($\beta_{\epsilon}=1-0.5\epsilon^{2}+0.806\epsilon^{4})$), and (iii)
$O(\epsilon^{2})$ corrections to $Q_{3}$ ($Q_{3}\to Q_{3,\epsilon}$) and to
the argument of $\cos{3f(\omega_{m})}$ in (35). We verified that the last
correction changes $\cos{3f(\omega_{m})}$ to $\cos{3f_{\epsilon}(\omega_{m})}$
with the same $f_{\epsilon}$ as in (33). This is the strong indication that
the series contain the same fully renormalized $f_{\epsilon}(\omega_{m})$ in
each term. Combining the results, we obtain, for $\omega_{m}\ll{\bar{g}}$,
$D(\omega_{m})=D_{\epsilon}(\omega_{m})$, where

$D_{\epsilon}(\omega_{m})=2\epsilon\left(\cos{f_{\epsilon}(\omega_{m})}+Q_{3,\epsilon}\epsilon^{3}\cos{3f_{\epsilon}(\omega_{m})}+Q_{5,\epsilon}\epsilon^{5}\cos{5f_{\epsilon}(\omega_{m})}+...\right)~{}{\text{\text{sign}}}\omega_{m}$
(37)

We emphasize that a continuous set of solutions exists only for $\gamma=2$.
Applying the same perturbative analysis for $\gamma<2$, we find that the
expansion holds in $\epsilon^{2}({\bar{g}}/|\omega_{m}|)^{2-\gamma}$ and
breaks at a finite $\omega_{min}\sim{\bar{g}}\epsilon^{2/(2-\gamma)}$ (see
Appendix A for more detail). At smaller $\omega_{m}$, $\Delta(\omega_{m})$
saturates, and $D(\omega_{m})\propto 1/\omega_{m}$. The forms of
$D(\omega_{m})$ at $\omega_{m}<\omega_{min}$ and $\omega_{m}>\omega_{min}$
match only for a discrete set of $\epsilon=\epsilon_{n}$, which implies that
for $\gamma<2$ the solutions of the full non-linear gap equation form a
discrete set.

Because $f_{\epsilon}(\omega_{m})$ contains $\log{\omega^{2}_{m}}$, each
$D_{\epsilon}(\omega_{m})$ from (37) changes sign an infinite number of times
down to $\omega_{m}=0$, i.e., in our original classification the gap functions
from the set are different realizations of $n=\infty$. At $\omega_{m}=0$, each
$D_{\epsilon}(\omega_{m})$ has an essential singularity as neither
$\lim_{\omega_{m}\to 0}D_{\epsilon}(\omega_{m})$ nor $\lim_{\omega_{m}\to
0}1/D_{\epsilon}(\omega_{m})$ exist.

For a generic $\epsilon$, Eq. (37) is valid for $\omega_{m}<\bar{g}$. At
larger $\omega_{m}$, $D_{\epsilon}(\omega_{m})=D_{\epsilon}/|\omega_{m}|^{2}$.
We expect that for every $\epsilon$, the crossover to proper high-frequency
behavior can be achieved by fixing the phase factor $\phi_{\epsilon}$ in (33)
(see Paper I for a similar analysis for the linearized gap equation for
$\gamma<1$).

Next, we see from Eq. (34) that $\beta^{2}_{\epsilon}$ decreases with
increasing $\epsilon$. It is natural to expect that it vanishes at some
$\epsilon_{cr}=O(1)$. The expansion in (37) holds only as long as
$\beta_{\epsilon}$ is real, as there is no solution of the nonlinear gap
equation for imaginary $\beta_{\epsilon}$ (see Paper I for detailed discussion
on this). For $\epsilon\leq\epsilon_{cr}$, $\beta_{\epsilon}$ is small, and
the range, where $D(\omega_{m})$ oscillates, is confined to small
$\omega_{m}\leq{\bar{g}}e^{-\pi/\beta_{\epsilon}}$. By properly taking the
double limit $\epsilon\to\epsilon_{cr}$ and $\omega_{m}\to 0$, one can obtain
an infinite set of gap functions, which change sign a given number of times in
the immediate vicinity of $\omega_{m}=0$. At $\epsilon=\epsilon_{cr}$ all
these gap functions coincide with $\Delta_{0}(\omega_{m})$ at any
$\omega_{m}>0$. This agrees with the observation in Paper IV that as $\gamma$
increases towards $2$, the region, where $\Delta_{n}(\omega_{m})$ with finite
$n$ change sign, gets confined to progressively smaller $\omega_{m}$, while at
larger $\omega_{m}$, all $\Delta_{n}(\omega_{m})$ with $n=0,1,2...$ nearly
coincide. We illustrate this in Fig.3. For consistency with the notations in
previous sections, it is convenient to introduce
$\xi=(\epsilon_{cr}-\epsilon)/\epsilon$ and label the continuum set of the gap
functions by $\Delta_{\xi}(\omega_{m})$. Then the end point solutions
$\epsilon\to 0$ and $\epsilon=\epsilon_{cr}$ are $\Delta_{\infty}(\omega_{m})$
and $\Delta_{0}(\omega_{m})$.

Figure 3: The gap function $\Delta_{n}(\omega_{m})$ for $\gamma<2$ and
$\gamma=2$. For $\gamma<2$, $\Delta_{n}(\omega_{m})$ changes sign $n$ times.
As $\gamma$ gets close to $2$, the frequency region where these $n$ sign
changes happen, shrinks to progressively smaller $\omega_{m}=0$, and at
$\gamma=2-0$, $\Delta_{n}(\omega_{m})$ with finite $n$ collapse into
$\Delta_{0}(\omega_{m})$ at all $\omega_{m}>0$. The continuum set of
$\Delta_{\xi}(\omega_{m})$ at $\gamma=2-0$ emerges from
$\Delta_{n}(\omega_{m})$ with $n\to\infty$, and the continuous parameter $\xi$
is determined by how the double limit $n\to\infty$ and $\gamma\to 2$ is taken.
As the consequence, all $\Delta_{\xi}(\omega_{m})$ with $\xi>0$ change sign
infinite number of times between $\omega_{m}=0$ and $\omega_{m}\sim{\bar{g}}$.
The solution of the linearized gap equation is the $\xi\to\infty$ limit of
this set.

It is beyond the ability of the order-by-order expansion to determine the form
of $\Delta_{\xi}(\omega_{m})$ near $\xi=0$. On general grounds, we expect that
corrections to $f_{\epsilon}(\omega_{m})\to f_{\xi}(\omega_{m})$ in (33)
become relevant starting already from small frequencies, and that at $\xi=0$,
the gap function coincides with $\Delta_{0}(\omega_{m})$, which we found in
the previous section. A way to reproduce this behavior is to assume that at
$\xi\to 0$, the series for $D(\omega_{m})$ in (37) become geometrical
$[Q_{2n+1}\epsilon_{cr}^{2n+1}\approx(-1)^{n}]$. In this case,

$D_{\xi}(\omega_{m})\sim\frac{\cos{f_{\xi}(\omega_{m})}}{\alpha^{2}+\cos^{2}{f_{\xi}(\omega_{m})}}~{}{\text{\text{sign}}}\omega_{m},$
(38)

where $\alpha\sim\xi^{2}$ and
$f_{\xi}(\omega_{m})\sim\sqrt{\alpha}\log{{\bar{g}}/|\omega_{m}|}+f^{*}(\omega_{m})$,
where $f^{*}(\omega_{m})$ is a regular function of $\omega_{m}$, which at low
frequencies reduces to $\pi/2+O(\omega_{m}/{\bar{g}})$. For any $\xi>0$, this
$D_{\xi}(\omega_{m})$ changes sign an infinite number of times, but at
$\xi=0$, $D_{\xi=0}(\omega_{m})\sim{\bar{g}}/\omega_{m}$, as we expect. We
also note that between the nodes (the vortex points), $D_{\xi}(\omega_{m})$
from (38) is large, of order $1/\xi$. Extending this $D(\omega_{m})$ to
complex frequencies, $z=\omega^{\prime}+i\omega^{{}^{\prime\prime}}$, we find
that there exist anti-vortices at small $z$ in the lower frequency half-plane.
At $\xi=0$, vortices and anti-vortices annihilate at $z=0$, leaving a regular
gap function $\Delta_{0}(\omega_{m})$.

In Appendix C we consider the extended $\gamma-$ model with non-equal
interactions in the particle-particle and particle-hole channels and introduce
$M\neq 1$ as a measure of the difference of the two interactions. For the
extended model, there is a critical $M_{max}$, below which the ground state is
a non-Fermi liquid with $\Delta=0$. For $\gamma=2$, $M_{max}=0$. We obtain the
set of $\Delta_{\xi}(\omega_{m})$ at small $\omega_{m}$ at $M=0+$ and show
that all gap functions from the continuous set appear simultaneously with the
overall magnitude $M^{1/2}$.

We next analyze the condensation energy $E_{c}$. We define $E_{c}$ as the
difference between the actual ground state energy $E_{\Delta}$ at a finite
$\Delta(\omega_{m})$ and the would be ground state energy of the normal state,
$E_{\Delta=0}$. The expression for $E_{c}$ for $\gamma=2$ has been obtained
before: [35, 37, 38, 39] and we just copy it here:

$\displaystyle
E_{c}=-N_{0}\int_{0}^{\infty}d\omega_{m}\omega_{m}\frac{(\sqrt{1+D^{2}(\omega_{m})}-1)^{2}}{\sqrt{1+D^{2}(\omega_{m})}}$
$\displaystyle-
N_{0}{\bar{g}}^{2}\int_{0}^{\infty}d\omega_{m}d\omega^{\prime}_{m}\frac{\left(\sqrt{1+D^{2}(\omega_{m})}-{\sqrt{1+D^{2}(\omega^{\prime}_{m})}}\right)^{2}}{\sqrt{1+D^{2}(\omega_{m})}\sqrt{1+D^{2}(\omega^{\prime}_{m})}}~{}\frac{\omega_{m}\omega^{\prime}_{m}}{(\omega^{2}_{m}-(\omega^{\prime}_{m})^{2})^{2}}$
(39)

This formula has been derived with the use of (4) and is therefore valid only
for the solutions of the gap equation. Both terms in (39) are negative, i.e.
any solution of the gap equation lowers the ground state energy compared to
the normal state.

Substituting (37) into (39), we find that $E_{c}=E_{c,\xi}$ is a continuous
function of $\xi$. At $\xi\gg 1$,

$E_{c,\xi}=-aN_{0}\frac{{\bar{g}}^{2}}{\xi^{4}}$ (40)

where $a=O(1)$. It is natural to expect that $|E_{c,\xi}|$ increases with
decreasing $\xi$ and reaches a maximum at $\xi=0$, see Fig. 4. 333The second
term in (39) diverges logarithmically at $\xi=0$ if we use
$D_{0}(\omega)\approx\Delta(0)/\omega$ at small frequencies. This divergence
comes from the putative normal state energy $E_{\Delta=0}$ while the ground
state energy $E_{\Delta}$ remains finite. For $\xi>0$, both $E_{\Delta=0}$ and
$E_{\Delta}$ have logarithmic singularities, which cancel out in
$E_{c}=E_{\Delta}-E_{\Delta=0}$.

Figure 4: (a) The condensation energy $E_{c}$ the solutions of the Eliashberg
gap equation for $\gamma<2$. $E_{c}=E_{c,n}$ is a discrete function of a
number of a solution, $n$. The largest condensation energy is for $n=0$. (b)
The condensation energy $E_{c,\xi}$ for $\gamma=2$. $E_{c,\xi}$ is a
continuous function of the parameter $\xi$. The condensation energy at $\xi=0$
is the accumulation point of all $E_{c,n}$ from $\gamma<2$ with finite
$n=0,1..$. Every other point on the curve $E_{c,\xi}$ comes from the limit
$n\to\infty$, and different $\xi$ correspond to different ways how the double
limit $n\to\infty$ and $\gamma\to 2$ is taken. In the limit $\xi\to\infty$,
$E_{c}$ is the condensation energy for infinitesimally small gap function
$\Delta_{\infty}(\omega_{m})$.

In the next two Sections we present corroborative evidence for the special,
critical behavior of the $\gamma$-model with $\gamma=2$ from the analysis of
the gap function on the real frequency axis and in the upper half-plane of
frequency.

##  IV Gap equation along the real frequency axis

As we said in the Introduction, the analysis of the gap equation for the
$\gamma$ model along real frequency axis should generally be more revealing
than the analysis along the Matsubara axis, because the pairing interaction on
the real axis
$V(\Omega)=(\cos(\pi\gamma/2)+i\text{sign}{(\Omega)}\sin(\pi\gamma/2))({\bar{g}}/|\Omega|)^{\gamma}$
is complex. The real part of the interaction becomes repulsive for $\gamma>1$,
and the imaginary part vanishes at $\gamma=2$ for any non-zero $\Omega$. This
makes the $\gamma=2$ case special.

We present the results for $\Delta(\omega)$ on the real axis in the same order
as in previous section: we first obtain the solution of the linearized gap
equation, which we label $\Delta_{\infty}(\omega)$, then analyze the solution
$\Delta_{0}(\omega)$, and then show that there is a one-parameter continuous
set of solutions $\Delta_{\epsilon}(\omega)$ in between
$\Delta_{\infty}(\omega)$ and $\Delta_{0}(\omega)$.

###  IV.1 Linearized gap equation in real frequencies

The linearized gap equation in real frequencies is obtained by taking the
limit $\Delta(\omega)\rightarrow 0$ in (9). We again introduce
$D_{\infty}(\omega)=\Delta_{\infty}(\omega)/\omega$ and re-write the gap
equation as

$D_{\infty}({\omega})=-\frac{{\bar{g}}^{2}}{\omega}\left[i\frac{\pi}{2}~{}\frac{dD_{\infty}(\omega)}{d\omega}\text{sign}\omega+\frac{D_{\infty}(\omega)}{\omega}+\int_{0}^{\infty}\frac{d{\omega}^{\prime}}{(|\omega|+\omega^{\prime})^{2}}~{}\Re
D_{\infty}(\omega^{\prime})\right],$ (41)

where $\Re$ stands for the real part. The $D_{\infty}({\omega})$ term in the
l.h.s. of (41) is the analog of $D_{\infty}(\omega_{m})$ in the l.h.s. of the
gap equation (11) on the Matsubara axis, and, like there, it originates from
the bare $\omega$ term in the fermionic Green’s function. Neglecting this
term, we find that the solution of (41) is

$D_{\infty}(\omega)=-2i\epsilon\cos{\left[\beta\left(\log{\left(\frac{\omega}{{\bar{g}}}\right)^{2}}-i\pi\text{sign}(\omega)\right)+\phi\right]},$
(42)

where $\beta=0.38187$ the same as in (13), and $\epsilon$ is infinitesimally
small. We note that this $D_{\infty}(\omega)$ can be obtained from
$D_{\infty}(\omega_{m})$, Eq. (13), by rotating from $i\omega_{m}$ to
$\omega+i0$. In explicit form,

$\displaystyle
D^{\prime}_{\infty}(\omega)=2\epsilon\sin{\left(\beta\log{\left(\frac{\omega}{{\bar{g}}}\right)^{2}}+\phi\right)}\sinh(\pi\beta)\text{sign}{\omega}$
$\displaystyle
D^{{}^{\prime\prime}}_{\infty}(\omega)=-2\epsilon\cos{\left(\beta\log{\left(\frac{\omega}{{\bar{g}}}\right)^{2}}+\phi\right)}\cosh(\pi\beta)$
(43)

Observe that $D^{\prime}_{\infty}(-\omega)=-D^{\prime}_{\infty}(\omega)$ and
$D^{{}^{\prime\prime}}_{\infty}(-\omega)=D^{{}^{\prime\prime}}_{\infty}(\omega)$,
as it should be. The relation

$\int_{0}^{\infty}dx~{}\frac{x^{2i\beta}}{(x+1)^{2}}=\frac{2\pi\beta}{\sinh(2\pi\beta)}=\frac{1}{\sinh^{2}(\pi\beta)},$
(44)

is useful for the verification that $D_{\infty}(\omega)$ satisfies Eq. (41)
without $D_{\infty}({\omega})$ in the l.h.s. Using another relation

$\int_{0}^{\infty}dx~{}\frac{x^{i\beta}}{x-1}=i\pi\coth(\pi\beta),$ (45)

one can verify that $D^{\prime}_{\infty}$ and $D^{{}^{\prime\prime}}_{\infty}$
satisfy KK relations:

$\frac{2}{\pi}\int
dx\frac{D^{\prime}_{\infty}(x)x}{x^{2}-\omega^{2}}=-D^{{}^{\prime\prime}}_{\infty}(\omega),\quad\frac{2\omega}{\pi}\int
dx\frac{D^{{}^{\prime\prime}}_{\infty}(x)}{x^{2}-\omega^{2}}=D^{{}^{\prime}}_{\infty}(\omega),$
(46)

where the integrals are principle values.

Figure 5: Real and imaginary part (a) and the phase $\eta_{\infty}(\omega)$
(b) of $D_{\infty}(\omega)$. The periodicity of oscillation is set by
$[(\omega/{\bar{g}})^{2}+\log(\omega/{\bar{g}})^{2}]/\pi$.

We next consider $\rvert\omega\rvert>{\bar{g}}$. To obtain
$D_{\infty}(\omega)$ in this region, we take as an input the exact solution on
the Matsubara axis and analytically continue it to the real axis. By
construction, this can be done by replacing $\omega_{m}$ by $(-iz)$ \- the
function $\Delta_{\infty}(z)$ is guaranteed to be analytic in the upper half-
plane of frequency. However, because we don’t have the exact analytical
expression for $D_{\infty}(\omega_{m})$ for an arbitrary $\omega_{m}$, we have
to replace $\omega_{m}$ by $-i(\omega+i0)$ in Eq. (15) and obtain
$\Delta_{\infty}(\omega)$ by integrating over $k$. For small
$\omega<{\bar{g}}$, we find, after this integration, series of corrections to
(42) in powers of $\omega/{\bar{g}}$. For large $\omega>{\bar{g}}$, the
largest contribution to $\Delta_{\infty}(\omega)$ comes from the continuation
of the universal oscillating term $\Delta_{\infty;u}$, Eq. (19). Upon rotation
to the real axis, this term splits into two. One remains exponentially small,
but in the other the exponential factor cancels out. As a result, on the real
axis we have (see Appendix B for details)

$D_{\infty;u}(\omega)\sim\sqrt{2}\epsilon
e^{{i\over\pi}\left[\left({\omega\over{\bar{g}}}\right)^{2}+\log\left({\omega\over{\bar{g}}}\right)^{2}\right]}.$
(47)

Other contributions contain powers of ${\bar{g}}/|\omega|$ and are smaller.
Neglecting them, we obtain $D_{\infty}(\omega)=D_{\infty;u}(\omega)$ at
$\omega\gg{\bar{g}}$.

Comparing this form with (42), we see that both $D^{\prime}_{\infty}(\omega)$
and $D^{{}^{\prime\prime}}_{\infty}(\omega)$ continue oscillating at
$\omega>{\bar{g}}$, but with the period set predominantly by
$(\omega/{\bar{g}})^{2}$ rather than by $\log{(\omega/{\bar{g}})^{2}}$. In
Fig. 5, we plot real and imaginary parts of $D_{\infty}(\omega)$ and the phase
of the gap, $\eta_{\infty}(\omega)$, defined via
$D_{\infty}(\omega)=|D_{\infty}(\omega)|e^{i\eta_{\infty}(\omega)}$, or,
equivalently, via
$\eta_{\infty}(\omega)={\text{Im}}\log{D_{\infty}(\omega)}={\text{Im}}\log{\Delta_{\infty}(\omega)}$.
We see that the phase winds up an infinite number of times between $\omega=0$
and $O(\bar{g})$ and between $O(\bar{g})$ and $\infty$. Oscillations at
$\omega<{\bar{g}}$ are directly related to oscillations of
$\Delta_{\infty}(\omega_{m})$ on the Matsubara axis, and there is one-to-one
correspondence between each phase winding by $2\pi$ on a real axis and a
vortex on the Matsubara axis. Oscillations and phase winding at
$\omega>{\bar{g}}$ are present on the real axis, but not on the Matsubara
axis. It is natural to relate this discrepancy to the fact that the pairing
interaction is attractive on the Matsubara axis, but on the real axis,
$V^{\prime}(\Omega)$ is repulsive, and a non-zero $D_{\infty}(\omega)$ comes
from $V^{{}^{\prime\prime}}(\Omega)\propto\delta(\Omega^{2})$ (see more on
this below).

###  IV.2 The function $D_{0}(\omega)$

We now consider the opposite limit of the real-axis partner of sign-preserving
$D_{0}(\omega_{m})$. At $\omega_{m}<{\bar{g}}$,
$D_{0}(\omega_{m})\approx\Delta_{0}(0)/\omega_{m}$, and $D_{0}(\omega)$ on the
real axis must also be close to $D_{0}(0)/\omega$ at $\omega<{\bar{g}}$. At
larger $\omega_{m}>{\bar{g}}$, we will see that $D_{0}(\omega_{m})$ and
$D_{0}(\omega)$ are very different: $D_{0}(\omega_{m})$ decays as
$1/\omega^{3}_{m}$, while $D_{0}(\omega)$ does not decay and oscillates in
sign.

The solution of the gap equation along the real axis for $\omega>{\bar{g}}$
has been found by Combescot [9], who build his analysis on earlier results by
Karakozov, Maksimov, and Mikhailovsky [7] and by Marsiglio and Carbotte [8].
We follow Ref. [9] below.

It is convenient to introduce $\phi_{0}(\omega)$ via
$D_{0}(\omega)=1/\sin{\phi_{0}(\omega)}$ and re-express the gap equation (9)
at $T=0$ as the equation on $\phi(\omega)$. The equation is

$\frac{d\phi_{0}(\omega)}{d\omega}=\frac{2}{\pi{\bar{g}}^{2}}\left[\omega
B(\omega)-A(\omega)\sin{\phi_{0}}\right]$ (48)

where $A(\omega)$ and $B(\omega)$ are given by Eq. (8). The initial condition
for $\phi_{0}$ is $\phi_{0}({\bar{g}})\approx{\bar{g}}/\Delta_{0}(0)=O(1)$,
consistent with $\phi_{0}(\omega)\approx\omega/\Delta_{0}(0)$ at
$\omega<{\bar{g}}$.

At $\omega\geq{\bar{g}}$, $B(\omega)\approx 1+{\bar{g}}^{2}/\omega^{2}$ and
$A(\omega)\approx-\alpha{\bar{g}}^{3}/\omega^{2}$, where $\alpha\approx 1.27$
(Ref. [9]) . The $A(\omega)$ term can then be neglected if $\phi_{0}(\omega)$
is real, as we will assume and then verify. Without this term, Eq. (48) can be
solved easily, and the result is

$\phi_{0}(\omega)\approx\frac{1}{\pi}\left(\log{\left(\frac{\omega}{{\bar{g}}}\right)^{2}}+\left(\frac{\omega}{{\bar{g}}}\right)^{2}+C\right)$
(49)

where $C={\bar{g}}/\Delta_{0}(0)-1/\pi$. We see that $\phi_{0}(\omega)$ is
real, as we anticipated. We note that this $\phi_{0}(\omega)$ coincides with
the argument of the exponent for $D_{\infty}(\omega)$ in (47)

The function

$D_{0}(\omega)=\frac{1}{\sin{\phi_{0}(\omega+i0)}}$ (50)

is a sign-changing function of $\omega$, whose real part almost diverges at a
set of frequencies where $\phi_{0}(\omega)=p\pi$, and $p=1,2..$ is an integer.
The imaginary component $D^{{}^{\prime\prime}}_{0}(\omega)$ is a set of
$\delta-$functions at these frequencies. We plot the gap function
$\Delta_{0}(\omega)=\omega D_{0}(\omega)$ in Fig.6.

Figure 6: $\Delta_{0}(\omega)=\omega/\sin{\phi_{0}(\omega)}$ for
$\phi_{0}(\omega)$ given by (49). The real part of the gap
$\Delta^{\prime}_{0}(\omega)$ diverges at the set of points where
$\phi_{0}(\omega)=p\pi$, $p=1,2...$. The imaginary part
$\Delta^{{}^{\prime\prime}}_{0}(\omega)$ is a set of $\delta-$functions at
these points. The behavior of $\Delta^{\prime}_{0}(\omega)$ has been obtained
in Refs. [7, 8, 9, 11].

To analyze the phase winding, we again introduce the phase factor via
$D_{0}(\omega)=|D_{0}(\omega)|e^{i\eta_{0}(\omega)}$ and consider how
$\eta_{0}(\omega)$ varies at $\omega\geq{\bar{g}}$. The imaginary component
$D^{{}^{\prime\prime}}_{0}(\omega)$ in (50) is infinitesimally small, except
in the vicinity of $\omega_{p}$, where $\phi_{0}(\omega_{p})=p\pi$. We use Eq.
(49) for $\phi_{0}(\omega)$ and express $D_{0}(\omega)$ near each such point
as

$D_{0}(\omega)\approx\frac{\pi{\bar{g}}^{2}}{2}\frac{(-1)^{p}\omega_{p}}{\omega^{2}_{p}+{\bar{g}}^{2}}\frac{1}{\omega-\omega_{p}+i\delta}.$
(51)

Then

$e^{i\eta_{0}(\omega)}=(-1)^{p}\frac{\omega-\omega_{p}-i\delta}{\sqrt{(\omega-\omega_{p})^{2}+\delta^{2}}}$
(52)

Figure 7: Variation of the phase of the gap $\eta_{0}(\omega)$
($\Delta_{0}(\omega)=|\Delta_{0}(\omega)|e^{i\eta_{0}(\omega)}$). We restrict
$\eta_{0}(\omega)$ to $(-\pi,\pi)$. Phase slips of $\eta_{0}(\omega)$ continue
up to infinite frequency.

We plot $\eta_{0}(\omega)$ in Fig.7. We see that the phase rapidly changes by
$\pi$ around each $\omega_{p}$. If we restrict $\eta_{0}(\omega)$ to
$(-\pi,\pi)$, we find that the phase jumps by $2\pi$ in between $\omega_{p}$
and $\omega_{p+1}$. The number of $\omega_{p}$ is infinite, hence the number
$2\pi$ jumps is also infinite. We reiterate that behavior has no analog the
Matsubara axis, where $D_{0}(i\omega_{m})$ is real and positive for
$\omega_{m}>0$, hence $\eta_{0}=0$. 444For a generic $\gamma$,
$\Delta_{0}(\omega_{m})$ and $\Delta^{{}^{\prime\prime}}_{0}(\omega)$ are
related by Cauchy formula:
$\Delta_{0}(\omega_{m})=(2/\pi)\int_{0}^{\infty}d\omega\Delta^{{}^{\prime\prime}}_{0}(\omega)\omega/(\omega^{2}+\omega^{2}_{m})$.
For $\gamma<2$, typical $\omega$ are of order $\omega_{m}$, and to reproduce
$\Delta_{0}(\omega_{m})\propto 1/|\omega_{m}|^{\gamma}$ one needs
$\Delta^{{}^{\prime\prime}}_{0}(\omega)=\sin(\pi\gamma/2)\text{
sign}\omega/|\omega|^{\gamma}$. The case $\gamma=2$ is an exception here
because $1/\omega^{2}_{m}$ dependence of $\Delta_{0}(\omega_{m})$ is obtained
by pulling $1/\omega^{2}_{m}$ out of the denominator in the Cauchy formula.
The remaining integral is determined by $\omega=O({\bar{g}})$ rather than
$O(\omega_{m})$. Because of this, the fact that $\Delta_{0}(\omega_{m})\propto
1/\omega^{2}_{m}$ at large frequencies does not imply that
$\Delta^{{}^{\prime\prime}}(\omega)$ must behave as $1/\omega^{2}$.

###  IV.3 The one-parameter set of gap functions

We follow the same strategy as in the analysis on the Matsubara axis and
expand the non-linear gap equation (9) in powers of $D^{2}$. We search for the
solution in the form

$D(\omega)=\sum^{\infty}_{j=0}D^{(2j+1)}(\omega)$ (53)

where $D^{(1)}(\omega)=D_{\infty}(\omega)$ and higher-order terms are obtained
by solving Eq. (9) iteratively. For $\omega<{\bar{g}}$, we use Eq. (42) for
$D_{\infty}(\omega)$. The computational steps are the same as in Sec. III.3,
and we obtain

$D_{\epsilon}(\omega)=-2i\epsilon\left(\cos{{\tilde{f}}_{\epsilon}(\omega)}+\epsilon^{2}Q_{3,\epsilon}\cos{3{\tilde{f}}_{\epsilon}(\omega)}+\epsilon^{4}Q_{5,\epsilon}\cos{5{\tilde{f}}_{\epsilon}(\omega)}+...\right)$
(54)

where $Q_{i,\epsilon}$ are the same as in (37) and

${\tilde{f}}_{\epsilon}(\omega)=\beta_{\epsilon}\left(\log{\left(\frac{\omega}{{\bar{g}}}\right)^{2}}-i\pi\text{sign}{\omega}\right)+\phi_{\epsilon}$
(55)

This $D_{\epsilon}(\omega)$ could also be obtained directly from (37) by
replacing $\log{\omega^{2}_{m}}$ by $\log{\omega^{2}}-i\pi$ in each term in
(37).

We recall that the continuous set exists for $\epsilon\leq\epsilon_{cr}$. For
any $\epsilon<\epsilon_{cr}$, $D(\omega)$ oscillates an infinite number of
times down to $\omega=0$. As $\epsilon$ approaches $\epsilon_{cr}$, log-
oscillations shift to progressively smaller frequencies. At
$\epsilon=\epsilon_{cr}$, $\beta_{\epsilon}$ vanishes and log-oscillations
disappear. The behavior of $D(\omega)$ at $\omega\to 0$ at
$\epsilon\to\epsilon_{cr}$ depends on how the double limit $\omega\to 0$ and
$\epsilon\to\epsilon_{cr}$ is taken.

Like we did in Sec. III.3, we introduce
$\xi=(\epsilon_{cr}-\epsilon)/\epsilon$ and re-express
$\Delta_{\epsilon}(\omega)$ as $\Delta_{\xi}(\omega)$. The two limits
$\epsilon=0$ and $\epsilon=\epsilon_{cr}$ now correspond to $\xi=\infty$ and
$\xi=0$, respectively. This brings the notations in line with the ones we used
in Secs. IV.1 and IV.2.

On the Matsubara frequency, all $\Delta_{\xi}(\omega_{m})$ behave in the same
way at $\omega_{m}>{\bar{g}}$: $\Delta_{\xi}\propto 1/\omega^{2}_{m}$. On the
real axis, the dependence on $\xi$ is more complex. To see this, we use the
solution of the linearized gap equation $D^{(1)}\propto
ie^{i\phi_{0}(\omega)}$ with $\phi_{0}(\omega)$, given by (49), and evaluate
$D^{(2n+1)}$ in order-by-order expansion of the non-liner gap equation in
$D^{2}$. Collecting the series, we obtain the closed form expression

$D_{\xi}(\omega)=\frac{-2ie^{i\phi_{0}(\omega)}}{1+\xi-e^{2i\phi_{0}(\omega)}/(1+\xi)}=\frac{1}{\sin[\phi_{0}(\omega)+i\log(1+\xi)]}$
(56)

This expression can be equivalently obtained by solving Eq. (48) for
$\phi(\omega)$ with the initial condition
$\phi({\bar{g}})={\bar{g}}/\Delta_{0}(0)+i\log{(1+\xi)}$.

The parameter $\xi$ runs between $0$ and $\infty$. For $\xi=0$, Eq. (56)
yields $D_{0}(\omega)=1/\sin{\phi_{0}(\omega)}$, which agrees with (50) (one
should add $i0$ to $\omega$ in this case). At $\xi\to\infty$ we recover, by
construction, the solution of the linearized gap equation. For any $\xi$,
including $\xi=0$, $D_{\xi}(\omega)$ oscillates up to an infinite frequency,
and its phase $\eta_{\xi}(\omega)$ winds up by an infinite number of $2\pi$
between $\omega\sim{\bar{g}}$ and $\omega=\infty$.

We see therefore that in both limits $\omega\ll{\bar{g}}$ and
$\omega\gg{\bar{g}}$, the solutions of the non-linear gap equation form a
continuous one-parameter set, Eqs. (54) and (56). We conjecture that for any
value of $\xi$, one can use a free phase factor $\phi_{\xi}$ in (55) to merge
small-$\omega$ and large-$\omega$ expressions into a single $D_{\xi}(\omega)$.

###  IV.4 density of states

The fermionic DoS is defined as
$N(\omega)=(-N_{0}/\pi){\text{Im}}G_{l}(\omega)$, where $N_{0}$ is the DoS in
the normal state and

$G_{l}(\omega)=-i\pi\sqrt{\frac{1}{1-D^{2}(\omega)}}$ (57)

is a retarded Green’s function, integrated over the dispersion.

In a BCS superconductor,
$N(\omega)\propto\text{Re}\omega/\sqrt{\omega^{2}-\Delta^{2}}$ vanishes at
$\omega<\Delta$, has an integrable singularity at $\omega=\Delta+0$, and is
non-zero for all $\omega>\Delta$ because quasiparticle states in a BCS
superconductor form a continuum
$\omega=\pm\sqrt{\Delta^{2}+(\epsilon_{k}-\mu)^{2}}$. In our case, the form of
$N(\omega)=N_{\xi}(\omega)$ strongly depends on $\xi$. At small
$\omega<{\bar{g}}$, $N_{\xi}(\omega)$ remains finite down to $\omega=0$ for
all $\xi>0$. In this respect, all such solutions describe gapless
superconductivity. The gap function $\Delta_{0}(\omega)$ tends to a finite
$\Delta_{0}(0)$ at small $\omega$, and the corresponding $N_{0}(\omega)$
vanishes, like in BCS superconductor. We show this in Fig.8(a).

Figure 8: DoS $N_{\xi}(\omega)$ for (a) $\omega<\bar{g}$ and (b)
$\omega>\bar{g}$ and for different $\xi$. For all $\xi>0$, $N_{\xi}(\omega)$
remains finite down to $\omega=0$ (a gapless superconductivity). For $\xi=0$,
the DoS $N_{0}(\omega)$ vanishes at small $\omega$ and has a set of
$\delta-$functional peaks at $\omega>{\bar{g}}$. In panel (c) we present the
schematic plot of the DoS at all frequencies.

At $\omega>{\bar{g}}$, $\Delta_{\xi}(\omega)$ is given by (56), and
$N_{\xi}(\omega)=N_{0}\Im\tan[\phi_{0}(\omega)+i\log{(1+\xi)}]$. For $\xi>0$,
$N_{\xi}(\omega)$ oscillates around $N_{0}$ up to $\omega=\infty$. The
amplitude of the oscillations increases with decreasing $\xi$. For $\xi=0$,
$N_{0}(\omega)=N_{0}\delta/(\cos^{2}\phi_{0}(\omega)+\delta^{2})$, where
$\delta=0+$. This DoS consists of a set of $\delta-$functions at the points
$\omega_{k}$, for which $\phi_{0}(\omega_{k})=\pi/2+k\pi$ ($k$ is an integer)
(Refs [9, 11, 12]). We show this in Fig.8(b). In Fig.8(c) we show
$N_{\xi}(\omega)$ in the whole range of frequencies.

The function $N_{0}(\omega)$ is the true DoS at $T=0$, as the $\xi=0$ solution
has the lowest condensation energy. It is different from the DoS in a BCS-type
superconductor, which is non-zero at all $\omega>\Delta$ and approaches
$N_{0}$ at $\omega\to\infty$. We emphasize that a qualitative distinction
holds only for $\gamma=2$. For smaller $\gamma$, the DoS for the $n=0$
solution evolves as a function of frequency, but still remains non-zero at all
$\omega>\Delta$ and approaches $N_{0}$ at infinite $\omega$ (see Paper IV).

In a BCS superconductor, a continuous $N(\omega)$ at $\omega>\Delta$ is the
consequence of the fact that fermionic energy $E_{k}$ is a continuous function
of the normal state dispersion $\epsilon_{k}$,
$E_{k}=\sqrt{\epsilon^{2}_{k}+\Delta^{2}}$. The form of $N_{0}(\omega)$ as a
set of $\delta$-functional peaks raises the issue whether fermionic energies
get quantized at $\gamma=2$. To address this issue, we compute the total
weight of each level: $N_{k}=(1/2\pi)\int N_{0}(\omega)/N_{0}$, where the
integration is confined to the vicinity of $\omega_{k}$. Using
$\phi_{0}(\omega)\approx\omega^{2}/\pi$, we obtain $N_{k}=1/\sqrt{8(1+2k)}$.
We see that $N_{k}<1$ for all $k$. Because of this, $\omega_{k}$ cannot be
viewed as true quantized fermionic energy levels, as a fermion is necessary
distributed between $\omega_{k}$ with different $k$.

##  V Gap function in the upper frequency half-plane.

Comparing $D_{\xi}(\omega)$ and $D_{\xi}(\omega_{m})$, we see that they are
similar at small frequencies, but very different at
$\omega,\omega_{m}>{\bar{g}}$. Indeed, on the real axis, the phase
$\eta_{\xi}(\omega)$ winds up by an infinite number of $2\pi$ between
$\omega=O({\bar{g}})$ and $\omega=\infty$, while near the Matsubara axis,
$\eta_{\xi}(\omega_{m})=0$ in this frequency range. The discrepancy implies
that phase winding must end somewhere between the real and the Matsubara axis.
We now argue that there is a set of vortices in the upper frequency half-
plane, at $|z|\geq{\bar{g}}$ and the phase winding drops by $2\pi$ each time
the axes of $z$ passes through a vortex upon rotation away from the real axis.

We use the Cauchy relation

$\Delta(z)=\frac{2}{\pi}\int_{0}^{\infty}dx\frac{x\Delta^{{}^{\prime\prime}}(x)}{x^{2}-z^{2}}$
(58)

to extend the gap function $\Delta(x)=xD(x)$ from the real axis to complex
$z=\omega^{\prime}+i\omega^{{}^{\prime\prime}}$ with
$\omega^{{}^{\prime\prime}}>0$. We use Eq. (56) for the gap function as we
expect vortices to be present at $|z|>{\bar{g}}$. Like before, we first
consider the cases $\xi=0$ and $\xi\to\infty$, and then extend the analysis to
arbitrary $\xi$.

###  V.1 case $\xi=0$

Using the expansion near $\phi(\omega_{p})=p\pi$, Eq. (51), we approximate
$\Delta^{{}^{\prime\prime}}_{0}(\omega)$ as

$\Delta^{{}^{\prime\prime}}_{0}(\omega)\approx\frac{\pi^{2}{\bar{g}}^{2}}{2}\sum_{p=1}^{\infty}\frac{(-1)^{p+1}\omega^{2}_{p}}{\omega^{2}_{p}+{\bar{g}}^{2}}\delta(\omega-\omega_{p}).$
(59)

Substituting into (58), we obtain

$\Delta_{0}(z)=\pi{\bar{g}}^{2}\sum_{p=1}^{\infty}\frac{(-1)^{p+1}\omega^{3}_{p}}{\omega^{2}_{p}+{\bar{g}}^{2}}\frac{1}{\omega^{2}_{p}-z^{2}}.$
(60)

Here $\omega_{p}$ is a solution of $\phi_{0}(\omega_{p})=\pi p$, where
$\phi_{0}(\omega)$ is given by (49). We verified numerically that KK relations
on the real axis are satisfied, i.e., if we use (59), we reproduce with high
accuracy $\Delta^{\prime}(\omega)$. On the Matsubara axis, $z=i\omega_{m}$,
Eq. (60) yields, at $\omega_{m}\gg{\bar{g}}$,

$\Delta_{0}(\omega_{m})=a\frac{{\bar{g}}^{3}}{\omega^{2}_{m}}$ (61)

where
$a=\pi\sum_{p=1}^{\infty}\left[(-1)^{p+1}\omega^{3}_{p}/((\omega^{2}_{p}+{\bar{g}}^{2}){\bar{g}})\right]$.
Approximating $\omega_{p}$ by ${\bar{g}}\pi\sqrt{p}$, we find $a=2.56$. The
number is somewhat larger than $1.27$, obtained by solving the gap equation on
the Matsubara axis (Ref. [9] and Sec. III.2). The difference likely comes from
subleading terms in $\phi(\omega)$.

We plot $\Delta_{0}(z)$ for a generic $z$ in the upper half-plane in Fig.9

Figure 9: $\Delta_{0}(z)$ in the upper half plane. Upper panel:
$\log{|\Delta_{0}(z)|}$. Blue spots mark the locations of dynamical vortices,
where $|\Delta_{0}(z)|=0$. Lower panel: the phase of the gap $\eta_{0}(z)$
defined via $\Delta_{0}(z)=|\Delta_{0}(z)|e^{i\eta_{0}(z)}$. The phase slips
by $2\pi$ upon crossing a white line in the direction from near-white to dark-
blue. The white lines are the locations of points where
$\Delta^{{}^{\prime\prime}}_{0}(z)=0$ and $\Delta^{\prime}_{0}(z)<0$.

We clearly see that there is a set of points, where
$\Delta^{\prime}_{0}(z)=\Delta^{{}^{\prime\prime}}_{0}(z)=0$. These points are
the centra of dynamical vortices with anti-clockwise circulation $2\pi$. The
vortices are located along a particular line in the complex plane. The set
extends to an infinite frequency, i.e., the number of vortices is infinite.
This is consistent with an infinite phase winding along the real axis. We
verified that if we use a more accurate form of $\omega_{p}$, the positions of
the vortices shift a bit, but their number remains infinite.

To see how the winding number changes once we rotate from the real to the
Matsubara axis, we introduce $z=|z|e^{i\psi}$ ($\psi=0$ along the positive
real semi-axis and $\pi/2$ along the Matsubara axis) and check the winding of
the phase of $\Delta_{0}(z)=|\Delta_{0}(z)|e^{i\eta_{0}(z)}$ between
$|z|\sim{\bar{g}}$ and $|z|\to\infty$ along the directions in the upper
frequency half-plane, specified by $\psi$. We show the results in Fig.10.

Figure 10: Phase variation $\eta_{0}(|z|,\psi)$ along different paths
specified by $\psi$, defined via $z=|z|e^{i\psi}$. Along real axis, $\psi=0$;
along Matsubara axis, $\psi=\pi/2$. Along the real axis the phase
$\eta_{0}(\omega)$ winds up an infinite number of times, i.e., the winding
number (the number of $2\pi$ phase slips) is infinite. For a finite $\psi$,
phase winding ends at some finite $|z|$, and the winding number becomes
finite.

We see that for any $\psi>0$, the phase $\eta_{0}(z)$ winds for $|z|$ below a
certain value, and then saturates. At larger $|z|$, both
$\Delta^{\prime}_{0}(z)$ and $\Delta^{{}^{\prime\prime}}_{0}(z)$ scale as
$1/|z|^{2}$ with no oscillations. Counting the total phase winding
$\delta\eta_{0}$ between $|z|=O({\bar{g}})$ and $|z|=\infty$, we see that
$\delta\eta_{0}=2\pi s$, where $s$ is an integer. It decreases by one every
time the direction set by $\psi$ passes through a vortex. The winding vanishes
at some $\psi\leq\pi/4$.

###  V.2 case $\xi=\infty$

We next consider the opposite limit $\xi=\infty$. The form of
$\Delta_{\infty}(z)$ can be obtained starting from (15) and replacing
$\omega^{2}_{m}$ by $|z|^{2}e^{i(2\psi-\pi)}$. This gives

$\Delta_{\infty}(z)\propto\int_{0}^{\infty}dk\left(b_{k}e^{-ik\log{|z|^{2}/{\bar{g}}^{2}}+(2\psi-\pi)k}+b_{-k}e^{ik\log{|z|^{2}/{\bar{g}}^{2}}-(2\psi-\pi)k}\right),$
(62)

where $b_{k}$ is defined in (16). We obtain $\Delta_{\infty}(z)$ by numerical
integration. We plot its phase $\eta_{\infty}(z)$ in Fig. 11. We again see
that there is an infinite array of vortices. The array extends to an infinite
frequency, where it approaches the real axis. The vortex arrangement in Fig.
11 is remarkable similar to that

Figure 11: The phase $\eta_{\infty}(z)$ of $\Delta_{\infty}(z)$ in the first
quarter of the complex plane of frequency ($\omega^{\prime}>0$,
$\omega^{{}^{\prime\prime}}>0$).

in Fig. 9 for $\xi=0$. Moreover, if we approximate $\phi(\omega)$ by the
leading term $(\omega/{\bar{g}})^{2}/\pi$, we find that the positions of the
vortices are at the same $z_{i}$ in both cases. We can see this by comparing
Fig.12(a) where $\xi\to\infty$ and Fig.12(c) where $\xi\to 0$. The gap
function $\Delta_{\xi}(z)$ are very similar in these two cases, despite that
the overall factors are different. The vortex positions for these two cases
are almost identical, as can be seen from Fig.12(d).

Figure 12: (a)-(c)Gap functions $\Delta_{\xi}$ defined in (56) for different
$\xi$ in the frequency upper half plane. Here we take
$\phi(\omega)\approx(\omega/\bar{g})^{2}/\pi$ (d) Comparison of the vortex
positions for different $\xi$ obtained by approximating $\phi(\omega)$ by
$(\omega/\bar{g})^{2}/\pi$. The results suggest that the positions of the
vortices almost do not depend on the value of $\xi$.

###  V.3 arbitrary $\xi$

The same infinite array of vortices exists for all $0<\xi<\infty$. As an
example, in Fig.12(b) we show the gap function for $\xi=1$. We clearly see
that there is an infinite array of vortices, similar to the ones for
$\xi=\infty$ and $\xi=0$, and the positions of vortices are almost
indistinguishable, see Fig.12(d). Analytically, if we use Eqs. (56) and (58),
we find that the positions of vortices are independent on $\xi$.

###  V.4 Essential singularity

There is another consequence of the existence of an infinite array of vortices
– each the gap function $\Delta_{\xi}(z)$ has an essential singularity at
$|z|=\infty$. Indeed, one can reach $|z|=\infty$ from the set of vortex
points, where $\Delta_{\xi}(z)=0$, or from the real axis, where
$\Delta_{\xi}(\omega)$ oscillates, and the amplitude of the oscillations does
not vanish at $\omega\to\infty$, hence neither
$\lim_{|z|\to\infty}\Delta_{\xi}(|z|)$ nor
$\lim_{|z|\to\infty}1/\Delta_{\xi}(|z|)$ exist. We emphasize that an essential
singularity is only present for $\gamma=2$. For smaller $\gamma$, phase
winding and associated vortices exist only at $|z|$ smaller than a certain,
$\gamma-$dependent value. At larger $|z|$, $\Delta(z)$ scales as
$1/|z|^{\gamma}$ and vanishes at $|z|=\infty$ no matter how this limit is
reached.

Further, for $\gamma=2$, the very existence of a non-zero $\Delta_{\xi}(z)$
for a generic $z$ away from vortex points, is ultimately related to an
essential singularity at $|z|=\infty$. The argument is that the set of vortex
points is complete, hence one can analytically continue the gap function from
this set to the upper half-plane of frequency in the same way as $\Delta(z)$
is obtained from a discrete set of Matsubara points $\omega_{m}=\pi T(2m+1)$
in standard diagrammatic calculations for interacting fermions. If this
analytical continuation was unique, we would obtain $\Delta(z)\equiv 0$,
because $\Delta(z)=0$ at the vortex points. For a non-zero $\Delta_{\xi}(z)$,
the extension must be multi-valued. A rigorous mathematical argument is that
this is the case when the end point of the set, $|z|=\infty$, goes outside the
domain of analyticity. This is exactly what we have because of an essential
singularity at $|z|=\infty$.

We conjecture that the multi-value nature of the extension is the reason why
the set of $\Delta(\omega)$ is a continuous one at $\gamma=2$. This is
plausible, particularly if the vortices are at the same $z_{i}$ for all $\xi$,
as Fig.12 seems to indicate. However, at the moment, we cannot prove this.

##  VI Finite $\omega_{D}$

###  VI.1 Gap equation at a finite $\omega_{D}$

We now consider the case when the bosonic mass is small but finite. By analogy
with the phonon case we call this mass $\omega_{D}$. On the Matsubara axis,
$\Delta_{0}(\omega_{m})$ changes little compared to the case $\omega_{D}=0$.
The set of $\Delta_{n}(\omega_{m})$ still exists at small $\omega_{D}$, but
becomes discrete and holds up to a finite $n_{max}$. In particular, there is
no solution of the linearized gap equation at $T=0$ for any non-zero
$\omega_{D}$. The value of $n_{max}$ can be estimated by noticing that if we,
e.g., depart from the solution on the Matsubara axis at $\omega_{D}=0$ and
compute corrections due to finite $\omega_{D}$, these corrections increase at
small $\omega_{m}$ and become $O(1)$ at $\omega_{m}\sim\omega_{D}$. A simple
experimentation shows that this sets $n_{max}$ at

$n_{max}\sim\frac{{\bar{g}}}{\omega_{D}}$ (63)

On the real axis, the gap equation still has the form $D(\omega)\omega
B(\omega)=A(\omega)+C(\omega)$, and $A(\omega)$ and $B(\omega)$ remain the
same as in (8), up to irrelevant small corrections. However, $C(\omega)$
changes to

$C(\omega)=-i\frac{\pi{\bar{g}}^{2}}{2\omega_{D}}\frac{D(\omega-\omega_{D})-D(\omega)}{\sqrt{1-D^{2}(\omega-\omega_{D})}}~{}\text{sign}\omega$
(64)

Expanding to first order in $\omega_{D}$ and introducing, as before,
$D(\omega)=1/\sin{\phi(\omega)}$, we obtain after straightforward algebra that
the gap equation reduces to

$\displaystyle{\dot{\phi}}-\frac{\omega_{D}}{2}\left((\dot{\phi})^{2}\tan{\phi(\omega)}+\ddot{\phi}\right)=\frac{2}{\pi{\bar{g}}^{2}}\left[\omega
B(\omega)-A(\omega)\sin{\phi(\omega)}\right]+...$ (65)

where dots stand for the terms with higher powers of $\omega_{D}$. A similar
equation at a finite $T$ instead of finite $\omega_{D}$ has been obtained by
Combescot [9]

For definiteness, let’s consider the case $\xi=0$. At $\omega\geq{\bar{g}}$,
$B(\omega)$ and $A(\omega)$ from (8) can be approximated by $B(\omega)\approx
1+{\bar{g}}^{2}/\omega^{2}$ and
$A(\omega)\approx-1.27{\bar{g}}^{3}/\omega^{2}$. To understand the effect of
$\omega_{D}$ we use as an input the solution at $\omega_{D}$ =0,
$\phi(\omega)\approx\omega^{2}/(\pi{\bar{g}}^{2})+i\delta$. Substituting this
input into (65), expanding near $\omega=\pi{\bar{g}}/\sqrt{2}$, where
$\phi(\omega)=\pi/2$, expressing
$\phi(\omega)=\phi^{\prime}(\omega)+i\phi^{\prime\prime}(\omega)$, and solving
for $\phi^{{}^{\prime\prime}}(\omega)$, we find that it jumps to
$O(\omega_{D})$ once $\omega$ exceeds $\pi{\bar{g}}/\sqrt{2}$. The same
happens at all $\omega_{n}=\pi{\bar{g}}/\sqrt{2}(2n+1)^{1/2}$, where
$\tan{\phi^{\prime}(\omega)}=0$. After $n$ jumps,
$\phi^{{}^{\prime\prime}}(\omega)$ becomes

$\phi^{{}^{\prime\prime}}(\omega)=\frac{\pi\omega_{D}}{\sqrt{2}\bar{g}}\sum_{m=1}^{n}\sqrt{2m+1}\approx\frac{2\pi\omega_{D}}{3{\bar{g}}}n^{3/2}\approx\frac{2\omega_{D}}{3\pi^{2}}\frac{\omega^{3}}{{\bar{g}}^{4}}.$
(66)

A more accurate, non-perturbative analysis of (65) shows that
$\phi^{{}^{\prime\prime}}(\omega)$ appears slightly before
$\phi^{\prime}(\omega)$ reaches $\pi/2$. This smoothes up the jumps, but the
functional form of $\phi^{{}^{\prime\prime}}(\omega)$ in (66) remains intact.
When both $\phi^{\prime}$ and $\phi^{{}^{\prime\prime}}$ are nonzero,
$D(\omega)$ is a complex function of frequency:

$D^{\prime}_{0}(\omega)=\omega\frac{\sin{\phi^{\prime}(\omega)}\cosh{\phi^{{}^{\prime\prime}}(\omega)}}{\sin^{2}{\phi^{\prime}(\omega)}+\sinh^{2}{\phi^{{}^{\prime\prime}}(\omega)}}~{}{\text{sign}}\omega$
(67)

and

$D^{\prime\prime}_{0}(\omega)=-\omega\frac{\cos{\phi^{\prime}(\omega)}\sinh{\phi^{{}^{\prime\prime}}(\omega)}}{\sin^{2}{\phi^{\prime}(\omega)}+\sinh^{2}{\phi^{{}^{\prime\prime}}(\omega)}}$
(68)

At $\omega>3\pi^{2}{\bar{g}}^{4}/(2\omega_{D})$,
$\phi^{{}^{\prime\prime}}(\omega)$ becomes larger than one. At such
frequencies, both $D^{\prime}(\omega)$ and $D^{\prime\prime}(\omega)$
oscillate with progressively decreasing magnitudes, approximately as the real
and the imaginary parts of

$-2ie^{i\omega^{2}/(\pi\bar{g}^{2})}e^{-\frac{2\omega_{D}}{3\pi^{2}\bar{g}}(\omega/{\bar{g}})^{3}},$
(69)

and the phase $\eta_{0}(\omega)$ gradually winds up as $\omega$ increases. We
show this in Fig.13.

Figure 13: The gap function $\Delta(\omega)$ (a) and the variation of its
phase $\eta(\omega)$ at $\omega_{D}=0.02\bar{g}$. From (69).

This behavior holds as long as $|A(\omega)|\ll\omega$, i.e.,
$\omega<\omega_{max}$, where

$\omega_{max}\sim{\bar{g}}\left(\frac{\bar{g}}{\omega_{D}}\log{\frac{\bar{g}}{\omega_{D}}}\right)^{1/3}$
(70)

At even larger frequencies, the $A(\omega)$ term cannot be neglected, and the
forms of $\phi^{\prime}(\omega)$ and $\phi^{{}^{\prime\prime}}(\omega)$
change. We show the full numerical solution of Eq. (65) in Fig.14.

Figure 14: (a) Numerical solution of Eq. (65) for a complex $\phi(\omega)$ at
$\omega_{D}=0.02\bar{g}$. (b) and (c): Asymptotic forms of
$\phi^{\prime}(\omega)$ and $\phi^{{}^{\prime\prime}}(\omega)$. At
$\omega<\omega_{max}$, $\phi^{\prime}(\omega)$ increases as $\omega^{2}$,
while $\phi^{\prime\prime}(\omega)$ first displays a step-like behavior and
then increases as $\omega^{3}$, with $\omega_{D}$ in the prefactor. At
$\omega>\omega_{max}$, $\phi^{\prime}(\omega)$ saturates at$(2m-1/2)\pi$,
where $m=19$ for our chosen $\omega_{D}$, and $\phi^{\prime\prime}(\omega)$
increases logarithmically. The corrections to asymptotic values oscillate with
the period set by $\omega^{2}/\pi$. In the numerical solution, we neglected
$\ddot{\phi}$ term in (65) compared to $(\dot{\phi})^{2}$ and verified that
this is a valid approximation.

We see that at $\omega>\omega_{max}$, $\phi^{{}^{\prime\prime}}(\omega)$ keeps
increasing, while $\phi^{\prime}(\omega)$ saturates. A simple analysis shows
that Eq. (65) is satisfied, up to corrections of order $\omega_{D}$, if

$\phi^{{}^{\prime\prime}}(\omega)=3\log{\frac{\omega}{{\bar{g}}}}+0.45,~{}~{}\phi^{\prime}(\omega)=-\frac{\pi}{2}+2m\pi$
(71)

where $m$ is an integer. Substituting this complex $\phi(\omega)$ into
$\Delta_{0}(\omega)=\omega/\sin{\phi(\omega)}\approx-2i\omega
e^{i\phi(\omega)}$, we see that at $\omega>\omega_{max}$, the real part of the
gap function gradually decreases as

$\Delta^{\prime}_{0}(\omega)=\frac{1.27{\bar{g}}^{3}}{\omega^{2}}$ (72)

To obtain $\Delta^{{}^{\prime\prime}}_{0}(\omega)$ at these frequencies, we
need to keep the $\omega_{D}$ term in the l.h.s. of (65) and obtain the
correction to (71), which we label as $\tilde{\phi}$. Solving perturbatively
for ${\tilde{\phi}}(\omega)$ we obtain

${\tilde{\phi}}(\omega)=f\left(\frac{\omega}{\omega_{max}}\right)~{}e^{i\omega^{2}/(\pi{\bar{g}}^{2})}$
(73)

where $f(...)$ is a decreasing function of the argument. The $\omega^{2}$
oscillations of ${\tilde{\phi}}(\omega)$ are clearly visible in the numerical
results for $\phi^{\prime}$ and $\phi^{{}^{\prime\prime}}$ in Fig. 14.
Substituting (73) into (68), we obtain

$\Delta^{{}^{\prime\prime}}_{0}(\omega)\sim\frac{{\bar{g}}^{3}}{\omega^{2}}f\left(\frac{\omega}{\omega_{max}}\right)~{}\cos{\frac{\omega^{2}}{\pi{\bar{g}}^{2}}}$
(74)

One can verify that an integer $m$ in (71) determines the number of $2\pi$
variations of $\eta_{0}(\omega)$ on the real axis and, equivalently, the
number of vortices at complex $z_{i}$. The value of $m$ decreases one-by-one
as $\omega_{D}$ increases and $\omega_{max}$ decreases. That $m$ is finite
implies that there is no essential singularity at $|z|=\infty$. Indeed, at the
largest frequencies $\Delta(\omega)\propto 1/\omega^{2}$.

For completeness, we verified that higher-order terms in $\omega_{D}$, which
we neglected in the l.h.s. of (65), become important at frequencies
$\omega\sim{\bar{g}}^{2}/\omega_{D}$, which well exceed $\omega_{max}$ and are
therefore irrelevant to our purposes.

##  VII Dressed superfluid stiffness

In this section we analyze superfluid stiffness and thermal corrections to a
superconducting order parameter. As we discussed in the Introduction, we
consider the $\gamma=2$-model as the double limit $\omega_{D}\to 0$,
$E_{F}\to\infty$, such that Migdal-Eliashberg parameter
$\lambda_{E}={\bar{g}}^{2}N_{0}/\omega_{D}$ remains small ($N_{0}\sim 1/E_{F}$
is the DoS per unit volume in the normal state). Accordingly, in the analysis
below we keep $\omega_{D}$ small, but finite.

###  VII.1 Bare stiffness

A superfluid stiffness is the ratio of the excess energy $E_{\eta}$ due to
inhomogeneous variation of the phase of a superconducting order parameter
$\Delta(r)=\Delta e^{i\eta(r)}$ and $\int dr\left(\nabla\eta(r)\right)^{2}$:
$E_{\eta}=\rho_{s}\int dr\left(\nabla\eta(r)\right)^{2}$. In the momentum
space,

$E_{\eta}=\rho_{s}~{}\sum_{q}q^{2}\eta_{q}^{2}$ (75)

A way to compute $\rho_{s}$ is to choose $\eta_{q}=\delta_{q,q_{0}}$ and
extract $\rho_{s}$ as the prefactor for $q^{2}_{0}$ term in the particle-
particle bubble (the sum of GG and FF terms) at zero frequency and finite $q$
(see Refs. [31, 33, 42]

At $\omega_{D}/{\bar{g}}>1$ , the system is in a weak coupling limit, and
superfluid stiffness at $T=0$ is a fraction of the Fermi energy,
$\rho_{s}(T=0)=E_{F}/(4\pi)$ (Refs. [31, 33]). This stiffness is much larger
than $T_{c}$ [43, *Emery1995]. At $T>0$, $\rho_{s}(T)$ drops and vanishes at
$T_{c}$, but at weak coupling a drop of $\rho_{s}$ occurs only in the
immediate vicinity of $T_{c}$.

At small $\omega_{D}/{\bar{g}}$, strong mass renormalization
$m^{*}/m=1+{\bar{g}}^{2}/\omega^{2}_{D}$ changes the stiffness to

$\rho_{s}(T=0)\sim
E_{F}\frac{\omega_{D}\Delta(0)}{{\bar{g}}^{2}}\sim\frac{T_{p}}{\lambda_{E}}$
(76)

where $T_{p}\sim\Delta(0)$ is the onset temperature of the pairing. As long as
$\lambda_{E}\leq 1$, $\rho_{s}(T=0)>T_{p}$.

We now relate the stiffness to the strength of thermal phase fluctuations of
$\Delta(r)=\Delta e^{i\eta(r)}$. For this, consider the correlator

$\langle\eta(r)\eta(0)\rangle=\frac{\int D[\eta]\eta(r)\eta(0)e^{-\rho_{s}\int
dr\left(\nabla\eta(r)\right)^{2}/T}}{\int D[\eta]e^{-\rho_{s}\int
dr\left(\nabla\eta(r)\right)^{2}/T}}$ (77)

We assume that in equilibrium $\eta(r)=0$ and expand $\langle
e^{i\eta(r)}\rangle$ as $1-\langle\eta^{2}(r)\rangle/2$. Transforming (77) to
the momentum space, we obtain
$\langle\Delta(r)\rangle=\Delta(1-\langle\eta^{2}\rangle)$, where

$\langle\eta^{2}\rangle=\frac{1}{N}\sum_{q}\frac{\prod_{q^{\prime}}\int
d\eta_{q^{\prime}}\eta^{2}_{q}~{}e^{-\rho_{s}q^{2}\eta^{2}_{q^{\prime}}/T}}{\prod_{q^{\prime}}\int
d\eta_{q^{\prime}}~{}e^{-\rho_{s}q^{2}\eta^{2}_{q^{\prime}}/T}}$ (78)

where $N$ is the number of particles in the system. Evaluating the integrals,
we obtain the conventional result [45]

$\langle\eta^{2}\rangle=\frac{T}{\rho_{s}}\frac{1}{N}\sum_{q}\frac{1}{q^{2}}.$
(79)

We assume for simplicity that spatial dimension $D$ is larger than $2$, in
which case the sum converges. By order of magnitude we then have

$\langle\eta^{2}\rangle~{}\sim\frac{T}{\rho_{s}(T)}$ (80)

As long as $\rho_{s}(T)>T$, fluctuation corrections to the order parameter are
small. This does not hold in the immediate vicinity of the onset temperature
of the pairing, $T_{p}$, but as long as $\rho_{s}(0)\gg T_{c}$, the $T$ range,
where fluctuations are strong and destroy phase coherence, is quite narrow,
i.e., superconducting $T_{c}$ remains close to $T_{p}$. We see that this holds
even when $\omega_{D}$ is small and the reduction of $\rho_{s}$ by mass
renormalization is strong.

###  VII.2 Dressing of $\rho_{s}$ by soft longitudinal fluctuations

We now argue that in our case the expression for $\langle\eta^{2}\rangle$ is
different due to the presence of a continuum gapless spectrum of condensation
energy, $E_{c,\xi}$, where, we remind, $\xi$ runs between $0$ and $\infty$,
and the bottom of the spectrum is at $\xi=0$. We will need states near the
bottom of the continuum, at $\xi\ll 1$. For such states, we assume

$E_{c,\xi}=E_{c,0}+b_{1}N_{0}N\bar{g}^{2}\xi^{2},$ (81)

where $b_{1}=O(1)$ and $N$ is the total number of particles. We will also need
superfluid stiffness $\rho_{s,\xi}$ for the states near the bottom of the
continuum. Evaluating the particle-particle susceptibility for a generic
$\Delta_{\xi}(\omega_{m})$ and extracting the $q^{2}$ term we obtain

$\rho_{s,\xi}\sim E_{F}\frac{\omega_{D}}{{\bar{g}}^{2}}\int
d\omega_{m}\frac{D^{2}_{\xi}(\omega_{m})}{1+D^{2}_{\xi}(\omega_{m})}$ (82)

For $\xi=0$, the integral is determined by $\omega_{m}\sim{\bar{g}}$, where
$D_{0}(\omega_{m})\sim 1$. This yields $\rho_{s,0}\sim
E_{F}\omega_{D}/{\bar{g}}\sim T_{p}/\lambda_{E}$, as in (76). For states with
$\xi>0$, the magnitude of $\Delta_{\xi}(\omega_{m})$ is reduced, and the
stiffness gets smaller. We assume that for the states near the bottom of the
continuum, the stiffness is obtained by expanding to first order in $\xi$:

$\rho_{s,\xi}=\rho_{s,0}\left(1-b_{2}\xi\right)$ (83)

where $b_{2}=O(1)$ is positive. The extra energy of a given state $\xi$ due to
phase variation is

$E_{\eta,\xi}=\rho_{s,\xi}\sum_{q}q^{2}\eta_{q}^{2}$ (84)

We assume (see reasoning below) that all states near the bottom of a continuum
contribute to the variation of the phase, i.e., the averaging in
$\langle(\eta_{q})^{2}\rangle$ is over both $\eta_{q}$ and $\xi$ with the
weight factor $e^{-E_{\xi}/T}$, where

$E_{\xi}=E_{c,\xi}+E_{\eta,\xi}=E_{0}+\delta E_{\xi},$ (85)

and

$\displaystyle E_{0}=E_{c,0}+\rho_{s,0}\sum_{q}q^{2}\eta_{q}^{2}$
$\displaystyle\delta
E_{\epsilon}=b_{1}N_{0}g^{2}\xi^{2}-b_{2}\rho_{s,0}\xi\sum_{q}q^{2}\eta_{q}^{2}$
(86)

If we neglected $\delta E_{\xi}$, we would obtain the same result as before:

$\langle\eta^{2}\rangle=\frac{T}{\rho_{s,0}}\frac{1}{N}\sum_{q}\frac{1}{q^{2}}$
(87)

Keeping $\delta E_{\xi}$ we find that $\langle\eta^{2}\rangle$ has an
additional overall factor, which we label as $I_{T}$. Dropping for simplicity
numerical prefactors $b_{1}$ and $b_{2}$, we obtain after integrating over
$\eta_{q^{\prime}}$

$I_{T}=\frac{\int d\xi e^{-Nf(\xi)\frac{1}{1-\xi}}}{\int d{\xi}e^{-Nf(\xi)}}$
(88)

where

$f(\xi)=\frac{N_{0}{\bar{g}}^{2}}{T}{\xi}^{2}-\frac{\xi}{2}$ (89)

We assume that the measure of the integration over $\xi$ is non-singular. The
linear in $\xi$ term in $f(\xi)$ comes from integration over
$\eta_{q^{\prime}}$ with $q^{\prime}\neq q$ (see Eq. 78). Each integration
over $q^{\prime}$ yields $1/\sqrt{1-\xi}$, and the product of the integrals
over all $q^{\prime}$ yields $1/(1-\xi)^{N/2}=e^{-(N/2)\log{(1-\xi)}}\approx
e^{(N/2)\xi}$.

At small $T$, the function $f(\xi)$ in (89) has a minimum at
${\xi}=T/(4N_{0}{\bar{g}}^{2})\sim T/(4\omega_{D}\lambda_{E})$. Then
$I_{T}=1/(1-T/(4\omega_{D}\lambda_{E}))$ and

$\langle\eta^{2}\rangle\sim\frac{T\lambda_{E}}{T_{p}}\frac{1}{1-\frac{T}{4\omega_{D}\lambda_{E}}}$
(90)

We see that the renormalizations coming from the low-energy states of the
continuum spectrum of the condensation energy hold in powers of
$T/\omega_{D}$. With these renormalizations, the fully dressed stiffness is

$\rho_{s}(T)=\frac{T_{p}}{\lambda_{E}}\left(1-\frac{T}{4\omega_{D}\lambda_{E}}\right)$
(91)

We see from that the value of $\rho_{s}(T)$ at $T\to 0$ and $\omega_{D}\to 0$
depends on the order of limits. At $T=0$, $\rho_{s}(0)=T_{p}/\lambda_{E}$ is
finite and exceeds $T_{p}$. At $\omega_{D}\to 0$ the corrections to stiffness
rapidly increase with $T$, and $\rho_{s}(T)$ becomes comparable to $T$ at
$T\sim\omega_{D}\lambda_{E}$. For the largest $\lambda_{E}\sim 1$, at which
our theory is valid, this holds at $T\sim\omega_{D}$. It is tempting to
associate this temperature with the actual $T_{c}$ above which the system
looses long-range phase coherence.

Further, there is an analogy between finite $\omega_{D}$ and finite
$2-\gamma$, as the two have similar effect on the gap function (see Paper IV).
Replacing $\omega_{D}$ by ${\bar{g}}(2-\gamma)$, we find that at
$\omega_{D}=0$ and $\gamma<2$, superconducting $T_{c}\sim{\bar{g}}(2-\gamma)$.

Before concluding this Section, we elaborate on our assumption that the
averaging over phase fluctuations should include low-energy states from the
continuum spectra of the condensation energy. Consider the case $\gamma<2$,
when the spectrum is still discrete and the $n=0$ solution has the lowest
condensation energy $E_{c,0}$. The energies $E_{c,n\geq 1}$ are close to
$E_{c,0}$, yet the solutions with different $n$ are topologically distinct as
$\Delta_{n}(\omega_{m})$ has $n$ vortices. These other states contribute to
the renormalization of the phase of $\Delta_{0}(\omega_{m})$ only if the
tunneling amplitude between the states $n=0$ and $n>0$ is non-zero, which
requires the barrier between $E_{c,0}$ and $E_{c,n}$ to be small. The height
of the barrier depends on the path along which a state without a vortex
transforms into a state with a vortex at some small $\omega_{m}$. A vortex can
either come from $\omega_{m}=\infty$, in which case the barrier is high, or
via a creation of a vortex-antivortex pair at $\omega_{m}=0$, in which case it
is low. For a generic $\gamma<2$, $\Delta_{n}(z)$ are regular at small $z$ in
the complex plane, hence one should not expect an anti-vortex nearby. However,
for $\gamma\to 2$, our candidate $\Delta_{\xi}(z)$, Eq. (38), possess anti-
vortices at small $z$ in the lower frequency half-plane. In this situation, it
is natural to expect that , the barriers between $E_{c,0}$ and $E_{c,n}$ with
$n>0$ are low, hence our reasoning is justified.

##  VIII Phase diagram of the $\gamma$ model

###  VIII.1 $\gamma=2$, finite $\omega_{D}$

In Fig. 15, we present the phase diagram for $\gamma=2$ for non-zero
$\omega_{D}$ and $T$. At $\omega_{D}=0$, the true transition temperature into
a SC state is zero, although the onset temperature for the pairing, $T_{p}$ is
finite. At finite $\omega_{D}$, $T_{c}$ is finite but much smaller than
$T_{p}$, at least for small $\omega_{D}$. In between $T_{c}$ and $T_{p}$ the
system displays pseudogap behavior: the spectral function and the DoS display
a peak at a finite frequency, but the spectral weight below the peak remains
finite. Close to $T_{p}$, the pairing is mainly induced by fermions with the
two lowest Matsubara frequencies $\pm\pi T$[46, 47, 48]. In this situation,
the position of the peak in the spectral function and the DoS increases
linearly with $T$, and the gap fills in as $T$ approaches $T_{p}$. In the $T$
range near $T_{c}$, fermions with all Matsubara frequencies contribute to the
pairing, and the positions of the maxima in the spectral function and the DoS
move to smaller frequencies as $T$ increases (gap closing behavior). We show
the DoS in the two regimes in Fig.18 below.

Figure 15: The phase diagram of the $\gamma$ model for $\gamma=2$ in variables
$(T/{\bar{g}},\omega_{D}/{\bar{g}})$, where $\omega_{D}$ is the mass of a
pairing boson. $T_{p}$ is the onset temperature of the pairing, and $T_{c}$ is
the actual superconducting transition temperature, below which the system
establishes phase coherence. In between the system displays pseudogap
behavior, in which fermionic pairs are formed, but there is no macroscopic
phase coherence. The dashed line separates the two regimes within the
pseudogap phase – the one at higher $T$, where the system behavior is chiefly
determined by fermions with the two lowest Matsubara frequencies $\pm\pi T$,
and the one at lower $T$, when fermions with all Matsubara frequencies
contribute to the pairing. In these two regimes the system displays gap
filling and gap closing behavior, respectively. Figure 16: The phase diagram
of the $\gamma$ model for a generic $\gamma<2$ at a finite $T$ and vanishing
$\omega_{D}$. For any $\gamma<2$, the true SC transition temperature $T_{c}$
is finite, but is smaller than the onset temperature for the paring, $T_{p}$.
In between $T_{c}$ and $T_{p}$, the system displays a pseudogap behavior.
There are two distinct behaviors in the pseudogap regime, like in Fig. 15:
close to $T_{p}$, the spectral function and the DoS display gap filling
behavior, while close to $T_{c}$, the behavior becomes more conventional and
the gap frequency shifts to a smaller value as $T$ increases.

###  VIII.2 $\omega_{D}=0$, $0<\gamma\leq 2$

In Fig. 16 we show the phase diagram for the $\gamma$ model with $0<\gamma\leq
2$, at $\omega_{D}=0$ and finite $T$. This phase diagram is based on the
results of this work and previous works (Papers I-IV). For $\gamma<2$, we
found earlier the largest condensation energy is for sign-preserving solution
of the gap equation ($n=0$ in our classification). Still, for any $\gamma>0$,
there exists an infinite set of topologically distinct solutions for the gap
(all with the same symmetry), labeled by integer $n$. This generates a
discrete spectrum of the condensation energy $E_{c,n}$. The spectrum is sparse
near the bottom at small $\gamma$, but becomes dense and flattens up at the
bottom as $\gamma$ approaches $2$. At $\gamma\leq 2$, the corrections to
superconducting order parameter from the states with $n\neq 0$ are small at
low $T$, but rapidly increase with increasing $T$ and destroy phase coherence
at $T_{c}\sim{\bar{g}}(2-\gamma)/\lambda_{E}$. For $\gamma\leq 2$, $T_{c}\ll
T_{p}$, and there exists a wide intermediate temperature range where the
system displays a pseudogap behavior. By continuity, we expect that the
pseudogap region to exist for all $\gamma>0$ albeit with a smaller width.

###  VIII.3 Properties of the pseudogap phase

####  VIII.3.1 toy model for $\gamma=2$

Let’s start with $\gamma=2$. At $T=0$ the DoS is the set of
$\delta-$functions, (Fig.17(a)) At a finite $T$, two new features appear.
First, $\Delta_{0}(\omega)$ decreases with increasing $\omega$ and displays no
oscillations above $\omega_{max}$, similar to the case with finite
$\omega_{D}$ discussed in VI. As a result, $\delta$-functional peaks in the
DoS at larger frequencies get broadened and eventually disappear. Second,
other $D_{\xi}(\omega)$ from the continuum spectrum of condensation energies
contribute to the DoS with Boltzmann factors. For all these solutions,
Im$\Delta_{\xi}(\omega)$ remains finite down to $\omega=0$. As a result, the
DoS also becomes non-zero at the smallest $\omega$ (this phenomenon is often
called a gapless superconductivity [49, 50, 34, 35]) We model both effects by
introducing a phenomenological
$\Delta(\omega)=\omega/\sin(ia+(\omega/{\bar{g}})^{2}(1+ib))$, where $a$ and
$b$ increase with $T$. We show the corresponding DoS in Fig. 17(b) and (c).

Figure 17: The density of states, $N(\omega)$, at different temperatures, for
a toy model with $\Delta_{0}(\omega)=\omega/\sin(ia+\omega^{2}(1+ib))$, where
$a$ and $b$ are two parameters, which increase with $T$. (a) The $T=0$ limit,
$a=b=10^{-4}$. The DoS has a set of $\delta$-functional peaks. (b) A finite
but small $T$, $a=b=0.05$. The first few peaks are well defined, but the peaks
at large frequencies get overdamped and disappear. (c) A higher temperature,
$a=b=0.25$. The peak at the smallest frequency is still present, at about the
same frequency as at $T=0$, but other peaks are washed out, and the spectral
weight below the peak increases, i.e., the DoS at low frequencies fulls in.

####  VIII.3.2 gap filling vs gap closing

We argue, based on earlier works[47, 46]), that there are two different
regimes of system behavior within the pseudogap phase. At low $T$, the
position of the peak in the DoS scales with $\Delta_{0}(0)$ and decreases as
$T$ increases (the gap “closes” with increasing $T$). At higher $T$, the peak
in the DOS shifts to higher frequencies and the spectral weight below the peak
increases (gap “fills in” with increasing $T$). We illustrate this in Fig.18.
This last behavior is at least partly related to the fact that in some finite
range of $T$ below $T_{p}$, the gap function on the Matsubara axis is strongly
peaked at the first Matsubara frequencies $\pm\pi T$ (Refs. [47, 46]). On the
real axis the corresponding $\Delta(\omega)$ displays $\omega/T$ scaling. For
such $\Delta(\omega)$, the peak frequency in the DoS increases linearly with
$T$.

Figure 18: The temperature evolution of the DoS $N(\omega)$. For $\gamma<2$
(upper panel) there is a SC order at $T<T_{c}$. In this regime and in the
pseudogap state at $T\geq T_{c}$, the temperature variation of $N(\omega)$
resembles that in a conventional BCS superconductor, i.e. when $T$ increases,
the position of the maximum in $N(\omega)$ moves to a smaller frequency. At
larger $T$ within the pseudogap phase, $N(\omega)$ displays gap filling
behavior when the peak position increases with increasing $T$ and
$N(\omega=0)$ increases towards its normal state value. For $\gamma=2$ (lower
panel), $T_{c}=0$, but the two different regimes of pseudogap behavior are
present.

At a finite $\omega_{D}$ and/or $2-\gamma$, the “gap filling” behavior holds
in some range between the onset temperature of the pairing $T_{p}$ and a
finite superconducting $T_{c}$ (Fig.18). To estimate the crossover temperature
between the two regimes, we compare the actual $T_{p}$ with the one obtained
by neglecting the contributions from fermions with $\omega_{m}=\pm\pi T$. We
show the results in Fig. 19. We see that for $\gamma=2$ the onset temperature
without $\pm\pi T$ fermions is strongly reduced – it is about $1/7$ of the
actual $T_{p}\sim 0.18{\bar{g}}$. This implies that the “gap filling” behavior
holds in a wide range below $T_{p}$ and crosses over to “gap closing” behavior
only near $T_{c}$.

Figure 19:  The onset temperature of the pairing, obtained without including
Matsubara frequencies $\omega_{m}=\pm\pi T$. For $\gamma=2$, this temperature
is roughly $1/7$ of the actual $T_{p,0}$.

##  IX Conclusions

In this paper, we extended our earlier analysis of the $\gamma-$model to
$\gamma=2$. The $\gamma=2$ model describes, among other cases, the pairing,
mediated by an Einstein boson, in the limit when the bosonic mass $\omega_{D}$
tends to zero. On the real axis, the effective interaction in this limit
$V(\Omega)=-{\bar{g}}^{2}/\Omega^{2}$ is repulsive, and, at a first glance,
should not give rise to pairing. However, the same interaction on the
Matsubara axis, $V(\Omega_{m})={\bar{g}}^{2}/\Omega^{2}_{m}$, is attractive,
and earlier calculations on the Matsubara axis found that the onset
temperature of the pairing, $T_{p}$, tends to a finite value
$T_{p}=0.1827{\bar{g}}$ at $\omega_{D}\to 0$
($T_{p}=0.1827\omega_{D}\sqrt{\lambda}$ in terminology of Ref.[16], which is
the same expression because $\lambda={\bar{g}}^{2}/\omega^{2}_{D}$). The issue
we discussed in this paper is whether this $T_{p}$ is close to the actual
superconducting $T_{c}$, or $T_{c}$ is smaller, and there is a range of
pseudogap behavior between $T_{c}$ and $T_{p}$. We argued that the actual
$T_{c}$ scales with $\omega_{D}$ and is much smaller than $T_{p}$ when
$\omega_{D}/{\bar{g}}$ is small.

To prove this, we solved the non-linear gap equation at $T=0$ and
$\omega_{D}=0$ and found a continuum of solutions, governed by a single
parameter $\xi$ ($0\leq\xi\leq\infty$). This in turn gives rise to a continuum
spectrum of condensation energy, $E_{c,\xi}$, which can be viewed as a
continuum gapless spectrum of “longitudinal” gap fluctuations. An infinite set
of the gap functions and the condensation energies exists already for
$\gamma<2$, but is a discrete one. For $\gamma=2$, this spectrum becomes
continuous in a manner similar to how a discrete set of energy levels in a
finite size crystal becomes a continuous vibration spectrum when system size
becomes infinite. In our case, $1/(2-\gamma)$ plays the role of a system size.

Without the contribution from the gapless longitudinal branch, superfluid
stiffness $\rho_{s}(T=0)$ is larger than $T_{p}$, and thermal corrections to
superconducting order parameter scale approximately as $T/\rho_{s}(0)$ and
remain small at all $T<T_{p}$. However, upon including contributions from the
longitudinal branch, we found that thermal corrections become of order one
already at much smaller $T\sim\omega_{D}$. We identified this temperature with
the actual superconducting $T_{c}$. We emphasize that $T_{c}$ vanishes at
$\omega_{D}=0$, and the behavior of the stiffness depends on the order in
which the double limit $\omega_{D}\to 0$ and $T\to 0$ is taken. This strongly
suggests that the $\gamma=2$ model is critical at $T=0$. At smaller $\gamma$,
the ground state is not critical at $\omega_{D}=0$, and
$T_{c}\sim{\bar{g}}(2-\gamma)$. It is finite but at $\gamma\leq 2$ is still
much smaller than $T_{p}\sim{\bar{g}}$.

We presented collaborative evidence that the $\gamma=2$ model is critical,
from the analysis of the continuum set of gap functions along real frequency
axis and in the upper half-plane of frequency. We found that for each
solution, there is an infinite array of $2\pi$ vortices in the upper frequency
half-plane. The array of vortices stretches up to an infinite frequency, where
each gap function from the continuous set has an essential singularity. We
speculated that different gap functions from the continuous set are different
extensions from the array of vortices, onto the upper half-plane of frequency.

At a finite $\omega_{D}$, the set of gap functions becomes discrete and
contains only a finite number of solutions, all of which behave regularly in
the high-frequency limit. The the number of vortices also becomes finite.
Still, at small $\omega_{D}/{\bar{g}}$, the system behavior over a wide
frequency range mimics that at $\omega_{D}=0$.

We showed the phase diagram of the $\gamma=2$ model in variables $T$ and
$\omega_{D}$ in Fig. 15 and the phase diagram of the $\gamma$ model at
$\omega_{D}=0$ in in variables $T$ and $\gamma$ in Fig. 16 In both cases,
there is range of pseudogap behavior between the onset temperature of the
pairing $T_{p}$ and the actual $T_{c}$. In the pseudogap region, the bound
pairs are formed, but there is no macroscopic phase coherence. We argued that
in most of the pseudogap regime, the DoS and other observables display “gap
filling” behavior, in which the peak position remains at a finite frequency up
to $T_{p}$, while the states below the peak gradually fill in.

In the next (last) paper in the series we consider the behavior of the
$\gamma$ model for $\gamma>2$ and show that the new physics emerges at $T=0$,
which gives rise to a reduction and eventual vanishing of the superfluid
stiffness in the ground state.

###### Acknowledgements.

We thank I. Aleiner, B. Altshuler, E. Berg, D. Chowdhury, L. Classen, K.
Efetov, R. Fernandes, A. Finkelstein, E. Fradkin, A. Georges, S. Hartnol, S.
Karchu, S. Kivelson, I. Klebanov, A. Klein, R. Laughlin, S-S. Lee, G.
Lonzarich, I. Esterlis, D. Maslov, F. Marsiglio, I. Mazin, M. Metlitski, W.
Metzner, A. Millis, D. Mozyrsky, C. Pepan, V. Pokrovsky, N. Prokofiev, S.
Raghu, S. Sachdev, T. Senthil, D. Scalapino, Y. Schattner, J. Schmalian, D.
Son, G. Tarnopolsky, A-M Tremblay, A. Tsvelik, G. Torroba, E. Yuzbashyan, J.
Zaanen, and particularly R. Combescot and Y. Wang for useful discussions. The
work by A.V.C. and Y.M.W. was supported by the NSF DMR-1834856. Y.-M.W,
S.-S.Z.,and A.V.C acknowledge the hospitality of KITP at UCSB, where part of
the work has been conducted. The research at KITP is supported by the National
Science Foundation under Grant No. NSF PHY-1748958. A.V.C. also acknowledges
the hospitality of Stanford University, where some results of this work have
been obtained. His stay at Stanford has been supported through the Gordon and
Betty Moore Foundation’s EPiQS Initiative, Grants GBMF4302 and GBMF8686.

##  Appendix A Expansion in $D^{2}(\omega_{m})$ for $\gamma=2$ and $\gamma<2$

In this Appendix we present some details of the analysis of the non-linear gap
equation for $\gamma=2$ and elaborate on the claim in the main text that a
continuous set of gap functions exists only for $\gamma=2$, while for smaller
$\gamma$, the set is a discrete one.

###  A.1 $\gamma=2$

We begin with $\gamma=2$. Consider first the limit of small frequencies
$\omega_{m}\ll{\bar{g}}$. For such $\omega_{m}$, $\Delta(\omega_{m})$ in the
l.h.s. of the gap equation (4) can be neglected, as its inclusion leads to
terms with extra $(\omega_{m}/{\bar{g}})^{2}$. This approximation is
equivalent to neglecting $\omega_{m}$ compared to the self-energy
$\Sigma(\omega_{m})$ and is similar to the “no $\omega_{m}$” approximation,
used in the studies of SYK-type models [51, 52, 53]. The non-linear gap
equation at $T=0$ without $\Delta(\omega_{m})$ in the l.h.s reduces to

$\int
d\omega_{m^{\prime}}\frac{D(\omega_{m^{\prime}})-D(\omega_{m})}{\sqrt{1+D^{2}(\omega_{m^{\prime}})}}~{}\frac{\operatorname{sgn}\omega^{\prime}_{m}}{|\omega_{m}-\omega_{m^{\prime}}|^{2}}=0$
(92)

where, we recall, $D(\omega_{m})=\Delta(\omega_{m})/\omega_{m}$.

The linearized gap equation is obtained from (92) by neglecting
$D^{2}(\omega_{m^{\prime}})$ in the denominator. The exact solution of the
linearized gap equation is Eq. (13):

$D(\omega_{m})=2\epsilon\cos{\left(\beta\log{\left(\frac{|\omega_{m}|}{{\bar{g}}}\right)^{2}}+\phi_{\epsilon}\right)}~{}\operatorname{sgn}\omega_{m}$
(93)

where $\epsilon$ is an arbitrary overall factor, $\phi_{\epsilon}$ is yet
undetermined constant, and $\beta=0.38187$ satisfies
$\pi\beta\tanh(\pi\beta)=1$.

We now expand Eq. (92) in powers of $D^{2}$. We will be searching for the
solution in the form

$D_{\epsilon}(\omega_{m})=2\sum_{n=0}^{\infty}\epsilon^{2n+1}Q_{2n+1}\cos{\left((2n+1)\left(\beta_{\epsilon}\log{\left(\frac{|\omega_{m}|}{{\bar{g}}}\right)^{2}}\right)+\phi_{\epsilon}\right)}$
(94)

Substituting into (92) and collecting contributions at each order in
$\epsilon^{2n+1}$, we find that $D_{\epsilon}(\omega_{m})$ given by (94) does
satisfy Eq. (92), and that all integrals are ultra-violet convergent, i.e.,
there is no need for regularization. The calculations are lengthy, but
straightforward. We checked explicitly that $\beta_{\epsilon}$ is the same in
all terms in (94) and is related to the original $\beta$ by

$\beta_{\epsilon}=\beta\left(1-\frac{\epsilon^{2}}{2}+0.806\epsilon^{4}+...\right)$
(95)

The numerical coefficients are $Q_{3}=0.222+O(\epsilon^{2})$,
$Q_{5}=0.043+O(\epsilon^{2})$. We cited this result and Eqs. (94) and (95) in
Sec. III.3.

At larger frequencies, we need to keep $\Delta(\omega_{m})$ in the l.h.s. of
(4). In the opposite limit $\omega_{m}\gg{\bar{g}}$, the leading term in
$D_{\epsilon}(\omega_{m})$ is obtained by pulling $1/\omega_{m}^{2}$ from the
integrand in the r.h.s.. Then we obtain

$D_{\epsilon}(\omega_{m})=\frac{a_{\epsilon}}{\omega^{3}_{m}}$ (96)

where

$a_{\epsilon}=\frac{{\bar{g}}^{2}}{2}\int
d\omega_{m^{\prime}}\frac{D(\omega_{m^{\prime}})}{\sqrt{1+D^{2}(\omega_{m^{\prime}})}}~{}\operatorname{sgn}\omega_{m}$
(97)

Substituting this form of $D(\omega_{m})$ into the integrand in the r.h.s, we
find that the integral is ultra-violet convergent, i.e., the solution (97) is
self-consistent.

The two solutions have to merge at $\omega_{m}\sim{\bar{g}}$. For the
linearized gap equation (the limit $\epsilon\to 0$) we verified that this does
happen for a certain value of $\phi_{\epsilon}$ in (94). We conjecture that
the same holds for other $\epsilon$, i.e., for a certain $\phi_{\epsilon}$,
$D_{\epsilon}(\omega_{m})$ smoothly evolves between (94) and (97). We did
similar analysis in Paper I. There, we demonstrated that for arbitrary
$\phi_{\epsilon}$, $D_{\epsilon}(\omega)$ of Eq. (94) approaches the constant
at $\omega_{m}\to\infty$, while the desired term
($D_{\epsilon}(\omega_{m})\propto 1/|\omega_{m}|^{\gamma+1}$ for a generic
$\gamma$) is the subleading one. For a particular $\phi_{\epsilon}$, a
constant vanishes, and the high-frequency behavior becomes the expected one.

We see from (95) that $\beta_{\epsilon}$ decreases with increasing $\epsilon$,
while the overall magnitude of $\Delta(\omega_{m})$ increases. It is natural
to expect that $\beta_{\epsilon}=0$ at some critical $\epsilon=\epsilon_{cr}$.
We explored this in the main text.

Another way to argue for the existence of $\epsilon_{cr}$ is to depart from
the opposite limit $\epsilon\gg 1$, where $D(\omega_{m})$ is supposed to be
large. In this case, we introduce $\Xi(\omega_{m})=1/D(\omega_{m})$ and re-
express the gap equation as

$\int
d\omega_{m^{\prime}}\frac{\Xi(\omega_{m^{\prime}})-\Xi(\omega_{m})}{\sqrt{1+\Xi^{2}(\omega_{m^{\prime}})}}~{}\frac{1}{|\omega_{m}-\omega_{m^{\prime}}|^{2}}=0$
(98)

Note the absence of $\operatorname{sgn}\omega^{\prime}_{m}$ in the integrand.
At small $\Xi$, we neglect the $\Xi^{2}(\omega_{m^{\prime}})$ in the
denominator and search for the solution in the form
$\Xi(\omega_{m})=\mbox{sign}\omega_{m}|\omega_{m}/{\bar{g}}|^{b}$.
Substituting into (98) we find $b=\pm 1$, i.e.,

$\Xi(\omega_{m})=\left(A_{1}\frac{{\bar{g}}}{|\omega_{m}|}+A_{2}\frac{|\omega_{m}|}{{\bar{g}}}\right)$
(99)

The first term does not satisfy the normalization condition and has to be
discarded (see Paper I for the details on this). This leaves no parameter to
adjust in order to match with the behavior at high frequencies. This implies
that there is no solution for the gap at large $\epsilon$.

There is a similarity between this analysis and the analysis in Paper I, where
we considered the $\gamma$ model for $\gamma<1$ and extended it using a
continuous variable $N$ to make interactions in the particle-hole and
particle-particle channels non-equivalent. There, we found that there exists
$N_{cr}$, which separates oscillating and non-oscillating solutions, and only
oscillating solutions are compatible with high-frequency behavior. Here,
$\epsilon$ plays the same role as $N$. We illustrate this in Fig. 20

Figure 20: The comparison between the behavior of $\beta_{\epsilon}$ in the
$\gamma=2$ model and $\beta_{N}$ in the model with $\gamma<1$, extended to
$N>1$.

###  A.2 $\gamma<2$

We now extend this approach to $\gamma<2$. The gap equation for
$D(\omega_{m})$ at $\omega_{m}\ll{\bar{g}}$ has the same form as in (92), only
$|\omega_{m}-\omega_{m^{\prime}}|^{2}$ in the denominator is replaced by
$|\omega_{m}-\omega_{m^{\prime}}|^{\gamma}$. The solution of the linearized
equation for $D(\omega_{m})$ is

$D(\omega_{m})=2\epsilon\left(\frac{{\bar{g}}}{|\omega_{m}|}\right)^{1-\gamma/2}\cos{\left(\beta_{\gamma}\log{\left(\frac{|\omega_{m}|}{{\bar{g}}}\right)^{\gamma}}+\phi\right)}~{}\mbox{sign}\omega$
(100)

where $\beta_{\gamma}$ is some regular function of $\gamma$. As before, we
search for the solutions in the form

$D_{\epsilon}(\omega_{m})=2\sum_{n=1}^{\infty}\left(\epsilon\left(\frac{{\bar{g}}}{|\omega_{m}|}\right)^{1-\gamma/2}\right)^{2n+1}Q_{2n+1}\cos\left[(2n+1)\left(\beta_{\gamma,\epsilon}\log{\left(\frac{|\omega_{m}|}{{\bar{g}}}\right)^{\gamma}}+\phi_{\epsilon}\right)\right]$
(101)

Substituting into (100), we find that the integrals that determine $Q_{2n+1}$
now contain infra-red divergencies. The only way to eliminate the divergencies
is to assume that $\Delta_{\epsilon}$ tends to a finite value at
$\omega_{m}\to 0$. But this is only possible for a discrete set of finite
$\epsilon$. We also note in passing that because the actual expansion
parameter is $\epsilon({\bar{g}}/|\omega_{m}|)^{1-\gamma/2}$, the expansion of
$\beta_{\gamma,\epsilon}$ in powers of $\epsilon$ yields
$\beta_{\gamma,\epsilon}=\beta_{\gamma}(1+a(\epsilon({\bar{g}}/|\omega_{m}|)^{1-\gamma/2})^{2})$.
For $a\neq 0$, this gives rise to additional terms, which are not matched by
the terms in Eq. (101). The only option then is to set $a=0$, i.e., leave
$\beta_{\gamma,\epsilon}$ equal to its bare value $\beta_{\gamma}$.

The outcome is that the continuous set of gap functions exists only for
$\gamma=2$. For smaller $\gamma$, this set is discrete. We also emphasize that
the distinction between $\gamma=2$ and $\gamma<2$ holds only at $T=0$. At any
finite $T$, the set of gap functions is a discrete one for all $\gamma\leq 2$,
and the solutions with different $n$ from the set vanish at different
temperatures $T_{p,n}$. In Fig. 21 we show the results of high-accuracy
numerical calculation of $T_{p,1}$ for $\gamma=2$. We see that $T_{p,1}\simeq
3.6827\times 10^{-3}{\bar{g}}$ is finite.

Figure 21: Determination of the temperature $T_{p,1}$, at which the $n=1$
solution develops in the $\gamma=2$ model. $M_{c}\sim 10^{3}$ is the largest
value of the Matsubara number, used in this numerical calculation.
Extrapolating $M_{c}$ to $\infty$ yields a finite value $T_{p,1}\simeq
3.6827\times 10^{-3}{\bar{g}}$.

##  Appendix B The gap function $\Delta_{\infty}$

The exact solution of the linearized gap equation at zero temperature has been
derived for $0<\gamma<1$ in Paper I and $1<\gamma<2$ in Paper IV. Here we
extend the analysis of Paper IV to $\gamma=2$.

We first solve for the gap function $\Delta_{\infty}(\omega_{m})$ along the
Matsubara axis. Following the same computational steps as in the analysis for
$\gamma<2$, we obtain
$D_{\infty}(\omega_{m})=\Delta_{\infty}(\omega_{m})/\omega_{m}$ in the form

$D_{\infty}(\omega_{m})=\epsilon\frac{{\bar{g}}}{\omega_{m}}\int_{-\infty}^{\infty}dkb_{k}e^{-ik\log{(\omega_{m}/{\bar{g}})^{2}}},$
(102)

where $\epsilon$ is an infinitesimal number,

$b_{k}=\frac{e^{-iI_{k}}}{\left[\cosh(\pi(k-\beta))\cosh(\pi(k+\beta))\right]^{1/2}},$
(103)

and

$I_{k}=\frac{1}{2}\int_{-\infty}^{\infty}dk^{\prime}\log{|\epsilon_{k^{\prime}}-1|}\tanh{\pi(k^{\prime}-k+i\delta)},$
(104)

Here $\epsilon_{k^{\prime}}=\pi k^{\prime}\tanh(\pi k^{\prime})$ and
$\beta\simeq 0.38187$ is the solution of $\pi\beta\tanh{(\pi\beta)}=1$. We
cited these results in Eqs. (15), (16) and (17) in the main text.

The integrals (104) and (102) can be computed numerically. We showed the
result for $D_{\infty}(\omega_{m})$ in Fig.1 in the main text. The function
$D_{\infty}(\omega_{m})$ oscillates at $\omega_{m}<{\bar{g}}$ and decays as
$1/|\omega_{m}|^{3}$ at $\omega_{m}>{\bar{g}}$.

###  B.1 Series expansion

The integral in Eq. (104) can be evaluated by closing the integration contour
along an infinite arc in the complex plane of frequency. For
$\rvert\omega_{m}\rvert<{\bar{g}}$, the arc must be in the upper half-plane,
and for $|\omega_{m}|>{\bar{g}}$, in the lower half-plane. The integral is
equal to the sum of the contributions from each pole of the function $b_{k}$
in the upper or lower half-plane. The position of these poles are obtained
from the representation of $b_{k}$ as an infinite product of the Gamma-
functions: (see Papers I and IV for details on this)

$b_{k}=\frac{\Gamma(1-ik)}{\Gamma(1+ik)}\Gamma\left(\frac{1}{2}+i(k+\beta)\right)\Gamma\left(\frac{1}{2}+i(k-\beta)\right)\prod_{m=1}^{\infty}\frac{\Gamma\left(\frac{1}{2}+i(k-i\beta_{m})\right)\Gamma\left(1+m-ik\right)}{\Gamma\left(\frac{1}{2}-i(k+i\beta_{m})\right)\Gamma\left(1+m+ik\right)}$
(105)

Here $\beta_{m}>0$ are the solutions of $\pi\beta_{m}\tan({\pi\beta_{m}})=-1$.
There is an infinite set of such $\beta_{m}$, specified by an integer
$m=0,1,2..$. Each $\beta_{m}$ is located within the interval
$1/2+m<\beta_{m}<3/2+m$. Viewed as a function of complex $k$, $b_{k}$ has
poles from individual $\Gamma$-functions in the upper half-plane at
$k=\pm\beta+i(n+1/2)$ and $k=i\beta_{m}+i(n+1/2)$ where $n,m=0,1,2,…$, and in
the lower half-plane, at $k=-i(n+1)$ and $k=-i(1+m+n)$, where $n=0,1,2,...$
and $m=1,2,...$.

###  B.2 $\rvert\omega_{m}\rvert<{\bar{g}}$

For $\rvert\omega_{m}\rvert<{\bar{g}}$, the relevant poles are at
$k=\pm\beta+i(n+1/2)$ and at $k=i\beta_{m}+i(n+1/2)$, $n=0,1,2,...$. This
yields series expansion for $D_{\infty}(y)$ with
$y=(\rvert\omega_{m}\rvert/{\bar{g}})^{2}$ in the form

$D_{\infty}(y)={\text{R}e}~{}\sum_{n=0}^{\infty}e^{(i\beta\log{y}+\phi)}C^{<}_{n}~{}y^{n}+\sum_{n,m=0}^{\infty}D^{<}_{n,m}y^{(n+\beta_{m})}$
(106)

The leading term in (106) at small $\omega_{m}$ comes from the contribution of
the poles at $k=\pm\beta+i/2$

$D_{\infty}(y)=C^{<}_{0}\cos{\left(\beta\log{y}+\phi\right)}$ (107)

The first subleading term comes from the contribution of the pole at
$k=i(1/2+\beta_{0})$ and scales as $y^{\beta_{0}}$, where $\beta_{0}\simeq
0.89$.

In the direct perturbation expansion in $y$, the series in $y^{n}$ (the first
term in (106)) come from fermions with internal $y^{\prime}\sim y$ and form
the “local” series. The second term in (106)) is the sum of contributions from
fermions with $y^{\prime}=O(1)$, which for $y\ll 1$ can we regarded as “non-
local”. The total $D_{\infty}(y)=D_{\infty,L}(y)+D_{\infty,NL}(y)$.

The coefficients $C^{<}_{n}$ in (106) can be obtained analytically, as we
already found in Papers I and III for $\gamma\leq 1$ and Paper IV for
$1<\gamma<2$. At $\gamma=2$, the result takes a very simple form

$C^{<}_{n}=C^{<}_{0}{i^{n}\over\beta^{n}n!}$ (108)

Substituting this into Eq. (106), we find that the first term (the local
contribution) becomes

$D_{\infty,L}(y)\propto\cos{\left[\beta(\log{y}-y)+\phi\right]}.$ (109)

It oscillates with the periodicity set by $\beta\log y$ for $y\ll 1$, i.e.,
$\rvert\omega_{m}\rvert\ll{\bar{g}}$, which is the right behavior of the gap
function at small frequencies, see (107).

We note, however, that the first subleading term in (109) scales as
$y~{}\sin({\beta\log{y}})$. This contribution is smaller than the actual
subleading term, which scales as $y^{0.89}$ and does not oscillate. This
implies that, besides the leading term, the form of $\Delta_{\infty}(y)$ is
determined by non-local corrections.

###  B.3 $\rvert\omega_{m}\rvert>{\bar{g}}$ and logarithmic correction

For $y>1$, i.e., $\rvert\omega_{m}\rvert>{\bar{g}}$, relevant poles are in the
lower half-plane. According to Eq. (105), a pole at $k=-i(n+1)$
($n=0,1,2,...$) is of order $n+1$, namely a simple pole at $n=0$, a double
pole at $n=1$, etc. The leading term in the limit of
$\rvert\omega_{m}\rvert\to\infty$ is the contribution from a simple pole at
$k=-i$ ($n=0$), This contribution accounts for $1/y$ behavior of
$\Delta_{\infty}(y)$ at large $y$. However, the subleading terms from the rest
poles contain extra logarithms on top of powers of $1/y$:

$\Delta_{\infty}(y)=\sum_{n=0}^{\infty}{\tilde{C}}^{>}_{n}~{}y^{-2(1+n)}\left(\log{y}\right)^{n},$
(110)

To demonstrate the presence of the logs, consider as an example the
contribution from the double pole at $k=-2i$. We shift $\gamma$ to $2-\delta$,
$\delta>0$ and then take the limit $\delta\to 0$. The expression of $b_{k}$
for $\gamma\leq 2$ is presented in Paper IV. Using it, we find that a double
pole splits into two simple poles at $z_{1}=-2i$ and $z_{2}=-(2+\delta/2)i$.
In the neighborhood of the two poles, the function $b_{k}$ takes the form
$\sim 1/(z-z_{1})/(z-z_{2})$. The contribution from the these two poles is
obtained by circling out a loop ${\cal C}$ enclosing $z_{1}$ and $z_{2}$.
Evaluating the integral and taking the limit $\delta\to 0$, we obtain

$\displaystyle\sqrt{y}\lim_{\delta\to
0^{+}}\frac{{\bar{g}}}{\omega_{m}}\oint_{{\cal
C}}dz{1\over(z+2i)(z+(2+\delta/2)i)}e^{-iz\log(y^{1-\delta/2})}=2\pi\frac{\log{y}}{y^{2}}.$
(111)

Similarly, the triple pole at $k=-3i$ gives rise to $(\log{y})^{2}/y^{3}$,
etc. Collecting the contributions from every pole on the lower-half plane, we
obtain (110).

###  B.4 The universal oscillating term at large $y$

We now show that the high-frequency form of $D_{\infty}(\omega_{m})$ contains
an additional oscillating contribution. This contribution is exponentially
small on the Matsubara axis, but, as we will see, it becomes the dominant one
on the real axis. To extract this contribution, we note that for large
$\rvert\omega_{m}\rvert/{\bar{g}}$, the argument of the cosine function,
$I_{k}+k\log{y}$, passes through extremum at $k\sim k_{*}=y/\pi$. Expanding
around this point and evaluating the Gaussian integral, we obtain the
universal piece $D_{\infty;u}(y)$ in the form

$D_{\infty;u}(y)=2\sqrt{2}\epsilon
e^{-y}\cos{\left[\frac{(\pi^{2}-2)}{2\pi}y+\frac{\pi}{4}\right]}.$ (112)

We see that $D_{\infty;u}(y)$ is exponentially small, yet this oscillating
term is present. The total $D_{\infty}(y)$ is the sum of (110) and (112).

###  B.5 $D_{\infty}(y)$ along real axis

Let’s now transform from Matsubara to real axis. We use $\omega$ instead of
$y$ for better transparency. By construction, the gap function along the real
axis is obtained by replacing $i\omega_{m}\to\omega+i0^{+}$ in the integrand
in the r.h.s. of Eq. (102). Under this transformation,
$\log(|\omega_{m}|/{\bar{g}})^{2}$ transforms into
$\log(\rvert\omega\rvert/{\bar{g}})^{2}-i\pi$. The integral in Eq. (102)
splits into two parts:

$D_{\infty}(\omega)=\epsilon\frac{{\bar{g}}}{\omega}\int_{0}^{\infty}dk\frac{e^{-\pi
k}e^{-iI_{k}-ik\log(\rvert\omega\rvert/{\bar{g}})^{2}}+e^{\pi
k}e^{iI_{k}+ik\log(\rvert\omega\rvert/{\bar{g}})^{2}}}{\sqrt{\cosh(\pi(k-\beta))\cosh(\pi(k+\beta))}}.$
(113)

Evaluating each integral by expanding near the point where $I_{k}\pm
k\log(\rvert\omega\rvert/{\bar{g}})^{2}$ passes through extremum and
approximating the denominator in (113) by its form at large $k$, we find that
the first term is small in $e^{-2\pi k}$, while in the second term the
exponential factor cancels out. Ignoring the first term, we obtain

$D_{\infty;u}(\omega)\approx\sqrt{2}\epsilon
e^{{i\over\pi}\left[\left({\omega\over{\bar{g}}}\right)^{2}+\log\left({\omega\over{\bar{g}}}\right)^{2}\right]}.$
(114)

Other contributions to $D_{\infty}(\omega)$ contain powers of
${\bar{g}}/\omega$ and are smaller. As a result, on the real axis,
$D_{\infty}(\omega)\approx D_{\infty;u}(\omega)$ at $\omega\gg{\bar{g}}$.

The same calculation can be carried out for an arbitrary complex frequency
$z=\omega^{\prime}+i\omega^{{}^{\prime\prime}}$ in the upper frequency plane.
For this, one has to replace $i\omega_{m}$ by $z\equiv\rvert z\rvert
e^{i\psi}$ ($0<\psi<\pi$) in the integrand in the r.h.s. of (102). This
changes $\log(|\omega_{m}|/{\bar{g}})^{2}$ to $\log(\rvert
z\rvert/{\bar{g}})^{2}+i(2\psi-\pi)$ and gives $\Delta_{\infty}(z)$, which we
presented in Eq. (62) in the main text.

##  Appendix C Extended $\gamma-$ model

In Papers I-III and other works [54, 48, 55, 56, 47, 46, 57, 58, 59], we and
others extended the $\gamma$ model to in-equal interactions in the particle-
particle and particle-hole channel. This was done by adding a factor $1/N$ to
the interaction in the particle-particle channel and leaving the interaction
in the particle-hole channel intact. The advantage of extending the model to
$N\neq 1$ is that superconducting order in the ground state exists for
$N<N_{cr}$, while for larger $N$ the ground state is a non-Fermi liquid. By
analyzing the gap equation near this point, one can obtain useful information
about how a discrete set of solutions emerges. In Paper III we argued that the
extension to $N\neq 1$ makes sense for $\gamma<1$, while for $\gamma\geq 1$,
the model with $N\neq 1$ possesses singularities, not present in the original
$\gamma$ model. We proposed another way to extend the model with $\gamma>1$,
which is free from singularities. The idea is to first explicitly cancel out
singularities in the original $\gamma$ model with $\gamma>1$, and only then
extend the model to $M\neq 1$ by making interactions in the particle-particle
and particle-hole channel in-equivalent. The extended model is then free from
singularities, and one can obtain critical $M_{cr}$, where superconducting
order disappears at $T=0$ (by our construction, it exists at $M>M_{cr}$). In
this Appendix, we analyze the extended $\gamma$ model for $\gamma=2$. We show