# Interplay between superconductivity and non-Fermi liquid at a quantum-
critical point in a metal.
IV: The $\gamma$ model and its phase diagram at $1<\gamma<2$.

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, USA
Andrey V. Chubukov  School of Physics and Astronomy and William I. Fine
Theoretical Physics Institute, University of Minnesota, Minneapolis, MN 55455,
USA

###### Abstract

In this paper we continue with our analysis of the interplay between the
pairing and the non-Fermi liquid behavior in a metal for a set of quantum-
critical (QC) systems with an effective dynamical electron-electron
interaction $V(\Omega_{m})\propto 1/|\Omega_{m}|^{\gamma}$, mediated by a
critical massless boson (the $\gamma$-model). In previous papers we considered
the cases $0<\gamma<1$ and $\gamma\approx 1$. We argued that the pairing by a
gapless boson is fundamentally different from BCS/Eliashberg pairing by a
massive boson as for the former there exists not one but an infinite discrete
set of topologically distinct solutions for the gap function
$\Delta_{n}(\omega_{m})$ at $T=0$ ($n=0,1,2...$), each with its own
condensation energy $E_{c,n}$. Here we extend the analysis to larger
$1<\gamma<2$. We argue that the discrete set of solutions survives, and the
spectrum of $E_{c,n}$ get progressively denser as $\gamma$ increases towards
$2$ and eventually becomes continuous at $\gamma\to 2$. This increases the
strength of ”longitudinal” gap fluctuations, which tend to reduce the actual
superconducting $T_{c}$ compared to the onset temperature for the pairing and
give rise to a pseudogap region of preformed pairs. We also detect two
features on the real axis, which develop at $\gamma>1$ and also become
critical at $\gamma\to 2$. First, the density of states evolves towards a set
of discrete $\delta-$functions. Second, an array of dynamical vortices emerges
in the upper frequency half-plane, near the real axis. We argue that these two
features come about because on a real axis, the real part of the dynamical
electron-electron interaction,
$V^{\prime}(\Omega)\propto\cos(\pi\gamma/2)/|\Omega|^{\gamma}$, becomes
repulsive for $\gamma>1$, and the imaginary
$V^{{}^{\prime\prime}}(\Omega)\propto\sin(\pi\gamma/2)/|\Omega|^{\gamma}$,
gets progressively smaller at $\gamma\to 2$. We speculate that the features on
the real axis are consistent with the development of a continuum spectrum of
the condensation energy, for which we used $\Delta_{n}(\omega_{m})$ on the
Matsubara axis. We consider the case $\gamma=2$ separately in the next paper.

##  I  Introduction.

In this paper, we continue our analysis of the competition between non-Fermi
liquid (NFL) physics and superconductivity (SC) near a quantum-critical point
(QCP) of an itinerant correlated electron system, whose low-energy properties
are described by an effective, momentum-integrated dynamical four-fermion
interaction $V(\Omega_{m})$. We consider a set of quantum-critical systems,
for which the effective interaction has the form
$V(\Omega_{m})={\bar{g}}^{\gamma}/|\Omega_{m}|^{\gamma}$ (the $\gamma$-model).
This interaction gives rise to NFL behavior in the normal state, with
$\Sigma(\omega_{m})\propto\omega^{1-\gamma}_{m}$, and at the same time
mediates an attraction in at least one pairing channel, specific to the
underlying microscopic model. The two tendencies compete with each other as a
NFL self-energy reduces the pairing strength, while the feedback from the
pairing reduces fermionic self-energy. In Ref. Abanov and Chubukov (2020) we
listed quantum-critical systems, whose low-energy physics is described by the
$\gamma$-model with different $\gamma$ and presented a long list of references
to earlier publications on this subject. In that and subsequent papers, Refs.
Wu  _et al._ (2020a, b), hereafter referred to as Papers I-III, respectively,
our group analyzed the behavior of the $\gamma$-model for $0<\gamma<1$ (at
$T=0$ in Paper I and at a finite $T$ in Paper II), and at around $\gamma=1$
(in Paper III). In Paper I, we found that the system does become unstable
towards pairing, i.e., in the ground state the pairing gap
$\Delta(\omega_{m})$ is non-zero. However, in a qualitative distinction with
BCS/Eliashberg theory of a pairing by a massive boson, we found an infinite
discrete set of solutions for the gap function, $\Delta_{n}(\omega_{m})$,
specified by an integer $n=0,1,2...$. All solutions have the same spatial gap
symmetry, but are topologically distinct in the sense that
$\Delta_{n}(\omega_{m})$ changes sign $n$ times as a function of Matsubara
frequency (we discuss topological aspects in this paper). A gap function
$\Delta_{n}(\omega_{m})$ remains approximately $\Delta_{n}(0)$ at
$\omega_{m}\leq\Delta_{n}(0)$, oscillates $n$ times at
$\Delta_{n}(0)<\omega_{m}<{\bar{g}}$, and at larger frequencies decays as
$1/|\omega_{m}|^{\gamma}$. The magnitude of $\Delta_{n}(0)$ decreases with $n$
and at large $n$ scales as $\Delta_{n}(0)\propto e^{-An}$, where $A$ is a
function of $\gamma$. In the limit $n\to\infty$, $\Delta_{\infty}$ is
infinitesimally small and is the solution of the linearized gap equation. In
Paper II, we found that $\Delta_{n}(\omega_{m})$ emerge at different onset
temperatures $T_{p,n}$. The onset temperature and the condensation energy at
$T=0$, $E_{c,n}$, is the largest for $n=0$, yet the existence of an infinite
set of minima of the free energy is a qualitatively new feature of the pairing
at a QCP. In Paper III, we argued that the system behavior is continuous
through $\gamma=1$, despite that the fermionic self-energy and the pairing
vertex diverge at $T=0$ for $\gamma\geq 1$. The divergencies cancel out in the
gap equation in a way similar to cancellation of contributions from non-
magnetic impurities in a BCS superconductor. As the consequence, the gap
functions $\Delta_{n}(\omega_{m})$ and the condensation energies $E_{c,n}$
still form a discrete set for $\gamma\geq 1$, in which $E_{c,0}$ is the
largest. 111This does not hold for an extended $\gamma-$model with different
interaction strength in the particle-hole and the particle-particle channels.
In the latter case $E_{c,n}$ form a discrete set for $\gamma<1$, but the set
becomes dense at $\gamma\to 1$ and transforms into a continuous one for
$\gamma\geq 1$. This gives rise to a gapless branch of longitudinal gap
fluctuations.

In this paper we study the $\gamma$ model for $\gamma$ in the range
$1<\gamma<2$. Our goal is to understand whether the presence of an infinite
number of gap functions has physical consequences. The most natural conjecture
is that the existence of the set gives rise to a new branch of ”longitudinal”
gap fluctuations, what increases the strength of phase fluctuations and may
create a pseudogap, preformed pair region between the onset temperature for
the pairing and the actual superconducting $T_{c}$. However, this is by no
means guaranteed first of all because each $E_{c,n}$ is proportional to the
total number of particles in the system, and potential barriers between minima
with different $n$ are infinitely high. Strong fluctuating effects are only
possible if the spectrum of condensation energies gets dense and eventually
becomes a continuous one.

We argue that this happens when $\gamma$ is larger than one. We show that the
spectrum of $E_{c,n}$ gets denser when $\gamma$ increases and for $\gamma\leq
2$ splits into two sub-sets, one for $n<n_{cr}$ and another for $n>n_{cr}$,
where $n_{cr}\sim 1/(2-\gamma)$. Condensation energies for $n<n_{cr}$ get
progressively closer to $E_{c,0}$ as $\gamma$ approaches $2$, while for
$n>n_{cr}$, $E_{c,n}$ depends on the ratio of two large numbers $n$ and
$n_{cr}$, which is a near-continuous quantity. In the limit $\gamma\to 2$,
$n_{cr}$ tends to infinity, and $E_{c,n}$ for all finite $n$ collapse into a
single $E_{c,0}$, while for $n\to\infty$, the condensation energy becomes a
continuous function of $n/n_{cr}$, whose value is determined by how the double
limit $n\to\infty$ and $n_{cr}\to\infty$ (i.e. $\gamma\to 2$) is taken. This
creates a continuous gapless spectrum of longitudinal fluctuations.

We present the numerical evidence for the collapse of the condensation energy
obtained by analysing how the functional form of the eigenfunctions
$\Delta_{n}(\omega_{m})$ with different finite $n$ evolve with $\gamma$ We
argue that the frequency range, where $\Delta_{n}(\omega_{m})$ changes sign
$n$ times, shrinks, and moves to progressively smaller frequencies at
$\gamma\to 2$, consistent our assertion that in this limit $E_{c,n}$ with
finite $n$ tend to $E_{c,0}$.

The transformation of the spectrum of $\Delta_{n}(\omega_{m})$ and $E_{c,n}$
from a discrete to a continuous one is a rather non-trivial phenomenon, and we
study it from different corners. We find corroborative evidence by analyzing
the gap function on the Matsubara axis. We then analyze how the gap function
evolves on the real axis, and in the upper half-plane of frequency. For a
generic complex $z=\omega^{\prime}+i\omega^{{}^{\prime\prime}}$
($\omega^{{}^{\prime\prime}}>0$), $\Delta_{n}(z)$ is a complex function:
$\Delta_{n}(z)=\Delta^{\prime}_{n}(z)+i\Delta^{{}^{\prime\prime}}_{n}(z)$, and
can define its phase as $\eta_{n}(z)={\text{Im}}\log{\Delta_{n}(z)}$. Our
reasoning to focus on the real axis for $\gamma$ between one and two stems
from the fact that the interaction there, $V(\Omega)\propto
e^{i\pi\gamma/2}/|\omega|^{\gamma}$, is complex, and its real part,
$V^{{}^{\prime}}(\Omega)=\cos(\pi\gamma/2)/|\omega|^{\gamma}$, becomes
repulsive for $\gamma>1$, while the imaginary part
$V^{{}^{\prime\prime}}(\Omega)=\sin(\pi\gamma/2)/|\omega|^{\gamma}$ vanishes
at $\gamma\to 2$. Meanwhile, on the Matsubara axis, $V(\Omega_{m})\propto
1/|\Omega_{m}|^{\gamma}$ remains attractive.

We show that this gives rise to a dichotomy between the behavior of the gap
function on the Matsubara and the real axis at frequencies larger than
${\bar{g}}$. Namely, on the Matsubara axis, all $\Delta_{n}(\omega_{m})$
decrease as $1/|\omega_{m}|^{\gamma}$ ($n$ sign changes of
$\Delta_{n}(\omega_{m})$ occur at $\omega_{m}<{\bar{g}}$). On the real axis,
$\Delta_{n}(\omega)$ display this power-law behavior only at large enough
$\omega>\omega_{cr}$, where $\omega_{cr}\sim{\bar{g}}/(2-\gamma)^{1/2}$, while
at ${\bar{g}}<\omega<\omega_{cr}$, $\Delta_{n}(\omega)$ for all $n$ oscillate,
and the phase of the gap winds up by $2\pi m$, where $m$ is an integer. The
value of $m$ depends on $\gamma$, but not on $n$, and approaches infinity when
$\gamma\to 2$. We trace these oscillations to the existence of an
exponentially small and seemingly irrelevant oscillating term in
$\Delta_{n}(\omega_{m})$ on the Matsubara axis. We argue that this term
becomes the dominant one in the gap function upon analytical continuation to
the real axis. We compute the density of state (DoS) and show that for
$\gamma>1$ it develops a set of maxima and minima, and for $\gamma\to 2$ it
evolves towards a set of discrete $\delta-$functions, i.e., for $\gamma\to 2$,
the energy levels available to paired quasiparticles effectively get
quantized.

We next consider the gap function in the upper frequency half-plane. We show
that the phase winding by $2\pi m$ along the real axis between ${\bar{g}}$ and
$\omega_{cr}$ implies the presence of $m$ dynamical vortices in the same range
of $|z|$ in the upper frequency half-plane. We identify the vortices as the
points where $\Delta(z)=0$ and show that each vortex crosses into the upper
frequency half-plane from the lower one at some $\gamma_{i}>1$. The number of
vortices is determined by $\gamma$ and is the same for all $n$. At $\gamma\to
2$, the number of vortices increases, the location of vortex points $z_{i}$
becomes independent on $n$, and the array of vortices (points where
$\Delta_{n}(z_{i})=0$) extends to $|z|\to\infty$. Simultaneously, oscillations
of $\Delta_{n}(\omega)$ on the real axis also extend to an infinite frequency.
This implies that at $\gamma\to 2$, the value of $\Delta_{n}(|z|\to\infty)$
depends on the path, i.e., there is an essential singularity at
$|z|\to\infty$.

The presence of an essential singularity is crucial by the following reason.
Once the set of vortices at $z_{i}$ becomes infinite and the points $z_{i}$
accumulate at $|z|\to\infty$, one can use it as an input and obtain
$\Delta_{n}(\omega_{m})$ by analytically continuing from the set onto the
Matsubara axis. If an essential singularity was not there, the analytic
continuation would be unique and one would obtain $\Delta_{n}(z)=0$
everywhere. A finite $\Delta_{n}(z)$ emerges only because of essential
singularity, whose presence also makes the analytic continuation non-unique.
The outcome here is that the gap functions $\Delta_{n}$, which form a
continuum spectrum on the Matsubara axis and give rise to a continuum,
massless spectrum of the condensation energies, originate from an essential
singularity at $|z|=\infty$ on the complex frequency half-plane. At
$\gamma<2$, the number of vortices remains finite, and the spectrum of
$\Delta_{n}(\omega_{m})$ remains discrete.

The structure of the paper is the following. In Sec. II we briefly review the
model and the gap equation. In Sec. we clarify our assertion that the gap
functions $\Delta_{n}(\omega_{m})$ with different $n$ are topologically
distinct. We show that each nodal point of $\Delta_{n}(\omega_{m})$ is a
dynamical vortex with circulation $2\pi$, such that $\Delta_{n}(\omega_{m})$
is a gap function with $n$ dynamical vortices and is a dynamical analog of a
nodal topological superconductor. This holds for any $\gamma>0$. In Sec. III
we analyze the structure of $\Delta_{n}(\omega_{m})$ along the Matsubara axis
for $1<\gamma<2$. In Sec. III.2 we present the exact solution of the
linearized gap equation, $\Delta_{\infty}(\omega_{m})$ for $\gamma>1$ and show
that it has a small oscillating component, not present for smaller $\gamma$.
In Sec. III.3 we analyze the sign-preserving solution
$\Delta_{n=0}(\omega_{m})$. We show that at vanishing $\omega_{m}$,
$\Delta_{0}(\omega_{m})$ tends to a finite $\Delta_{0}(0)\sim{\bar{g}}$, and
at $\omega_{m}>{\bar{g}}$, it decreases as $1/|\omega_{m}|^{\gamma}$. In Sec.
III.4 we analyze $\Delta_{n}(\omega_{m})$ with a finite $n$ and the
corresponding condensation energy $E_{c,n}$. We show that at
$\omega_{m}<{\bar{g}}$, $\Delta_{n}(\omega_{m})$ oscillates $n$ times as
function of $\log{|\omega_{m}|}$ and then saturates at $\Delta_{n}(0)$, which
decrease exponentially with $n$, while at $\omega_{m}>{\bar{g}}$, it decreases
as $1/|\omega_{m}|^{\gamma}$ for all $n$. We show that as $\gamma$ increases
towards $2$, the spectrum of $E_{c,n}$ progressively splits into two subsets,
one for $n<n_{cr}\sim 1/(2-\gamma)$ and the other for $n>n_{cr}$. In Sec.
III.4.2 we discuss the evolution of $\Delta_{n}(\omega_{m})$ and the
corresponding onset temperatures $T_{p,n}$ at $\gamma\to 2$.

In Sec.IV we analyze the gap equation along real frequency axis. We present
the gap equation in Sec. IV.1 and solve it for $n=\infty$ in Sec. IV.2, for
$n=0$ in Sec. IV.3, and for a finite $n$ in Sec. IV.4. We show a complex
$\Delta_{n}(\omega)$ oscillates $n$ times at $\omega<{\bar{g}}$. This holds
for all $\gamma$. Beyond that, for $\gamma>1$, $\Delta_{n}(\omega)$ with any
$n$ oscillates at ${\bar{g}}<\omega<\omega_{cr}$, and the phase of the gap
function winds up by an integer number of $2\pi$ in this region.

In Sec. V we analyze the gap function $\Delta(z)$ in the upper frequency half-
plane, at $z=\omega^{\prime}+i\omega^{{}^{\prime\prime}}$. We again consider
separately the cases $n=\infty$ (Sec. V.1), $n=0$ (Sec. V.2) and a finite $n$
(Sec. V.3). We show that besides $n$ vortices on the Matsubara axis, which are
present for all $\gamma$, another set of vortices appears for $\gamma>1$.
These new dynamical vortices are located near the real axis, at
${\bar{g}}<|z|<\omega_{cr}$, and cause oscillations on the real axis in the
same range of $\omega$. Each vortex from this new set moves into the upper
half-plane of frequency from the lower one at a particular discrete
$\gamma_{i}$. The number of these vortices increases as $\gamma\to 2$ and they
line up along a particular path in the upper frequency half-plane.

We present our conclusions in Sec. VI. Some technical details of calculations
are presented in the Appendix. The case $\gamma=2$ requires special
consideration and we will analyze it in the next paper, where we also discuss
in detail the behavior of the superfluid stiffness and the interplay between
the pair formation and the true superconductivity.

##  II $\gamma$-model

We consider itinerant fermions at a QCP, interacting by exchanging
fluctuations of a critical order parameter. At a QCP, the propagator of a soft
boson becomes massless and mediates singular interaction between fermions. We
assume that this interaction is attractive in at least one pairing channel and
that a pairing boson can be treated as slow mode compared to a fermion, i.e.,
at a given momentum $q$, typical fermionic frequency is much larger than
typical bosonic frequency. In this situation, the ratio of typical bosonic and
fermionic frequencies is a small parameter, analogous to the Migdal parameter
in the electron-phonon case. This small parameter allows one to neglect vertex
corrections and explicitly integrate feromionic and bosonic propagators over
the momentum components. As a result, the problem reduces to a set of coupled
integral equations for frequency dependent fermionic self-energy and the
pairing vertex, in which an input is the effective, singular, frequency
dependent interaction
$V(\Omega_{m})={\bar{g}}^{\gamma}/|\Omega_{m}|^{\gamma}$. The exponent
$\gamma$ is different for different microscopic systems. This model, nicknamed
the $\gamma-$model, has been introduced in Paper I and we refer a reader to
the list of references in that paper to earlier works on the justification of
the model and its relation to various microscopic quantum-critical systems
with momentum and frequency-dependent interaction.

The local interaction $V(\Omega_{m})$ contributes to the fermionic self-energy
$\Sigma(k_{F},\omega_{m})=\Sigma(\omega_{m})$ and the pairing vertex
$\Phi(\omega_{m})$. The coupled equations for $\Sigma(\omega_{m})$ and
$\Phi(\omega_{m})$ have the same structure as Eliashberg equations for a
dispersion-less phonon, albeit with the exponent $\gamma$, instead of $2$, and
for shortness we will be calling them “Eliashberg equations”.

The superconducting gap function $\Delta(\omega_{m})$ is defined as

$\Delta(\omega_{m})=\omega_{m}\frac{\Phi(\omega_{m})}{{\tilde{\Sigma}}(\omega_{m})}=\frac{\Phi(\omega_{m})}{1+\Sigma(\omega_{m})/\omega_{m}}$
(1)

At a finite $T$ the Eliashberg gap equation is

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

The gap function is defined up to an overall $U(1)$ phase factor, which we set
to zero. The gap equation for infinitesimally small $\Delta(\omega_{m})$

$\Delta(\omega_{m})={\bar{g}}^{\gamma}\pi
T\sum_{m^{\prime}}\frac{\Delta(\omega^{\prime}_{m})-\Delta(\omega_{m})\frac{\omega^{\prime}_{m}}{\omega_{m}}}{|\omega^{\prime}_{m}||\omega_{m}-\omega^{\prime}_{m}|^{\gamma}}.$
(3)

determines the onset temperature for the pairing. As in previous papers, we
label this temperature $T_{p}$. Within Eliashberg theory, it is the same as
superconducting $T_{c}$, but in the presence of fluctuations $T_{p}$ is
generally larger than $T_{c}$.

##  III Gap function on the Matsubara axis, $T=0$

In Paper I we showed that at $T=0$ and $\gamma<1$, Eq. (2) has an infinite
number of topologically distant solutions, $\Delta_{n}(\omega_{m})$, all with
the same spatial gap symmetry. A gap function with a given $n$ tends to a
finite value at $\omega_{m}=0$ and decays as $1/|\omega_{m}|^{\gamma}$ at
$\omega_{m}\geq{\bar{g}}$, but in between it changes sign $n$ times. In Paper
III we showed that each of the two terms in the r.h.s. of (2) (separate
contributions from the pairing vertex and the fermionic self-energy) diverges
for $\gamma\geq 1$ as $\int d\Omega_{m}/|\Omega_{m}|^{\gamma}$. However, the
two divergencies cancel out, and $\Delta_{n}(\omega_{m})$ at $T=0$ evolve
smoothly through $\gamma=1$. We now extend the analysis to $\gamma$ between
$1$ and $2$. We show that the forms of $\Delta_{n}(\omega_{m})$ remain
qualitatively the same as for smaller $\gamma$, yet the spectrum of
condensation energies become progressively more dense. We also argue that when
$\gamma$ increases towards $2$, the spectrum of condensation energy $E_{n,c}$
becomes more and more dense, and in the limit $\gamma\to 2$ all $E_{n,c}$ with
finite $n<n_{cr}$ become almost equal to $E_{0,c}$, while $E_{n>n_{cr},c}$
form a continuous spectrum.

###  III.1 Nodal points on the Matsubara axis as dynamical vortices

Before we discuss $\Delta_{n}(\omega_{m})$, we pause for a moment and
elaborate on the notion that $\Delta_{n}(\omega_{m})$ are topologically
different. We argue that each nodal point is a dynamical vortex, i.e.,
$\Delta_{n}(\omega_{m})$ is an $n-$vortex state – a dynamical analog of a
nodal topological superconductor.

For definiteness, let’s compare the behavior of sign-preserving
$\Delta_{0}(\omega_{m})$ and of $\Delta_{1}(\omega_{m})$, and
$\Delta_{2}(\omega_{m})$, which change sign once and twice between
$\omega_{m}=0$ and $\omega_{m}=\infty$, respectively. We show these functions
in the left panel of Fig.1. Suppose that $\Delta_{1}(\omega_{m})$ changes sign
at $\omega_{m}=\omega_{1}$. Near this frequency,
$\Delta_{1}(\omega_{m})=-c(\omega_{m}-\omega_{1})$, where $c>0$ for
consistency with the figure. Let us analytically continue
$\Delta_{1}(\omega_{m})$ to a near vicinity of the Matsubara axis, to
$z=\omega^{\prime}+i\omega^{{}^{\prime\prime}}$ (on the Matsubara axis,
$z=i\omega_{m}$). Because $\Delta_{1}(\omega_{m})$ is non-singular,
$\Delta_{1}(z)=\Delta_{1}(i\omega_{m}\to\omega^{\prime}+i\omega^{{}^{\prime\prime}})$.
For any non-zero $\omega^{\prime}$, $\Delta_{1}(z)$ is a complex function:
$\Delta_{1}(z)=\Delta^{\prime}_{1}(z)+i\Delta^{{}^{\prime\prime}}_{1}(z)$, and
we can introduce the phase of $\Delta_{1}(z)$ as

$\eta_{1}(z)={\text{Im}}[\log{\Delta_{1}(z)}].$ (4)

We plot the variation of $\eta_{1}(z)$ around $z=i\omega_{1}$ in the right
panel of Fig.1.

Figure 1: Left panel: The schematic forms of the first 3 solutions of the gap
equation $\Delta_{0}(\omega_{m})$, $\Delta_{1}(\omega_{m})$ and
$\Delta_{2}(\omega_{m})$. A gap function $\Delta_{n}(\omega_{m})$ changes sign
$n$ times at $\omega_{m}>0$. Right panel: Variation of the phase of
$\Delta_{1}(z)$ ($z=\omega^{\prime}+i\omega^{{}^{\prime\prime}}$) around the
nodal point at $\omega_{m}=\omega_{m1}$ . Anticlockwise circulation of the
phase around $i\omega_{m1}$ is $2\pi$. The same holds for all other nodal
points.

We see that the phase varies by $2\pi$ upon anticlockwise circulation around
$\omega_{1}$. This implies that the nodal point at $\omega_{1}$ is in fact the
center of a dynamical vortex. One can verify that $\Delta_{n}(\omega_{m})$
with a generic $n$ contains $n$ vortices, each with anticlockwise circulation
$2\pi$.

Figure 2: Analytic continuation of $\Delta_{n}(\omega_{m})$ in Fig.1 to the
upper half plane ($i\omega_{m}\to
z=\omega^{\prime}+i\omega^{{}^{\prime\prime}}$) using Padé approximant. Panel
(a): the accuracy of fitting $\Delta_{n}(\omega_{m})$ by Padé approximant.
Panel (b): the phase of the gap function in the upper frequency half-plane,
$\eta_{n}(z)$ = Im$[\log{\Delta(z)}]$. The locations of the vortices are
marked by red dots. Panel (c): the gap function $\Delta_{n}(\omega)$ along the
real axis. Panel (d): variation of $\eta_{n}(\omega)$ along the real axis.
Each vortex gives rise to $2\pi$ phase variation. For convenience of
presentation, we confined $\eta_{n}(\omega)$ to $(-\pi,\pi)$. In this
convention, $\eta_{n}(\omega)$ winds by $2\pi n$ due to vortices.

One can straightforwardly verify that $n$ vortices on the Matsubara axis give
rise to $2\pi n$ phase variation on the real axis, between $-\infty$ and
$\infty$. To see this, one should compute
$\int_{-\infty}^{\infty}d\omega\partial\eta_{1}(\omega)/\partial\omega$ by
closing the integration contour over the upper half-plane. The function
$\partial\eta_{1}(\omega)/\partial\omega$ is analytic in the upper half-plane
except for the nodal points where it has simple poles. Modifying the contour
to circle out each nodal point, one obtains $2\pi n$ phase variation from the
vortices. In addition, there is also the $\pi\gamma$ from the integral over
the upper half-circle, due to the fact that at the largest frequencies,
$\Delta_{n}(z)\propto 1/(-z^{2})^{\gamma/2}$. This form is consistent with
$\Delta_{n}(\omega_{m})\propto 1/|\omega_{m}|^{\gamma}$, as one can verify by
Cauchy relation.

In Fig. 2 we show $\Delta_{n}(z)$, where
$z=\omega^{\prime}+i\omega^{{}^{\prime\prime}}$, for $n=0,1$, and $2$,
obtained by analytic continuation from $\Delta_{n}(\omega_{m})$ in Fig.1 using
Padé approximants methodBeach  _et al._ (2000). We see that
$\eta_{n}(\omega+i0^{+})$ indeed winds by $2\pi n$.

We now analyze the full frequency dependence of $\Delta_{n}(\omega_{m})$ for
$1<\gamma<2$. We first consider separately the opposite limits $n=\infty$ and
$n=0$ and then discuss a finite $n$.

###  III.2 Frequency dependence of $\Delta_{n}(\omega_{m})$. The case
$n=\infty$

The gap function $\Delta_{\infty}(\omega_{m})$ is a potential solution of the
linearized gap equation

$\Delta_{\infty}(\omega_{m})=\frac{{\bar{g}}^{\gamma}}{2}\int
d\omega^{\prime}_{m}\frac{\Delta_{\infty}(\omega^{\prime}_{m})-\Delta_{\infty}(\omega_{m})\frac{\omega^{\prime}_{m}}{\omega_{m}}}{|\omega^{\prime}_{m}|}~{}\frac{1}{|\omega_{m}-\omega^{\prime}_{m}|^{\gamma}}.$
(5)

We found that the solution exists and obtained the exact form of
$\Delta_{\infty}(\omega_{m})$. We emphasize that the solution exists despite
that $T=0$ is not a critical point for superconductivity. We will see below
that non-linear gap equation at $T=0$ has solution(s) with a finite gap
magnitude. The solution exists only at a QCP, when the pairing is mediated by
a gapless boson.

We present technical details of the analysis of Eq. (5) in Appendix A and here
quote the result. It is convenient to introduce

$y=\left(\frac{|\omega_{m}|}{\bar{g}}\right)^{\gamma},$ (6)

The gap function $\Delta(\omega_{m})=\Delta_{\infty}(y)$

$\Delta_{\infty}(y)=\int_{-\infty}^{\infty}dkb_{k}e^{-ik\log{y}}.$ (7)

is the Fourier transform of $b_{k}$ given by

$b_{k}=\frac{e^{-i(I_{k}+k\log{(\gamma-1)})}}{\left[\cosh(\pi(k-\beta))\cosh(\pi(k+\beta))\right]^{1/2}}.$
(8)

Here

$I_{k}=\frac{1}{2}\int_{-\infty}^{\infty}dk^{\prime}\log{|\epsilon_{k^{\prime}}-1|}\tanh{\pi(k^{\prime}-k)},$
(9)
$\epsilon_{k^{\prime}}=\frac{1-\gamma}{2}\frac{\Gamma\left(\frac{\gamma}{2}\left(1+2ik^{\prime}\right)\right)\Gamma\left(\frac{\gamma}{2}\left(1-2ik^{\prime}\right)\right)}{\Gamma(\gamma)}\left(1+\frac{\cosh{\pi\gamma
k^{\prime}}}{\cos{\pi\gamma/2}}\right),$ (10)

and $\beta>0$ is the solution of $\epsilon_{\beta}=1$. We plot
$\epsilon_{k^{\prime}}$ in Fig. 3 (a). The value of $\beta$ depends on
$\gamma$ and evolves between $\beta\approx 0.79$ for $\gamma=1$ and
$\beta\approx 0.39$ for $\gamma=2$.

Figure 3:  Functions $\epsilon_{k}$ (a) and $\epsilon_{ik}$ (b) for
$\gamma=1.5$. The red dots show the solutions $\epsilon_{\beta}=1$ and
$\epsilon_{i\beta_{m}}=1$ ($m=0,1,2,...$).

The gap function $\Delta_{\infty}(y)$ can be straightforwardly obtained
numerically. We show the results in Fig. 4 for $\gamma=1.91$. At small $y$,
i.e., at $\omega_{m}<{\bar{g}}$, the gap function oscillates down to the
lowest frequencies with the period set by $\log{y}$. At $y>1$, i.e., at
$\omega>{\bar{g}}$, oscillations end, and $\Delta_{\infty}(y)$ gradually
decreases with increasing $y$. By practical reasons, we set the lowest
frequency in Fig. 4 at $10^{-18}{\bar{g}}$. There are 10 zeros of
$\Delta_{\infty}(\omega)$ (10 vortex points) above this frequency.

The limiting forms of $\Delta_{\infty}(y)$ at small and large $y$ can be
obtained analytically. At $y\ll 1$ we find

$\Delta_{\infty}(y)\sim y^{1/2}\cos{\left(\beta\log{y}+\phi\right)},$ (11)

where $\phi$ is a $\gamma-$dependent number. At $y\gg 1$

$\Delta_{\infty}(y)\sim\frac{1}{y}$ (12)

Figure 4: The exact solution $\Delta_{\infty}(\omega_{m})$ for $\gamma=1.91$.
To better display the distinctive power-law behavior, we use different scales
for $\omega>{\bar{g}}$ and $\omega<{\bar{g}}$.

The corrections to (11) hold in powers of $y$ and form two series, local and
non-local, depending on whether they come from fermions with running
frequencies of order $\omega_{m}$ or much larger frequencies, of order
${\bar{g}}$. The full $\Delta_{\infty}(y)$ for $y<1$ is

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

where the first term describes the sum of the leading term and the local
corrections, and the second term describes the non-local corrections. In the
latter, $\beta_{m}>0$ are the solutions of $\epsilon_{i\beta_{m}}=1$ for
imaginary argument(same as in Fig.3). There is an infinite set of such
$\beta_{m}$, specified by integer $m=0,1,2..$ and located at
$1/2+2m/\gamma<\beta_{m}<1/2+2(m+1)/\gamma$, see Fig. 3 (b).

For $y>1$, the corrections to (12) hold in powers of $1/y$ and also form two
series:

$\Delta_{\infty}(y)=\sum_{n=0}^{\infty}C^{>}_{n}~{}\left(\frac{1}{y}\right)^{(n+1)}+\sum_{n,m=0}^{\infty}D^{>}_{n,m}\left(\frac{1}{y}\right)^{(n+1+2(m+1)/\gamma)}$
(14)

Eq. (12) is the $C^{>}_{0}$ term in this series.

####  III.2.1 Distinction between $\Delta_{\infty}(y)$ at $\gamma<1$ and
$\gamma>1$.

The limiting forms of $\Delta_{\infty}(y)$ at small and large $y$ and the
general structure of $\Delta_{\infty}(\omega_{m})$ for $\gamma>1$ are the same
as we obtained earlier for $\gamma<1$. However, on a more careful look, we
found that there is a qualitative distinction between the two cases for $y>1$.
Namely, for $\gamma<1$, the series in (14) are convergent, while for
$\gamma>1$, series expansion in (14) holds only up to a certain $n_{max}\sim
y^{1/(\gamma-1)}$, while instead of power series with larger $n$, there
appears an oscillating term (see below). The value $n_{max}$ is large for
$y\gg 1$, such that this term has an exponentially small prefactor
$e^{-y^{1/(\gamma-1)}}$. Yet, we will see below that once $\Delta_{\infty}$ is
analytically continued from Matsubara to real axis, the exponential smallness
is reduced, particularly near $\gamma=2$, and the oscillating term becomes the
dominant contribution to $\Delta_{\infty}(\omega)$ over a wide frequency range
and gives rise to the emergence of a new vortex structure in the upper half-
plane of frequency.

The difference between $\Delta_{\infty}(y\gg 1)$ for $\gamma<1$ and $\gamma>1$
originates from change of the behavior of $\epsilon_{k^{\prime}}$ in Eq. (9).
At $\gamma<1$, $\epsilon_{k^{\prime}}$ vanishes at $k^{\prime}\to\infty$,
while at $\gamma>1$, $\epsilon_{k^{\prime}}$ diverges at large $k^{\prime}$ as

$\epsilon_{k^{\prime}}=(A_{\gamma}|k^{\prime}|)^{\gamma-1}(\gamma-1),~{}~{}A_{\gamma}=\gamma\left(\frac{\pi}{2\Gamma(\gamma)\cos\left[\pi(2-\gamma)/2\right]}\right)^{\frac{1}{\gamma-1}}$
(15)

Substituting this form into (9) and integrating over $k^{\prime}$, we find
that at large $k$,

$I_{k}+k\log{(\gamma-1)=-(\gamma-1)k\log{\frac{A_{\gamma}k}{e}}}$ (16)

The integral comes from $k^{\prime}\sim k$, what justifies using the form of
$\epsilon_{k^{\prime}}$ at large $k^{\prime}$. Substituting (16) into (8) and
approximating the denominator in (8) by its value at large $k$, we obtain

$b_{k}\approx e^{-\pi|k|}e^{i(\gamma-1)k\log{\frac{A_{\gamma}k}{e}}}$ (17)

Substituting this $b_{k}$ into (7), we find the contribution to
$\Delta_{\infty}(y)$ in the form

$\Delta_{\infty}(y)\sim\int_{k_{min}}^{\infty}e^{-\pi
k}\cos{\left((\gamma-1)k\log{\frac{y^{1/(\gamma-1)}e}{A_{\gamma}k}}\right)}$
(18)

where $k_{min}=O(1)$. The argument of the $\cos$ has a maximum at
$k=k^{*}=y^{1/(\gamma-1)}/A_{\gamma}$. For large $y$, this $k^{*}$ is well
inside the range of the integration. Expanding to quadratic order near the
maximum and substituting into (18), we find after some algebra (see Appendix
A) that the integral in (18) contains the universal contribution, which we
label as $\Delta^{u}_{\infty}$:

$\Delta^{u}_{\infty}(y)\sim
y^{\frac{1}{2(\gamma-1)}}e^{-\frac{\pi}{A_{\gamma}}y^{\frac{1}{(\gamma-1)}}}\cos{\left[y^{\frac{1}{(\gamma-1)}}~{}\frac{(\gamma-1)}{A_{\gamma}}\left(1-\frac{\pi^{2}}{2(\gamma-1)^{2}}\right)-\frac{\pi}{4}\right]}$
(19)

In the actual Matsubara frequencies this formula reads

$\Delta^{u}_{\infty}(\omega_{m})\propto|\omega_{m}|^{\frac{\gamma}{2(\gamma-1)}}e^{-\frac{\pi(|\omega_{m}|/{\bar{g}})^{\gamma/(\gamma-1)}}{A_{\gamma}}}\cos{\left[\left(\frac{|\omega_{m}|}{{\bar{g}}}\right)^{\gamma/(\gamma-1)}\frac{(\gamma-1)}{A_{\gamma}}\left(1-\frac{\pi^{2}}{2(\gamma-1)^{2}}\right)-\frac{\pi}{4}\right]}$
(20)

The oscillating contribution is exponentially small at large $y$. It vanishes
at $\gamma=1$ and does not exist at smaller $\gamma$. We will return to this
contribution in Sec. IV.2, where we discuss $\Delta^{u}_{\infty}(\omega)$ on
the real axis.

###  III.3 Frequency dependence of $\Delta_{n}(\omega_{m})$. The case $n=0$

We now consider the solutions of the non-linear gap equation (2) at $T=0$. We
begin with sign-preserving solution $\Delta_{0}(\omega)$. To set the stage, we
first quickly review the results for $\gamma<1$. For such $\gamma$, the term
with $m^{\prime}=m$ in the r.h.s. of (2) is non-divergent and can be neglected
as it renormalizes the prefactor for $\Delta_{0}(\omega_{m})$ in the r.h.s. by
at most $O(1)$. Without this term, we have

$\Delta_{0}(\omega_{m})=\frac{{\bar{g}}^{\gamma}}{2}\int_{-\infty}^{\infty}\frac{d\omega^{\prime}_{m}\Delta_{0}(\omega^{\prime}_{m})}{|\omega_{m}-\omega^{\prime}_{m}|^{\gamma}\sqrt{\Delta^{2}_{0}(\omega^{\prime}_{m})+(\omega^{\prime}_{m})^{2}}},$
(21)

At small $\omega_{m}$, $\Delta_{0}(\omega_{m})\approx\Delta_{0}(0)$, and we
have

$\Delta_{0}(0)={\bar{g}}^{\gamma}\int_{0}^{\infty}d\omega^{\prime}_{m}\frac{\Delta_{0}(\omega^{\prime}_{m})}{|\omega^{\prime}_{m}|^{\gamma}\sqrt{\Delta^{2}_{0}(\omega^{\prime}_{m})+(\omega^{\prime}_{m})^{2}}}$
(22)

The integral is determined by
$\omega^{\prime}_{m}\sim\Delta_{0}(\omega^{\prime}_{m})$. For such
$\omega^{\prime}_{m}$, we can approximate $\Delta_{0}(\omega^{\prime}_{m})$ by
$\Delta_{0}(0)$ up to $O(1)$ corrections. Evaluating the integral in (22), we
then obtain $\Delta_{0}(0)\sim{\bar{g}}$. In the opposite limit
$\omega_{m}\gg{\bar{g}}$, we can pull out $1/|\omega|^{\gamma}$ from the
r.h.s. of (21) and obtain

$\Delta_{0}(\omega_{m})=\left(\frac{{\bar{g}}}{|\omega_{m}|}\right)^{\gamma}\int_{0}^{\infty}d\omega^{\prime}_{m}\frac{\Delta_{0}(\omega^{\prime}_{m})}{\sqrt{\Delta^{2}_{0}(\omega^{\prime}_{m})+(\omega^{\prime}_{m})^{2}}}$
(23)

The integral is determined by $\omega^{\prime}_{m}=O({\bar{g}})$ and is
$O(\bar{g})$, such that at high frequencies,
$\Delta_{0}(\omega_{m})\sim{\bar{g}}({\bar{g}}/|\omega_{m}|)^{\gamma}$. The
low-frequency and high-frequency forms of $\Delta_{0}(\omega_{m})$ then match
at $\omega_{m}\sim{\bar{g}}$.

Now we do the same analysis for $\gamma>1$. For $\gamma>1$, there is an
identity

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

We use it to re-express Eq. (2) as

$\displaystyle\Delta_{0}(\omega_{m})\left[1-\frac{{\bar{g}}^{\gamma}}{2}\int_{-\infty}^{\infty}d\omega^{\prime}_{m}\frac{1-\frac{\omega^{\prime}_{m}}{\omega_{m}}}{|\omega_{m}-\omega^{\prime}_{m}|^{\gamma}}\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}}^{\gamma}}{2}\int_{-\infty}^{\infty}d\omega^{\prime}_{m}\frac{\Delta_{0}(\omega^{\prime}_{m})-\Delta_{0}(\omega_{m})}{|\omega_{m}-\omega^{\prime}_{m}|^{\gamma}\sqrt{\Delta^{2}_{0}(\omega^{\prime}_{m})+(\omega^{\prime}_{m})^{2}}}$
(25)

Each integral in (25) is infra-red convergent and determined by
$\omega^{\prime}_{m}\sim\Delta_{0}(\omega^{\prime}_{m})$.

In the limit $\omega_{m}\to 0$, Eq. (25) reduces to

$\displaystyle\Delta_{0}(0)\left[1-{\bar{g}}^{\gamma}(\gamma-1)\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}|^{\gamma}}\right]$
$\displaystyle={\bar{g}}^{\gamma}\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}|^{\gamma}}$
(26)

To estimate the magnitude of $\Delta_{0}(0)$, assume that
$\Delta_{0}(\omega^{\prime}_{m})$ remains approximately equal to
$\Delta_{0}(0)$ for relevant
$\omega^{\prime}_{m}\leq\Delta_{0}(\omega^{\prime}_{m})$.

The r.h.s. of (26) then vanishes, and we obtain

$1={\bar{g}}^{\gamma}(\gamma-1)\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}|^{\gamma}}=\frac{1-\gamma}{2\sqrt{\pi}}\left(\frac{\bar{g}}{\Delta_{0}(0)}\right)^{\gamma}\Gamma\left(\frac{1-\gamma}{2}\right)\Gamma\left(\frac{\gamma}{2}\right)$
(27)

We plot $\Delta_{0}(0)/\bar{g}$ as a function of $\gamma$ in Fig.5 (a).

Figure 5: Panel (a): the solution of Eq.(27) for $\Delta_{0}(\omega_{m}=0)$;
Panel (b): the onset temperature of the pairing $T_{p,0}$ vs $\gamma$. Panel
(c): the ratio $2\Delta_{0}(0)/T_{p,0}$ vs $\gamma$. The ratio is slightly
different from the one obtained numerically in Ref. Wu  _et al._ (2019a) as
here we use Eq.(27) for $\Delta_{0}(\omega_{m}=0)$.

We see that $\Delta_{0}(0)$ is of order ${\bar{g}}$ and weakly depends on
$\gamma$. One can further expand $\Delta_{0}(\omega_{m})$ in frequency and
find that the expansion is analytic and holds in powers of $\omega^{2}_{m}$.
For completeness, we also show the onset temperature $T_{p,0}$ and the ratio
$2\Delta_{0}(0)/T_{p,0}$, both as functions of $\gamma$. The results are
consistent with earlier calculations Marsiglio and Carbotte (1991); Wang  _et
al._ (2016); Lee  _et al._ (2018); Wu  _et al._ (2019a).

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

$\Delta_{0}(\omega_{m})\approx
Q_{\gamma}\left(\frac{{\bar{g}}}{|\omega_{m}|}\right)^{\gamma}$ (28)

where

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

The integral is again determined by frequencies at which
$\Delta_{0}(\omega^{\prime}_{m})\sim\omega^{\prime}_{m}$ and is of order
$\Delta_{0}(0)$. Then
$\Delta_{0}(\omega_{m})\sim\Delta_{0}(0)({\bar{g}}/|\omega_{m}|)^{\gamma}$. In
Fig.6 we show $Q_{\gamma}$, obtained by solving numerically the nonlinear gap
equation for various $\gamma$. We see that $Q_{\gamma}$ is indeed of order
$Q_{\gamma}\sim\Delta_{0}(0)\sim\bar{g}$.

Figure 6: $Q_{\gamma}$, defined in Eq. (29), as a function of $\gamma$. To
obtain this function, we solved numerically the non-linear gap equation for
sign-preserving $\Delta_{0}(\omega_{m})$ at a small temperature.

The high-frequency behavior of $\Delta_{0}(\omega_{m})$ in (28) is the same as
of $\Delta_{\infty}(\omega_{m})$. This is expected because at such frequencies
$\Delta_{0}(\omega_{m})\ll\omega_{m}$. This implies that
$\Delta_{0}(\omega_{m})$ has the same exponentially small oscillating term as
in Eq. (19). In Sec.IV.2 we will see how this term evolves upon analytical
continuation to the real axis.

###  III.4 Frequency dependence of $\Delta_{n}(\omega_{m})$. Finite $n$.

####  III.4.1 A generic $\gamma$ from $1<\gamma<2$.

To analyze the behavior of $\Delta_{n}(\omega_{m})$ at a finite $n$, one has
to solve the non-linear gap equation at $T=0$ for a class of functions which
change sign $n$ times. At $\omega_{m}>{\bar{g}}$, $\Delta_{n}(\omega_{m})$
decreases as $1/|\omega_{m}|^{\gamma}$ and has an exponentially small
oscillating piece, the same as in (19). At $\omega_{m}<{\bar{g}}$,
$\Delta_{n}(\omega_{m})$ oscillates $n$ times and saturates at the smallest
$\omega_{m}$ at $\Delta_{n}(0)$, which decreases with increasing $n$. To
estimate the magnitude of $\Delta_{n}(0)$ and, more generally, understand why
there is a discrete set of solutions, we use the solution of the linearized
gap equation, Eq. (11), as an input, and treat the non-linear terms by
expanding in powers of $\Delta_{n}(\omega_{m})/\omega_{m}$. The gap function,
obtained this way, is represented by the series

$\Delta_{n}(y)={\bar{g}}y^{1/2}\sum_{m=0}^{\infty}C_{n}^{2m+1}f_{m}(y),~{}~{}y=\left(\frac{|\omega_{m}|}{\bar{g}}\right)^{\gamma}$
(30)

where $f_{0}(y)=\cos{(\beta\log{y}+\phi)}$, $f_{m\geq 1}(y)$ are obtained in
the order-by-order expansion in $\Delta_{n}(\omega_{m})/\omega_{m}$, and
$C_{n}$ is yet unknown factor. Evaluating the few first $f_{m}(y)$, we obtain
that they are also oscillating functions of $\log{y}$ with $y-$dependent
prefactors of order $(1/y^{2m})^{1/\gamma-1/2}$. The perturbative expansion
holds as long as $C_{n}<y^{(1/\gamma-1/2)}$. In this range, one can
approximate the r.h.s. of (30) by the $m=0$ term for order-by-magnitude
estimates. The function $f_{0}(y)$ changes sign $n$ times between $y=O(1)$ and
$y_{min}\sim e^{-n\pi/\beta}$. It is natural to associate $y_{min}$ with the
lower boundary of the perturbative expansion. Doing this, we obtain a discrete
set of $C_{n}$:

$C_{n}\sim e^{\frac{-n\pi}{\beta}\left(1/\gamma-1/2\right)}$ (31)

For large $n$,
$\Delta_{n}(y)\approx{\bar{g}}y^{1/2}C_{n}\cos{(\beta\log{y}+\phi)}$. For
$n=O(1)$, $C_{n}=O(1)$, and one needs to keep the full series in (30). This
reasoning also yields $\Delta_{n}(0)\sim
y^{1/2}_{min}C_{n}\sim(y_{min})^{1/\gamma}\sim{\bar{g}}e^{-n\pi/(\beta\gamma)}$.

####  III.4.2 $\Delta_{n}(\omega_{m})$ for $\gamma$ close to $2$.

Eq. (31) shows that for a generic $\gamma$, $C_{n}=O(1)$ for $n=O(1)$, and
decays exponentially at larger $n$. However, for $\gamma$ close to $2$, the
dependence on $n$ in (31) is actually via the product $n^{*}=n(2-\gamma)\sim
n/n_{cr}$. In the limit $\gamma\to 2$, $n^{*}\to 0$ for any finite $n$. The
corresponding $C_{n}$ then become independent on $n$ and coincide with $C_{0}$
up to corrections, which become relevant only at the smallest $y$. At the same
time, at $n\to\infty$, $n^{*}$ becomes a continuous variable, whose value
depends on how the double limit $n\to\infty$ and $\gamma\to 2$ is taken. Then
$C_{n\to\infty}$ and $\Delta_{n\to\infty}$ become continuous functions of
$n^{*}$. 222This is similar to the behavior of $\Delta_{n}$ for $\gamma\approx
1$ in the extended $\gamma$ model, when the extension parameter $N\neq 1$. We
illustrate this in schematically in Fig.7.

Figure 7: A schematic dependence of $\Delta_{n}$ on $n$ at $\gamma\to 2$, at
some fixed $\omega_{m}$. All $\Delta_{n}$ with finite $n$ tend to the same
value, while the magnitude of $\Delta_{n\to\infty}$ depends on the ratio
$n^{*}=n/n_{cr}$, where $n_{cr}\sim 1/(2-\gamma)$.

We show the numerical evidence for this behavior in Figs.8 and 9. We use the
fact that each $\Delta_{n}(\omega_{m})$ terminates at its own $T_{p,n}$ and
analyze the functional forms of the eigenvalues $\Delta_{n}(\omega_{m})$ at
$T_{p,n}$, assuming that these forms do not change significantly between
$T_{p,n}$ and $T=0$. In Fig.8 we plot the largest frequency at which
$\Delta_{1,2}(\omega_{m})$ change sign We see that these frequencies get
progressively smaller as $\gamma$ approaches $2$, while at larger $\omega_{m}$
the functional forms of $\Delta_{1}(\omega_{m})$ and $\Delta_{2}(\omega_{m})$
are the same as of $\Delta_{0}(\omega_{m})$. We expect the same to hold for
other finite $n$.

Figure 8: The positions of the second zero of $\Delta_{2}(\omega_{m})$, as a
functions of $\gamma$. This frequency decreases as $\gamma\to 2$, implying
that sign change at the highest $\omega_{m}$ occurs at progressively smaller
frequency. The gap function $\Delta_{2}(\omega_{m})$ has been obtained by
solving the linearized gap equation at the corresponding $T_{p,2}$. We
conjecture that the functional form of $\Delta_{2}(\omega_{m})$ does not
change substantially between $T_{p,2}$ and $T=0$.

On a more careful look, we find that the shift of oscillations to smaller
frequencies is related to the reduction of the corresponding $T_{p,n}$.
Namely, we verified that the largest frequency, at which
$\Delta_{n}(\omega_{m})$ changes sign, scales with the corresponding
$T_{p,n}$. In Fig. 9 we show $T_{p,n}$ for $n=1,2$ as a function of $\gamma$.
We see that $T_{p,n}$ decreases as $\gamma$ increases toward $2$ and, within
our accuracy, vanishes at $\gamma=2$.

Now, each gap function $\Delta_{n}(\omega_{m})$ generates a certain
condensation energy $E_{c,n}$. For a generic $\gamma<2$, the spectrum of
$E_{c,n}$ is discrete, and $E_{c,0}$ is the largest by magnitude. Then, at low
$T\ll T_{p,0}$, only the $n=0$ solution matters, while the existence of other
$E_{c,n}$ affects the system behavior only at $T\leq T_{p,0}$. As $\gamma$
increases towards two, the spectrum of $E_{c,n}$ becomes denser, and
$(E_{c,0}-E_{c,n})/E_{c,0}$ progressively gets smaller for any finite $n$. At
$\gamma\to 2$, the spectrum of $E_{c,n}$ can be viewed as almost continuous
spectrum with a small gap, i.e., there emerges a branch of low-energy
”longitudinal” gap fluctuations. These fluctuations affect the system behavior
beginning at a progressively smaller $T$, as $\gamma$ approaches two. At
$\gamma=2-0$, the spectrum of condensation energies becomes a gapless
continuous function $E_{c}(n^{*})$ with $E_{c}(0)=E_{c,0}$ and
$E_{c}(\infty)=E_{c,\infty}$. In the next Paper V, where we specifically
analyze the case $\gamma=2$, we show there that gapless longitudinal
fluctuations give rise to singular downward renormalization of the stiffness.

Figure 9: The onset temperature $T_{p,1}$(left) and $T_{p,2}$(right) vs
$\gamma$. Within our accuracy, $T_{p,1}$ and $T_{p,2}<T_{p,1}$ vanish at
$\gamma=2$. Other $T_{p,n}$ with $n>2$ are smaller than $T_{p,2}$ and also
vanish at $\gamma\to 2$.

##  IV Gap function on the real axis

The analysis on the Matsubara axis is sufficient for the computation of the
condensation energy at $T=0$ and thermodynamic properties at $T>0$. For the
analysis of transport and spectroscopic properties, one needs to know the gap
function on the real frequency axis. In this section we obtain
$\Delta_{n}(\omega)$ for real $\omega$. The gap function on the real axis is
complex:
$\Delta_{n}(\omega)=\Delta^{\prime}_{n}(\omega)+i\Delta^{{}^{\prime\prime}}(\omega)$.
It is convenient to introduce its $\omega-$dependent phase as
$\Delta_{n}(\omega)=|\Delta_{n}(\omega)|e^{i\eta(\omega)}$. As usual,
$\eta(\omega)$ is defined up to a constant.

###  IV.1 Eliashberg gap equation on the real axis

Eliashberg gap equation in real frequencies cannot be obtained by simply
replacing $\omega_{m}$ by $-i\omega$ in (2) as
$\Delta_{n}(\omega_{m}\to-i\omega)$ would have branch cuts in the upper
frequency half-plane, while the actual $\Delta(z)$ must be an analytic
function.

Below we follow the approach, suggested in Refs Karakozov  _et al._ (1991);
Marsiglio and Carbotte (1991); Combescot (1995); Wu  _et al._ (2019b) for
electron-phonon problem, i.e., convert to real axis using the spectral
representation $f(i\omega_{m})=(1/\pi)\int
dxf^{{}^{\prime\prime}}(x)/(x-i\omega_{m})$ and, where possible, use
$\Delta(\omega_{m})$, which we treat as the known function. For our $\gamma$
model this procedure yields

$\Delta(\omega)B(\omega)=A(\omega)+C(\omega)$ (32)

where, as before, $D(\omega)=\Delta(\omega)/\omega$, and

$\displaystyle A(\omega)$ $\displaystyle=\pi
T\sum_{\omega_{m}>0}\frac{D(\omega_{m})}{\sqrt{1+D^{2}(\omega_{m})}}\left(V(\omega_{m}+i\omega)+V(\omega_{m}-i\omega)\right)$
(33) $\displaystyle B(\omega)$
$\displaystyle=1+\frac{i\pi}{\omega}T\sum_{\omega_{m}>0}\frac{1}{\sqrt{1+D^{2}(\omega_{m})}}\left(V(\omega_{m}+i\omega)-V(\omega_{m}-i\omega)\right)$
$\displaystyle C(\omega)$
$\displaystyle=\frac{i}{2}\int_{-\infty}^{\infty}d\Omega
V^{{}^{\prime\prime}}(\Omega)\frac{D(\omega-\Omega)-D(\omega)}{\sqrt{1-D^{2}(\omega-\Omega)}}\left(\tanh\frac{\omega-\Omega}{2T}+\coth\frac{\Omega}{2T}\right)$

The functions $A(\omega)$ and $B(\omega)$ are obtained by just rotating the
gap equation from $\omega_{m}$ to $-i\omega$, and the function $C(\omega)$
contains the extra term that cancels out a parasitic contribution from the
branch cuts.

The interaction $V(\Omega)$ on the real axis is
$V(\omega)=V^{\prime}(\Omega)+iV^{{}^{\prime\prime}}(\Omega)$, where

$V^{\prime}(\Omega)=\left(\frac{{\bar{g}}}{|\Omega|}\right)^{\gamma}\cos{\frac{\pi\gamma}{2}},~{}~{}V^{{}^{\prime\prime}}(\Omega)=\left(\frac{{\bar{g}}}{|\Omega|}\right)^{\gamma}\sin{\frac{\pi\gamma}{2}}{\text{sign}}\omega$
(34)

Observe that on the Matsubara axis,
$V(\Omega_{m})=({\bar{g}}/|\Omega_{m}|)^{\gamma}$ positive (attractive in our
notations), on a real axis $V^{{}^{\prime}}(\omega)$ is negative (repulsive)
for $\gamma>1$.

Below we focus on the case of zero temperature. At $T=0$, the expressions for
$A,B$ and $C$ reduce to

$A(\omega)=\bar{g}^{\gamma}\frac{2}{\pi}\sin\frac{\pi\gamma}{2}\int_{0}^{\infty}d\omega_{m}\frac{D(\omega_{m})}{\sqrt{1+D^{2}(\omega_{m})}}\int_{0}^{\infty}\frac{d\Omega}{\Omega^{\gamma-1}}\frac{\Omega^{2}-\omega^{2}+\omega^{2}_{m}}{[(\Omega-\omega)^{2}+\omega^{2}_{m}][(\Omega+\omega)^{2}+\omega^{2}_{m}]}$
(35)
$B(\omega)=\left(1+\bar{g}^{\gamma}\frac{4}{\pi}\sin\frac{\pi\gamma}{2}\int_{0}^{\infty}\frac{d\omega_{m}}{\sqrt{1+D^{2}(\omega_{m})}}\int_{0}^{\infty}\frac{d\Omega}{\Omega^{\gamma-1}}\frac{\omega_{m}}{[(\Omega-\omega)^{2}+\omega^{2}_{m}][(\Omega+\omega)^{2}+\omega^{2}_{m}]}\right)$
(36)
$C(\omega)=i{\bar{g}}^{\gamma}~{}\sin{\frac{\pi\gamma}{2}}~{}{\text{sign}}\omega\int_{0}^{|\omega|}\frac{d\Omega}{\Omega^{\gamma}}~{}~{}\frac{D(|\omega|-\Omega)-D(|\omega|)}{\sqrt{1-D^{2}(|\omega|-\Omega)}}$
(37)

At large $\omega$,

$B(\omega)\approx 1+\left(\frac{\bar{g}}{|\omega|}\right)^{\gamma}I_{\gamma},$
(38)

where

$I_{\gamma}=\frac{i}{2}\int_{0}^{\infty}\frac{(x-i)^{\gamma}-(x+i)^{\gamma}}{(x^{2}+1)^{\gamma}}$
(39)

In the two limits, $I_{1}=\pi/2,I_{2}=1$.

###  IV.2 Frequency dependence of $\Delta_{n}(\omega)$. The case $n=\infty$.

We verified that the exact $\Delta_{\infty}(\omega)$ can be obtained by
converting Eq. (7) from Matsubara to real axis by a rotation, i.e., by
replacing $\omega_{m}$ by $-i\omega+\delta$. Under such rotation,
$e^{-ik\log{y}}=e^{-ik\gamma\log{|\omega_{m}|/{\bar{g}}}}$ transforms into
$e^{-ik\gamma\log{|\omega|/{\bar{g}}}}e^{-(k\pi\gamma/2)\text{ sgn}\omega}$.
At small $\omega<{\bar{g}}$, $\Delta^{\prime}_{\infty}(\omega)$ and
$\Delta^{{}^{\prime\prime}}_{\infty}(\omega)$ oscillate as functions of
$\log|\omega|$ down to the lowest frequencies. In explicit form, we have in
this regime

$\Delta_{\infty}(\omega)=C_{\infty}\left(\frac{|\omega|}{\bar{g}}\right)^{\gamma/2}e^{-\frac{i\pi\gamma}{4}\text{sgn}\omega}\cos{\left(\beta\gamma\left(\log{\frac{|\omega|}{\bar{g}}}-i\frac{\pi}{2}\text{sgn}\omega\right)+\phi\right)}$
(40)

where $C_{\infty}$ is infinitesimally small. Separating real and imaginary
parts, we obtain

$\displaystyle\Delta^{\prime}_{\infty}(\omega)=C_{\infty}\left(\cos{\frac{\pi\gamma}{4}}\cos{D_{1}}\cosh{D_{2}}+\sin{\frac{\pi\gamma}{4}}\sin{D_{1}}\sinh{D_{2}}{\text{sign}}\omega\right)$
$\displaystyle\Delta^{{}^{\prime\prime}}_{\infty}(\omega)=C_{\infty}\left(\cos{\frac{\pi\gamma}{4}}\sin{D_{1}}\sinh{D_{2}}-\sin{\frac{\pi\gamma}{4}}\cos{D_{1}}\cosh{D_{2}}{\text{sign}}\omega\right)$
(41)

where

$\displaystyle
D_{1}=\beta\log{\left(\frac{|\omega|}{\bar{g}}\right)^{\gamma}}+\phi,$
$\displaystyle D_{2}=\frac{\pi\beta\gamma}{2}{\text{sign}}\omega$ (42)

We show $\Delta^{\prime}_{\infty}(\omega)$,
$\Delta^{{}^{\prime\prime}}_{\infty}(\omega)$, and the phase
$\eta_{\infty}(\omega)$ for $\omega<{\bar{g}}$ in Fig. 10 for $\gamma=1.91$.
We set the lowest frequency in at $10^{-17}{\bar{g}}$. In this case the total
phase variation up to $O({\bar{g}})$ is $10\pi$. This is consistent with 10
vortices on the Matsubara axis (Fig. 20) as each vortex gives rise to a phase
change of $\pi$ for positive $\omega$ on the real axis and another phase
change of $\pi$ for negative $\omega$.

Figure 10: The complex gap function $\Delta_{\infty}(\omega)$ along the real
frequency axis and its phase $\eta_{\infty}(\omega)$, constrained to
$(-\pi,\pi)$, for $\gamma=1.91$. $\Delta_{\infty}(\omega)$ oscillates an
infinite number of times between $\omega=0$ and $\omega/{\bar{g}}=O(1)$. For
better understanding of the total phase variation $\delta\eta$ , we count
oscillations starting from a small but finite $\omega_{min}\sim
10^{-17}\bar{g}$. For this choice there are 10 oscillations at
$\omega<{\bar{g}}$ (panel a). Accordingly, $\eta_{\infty}(\omega)$ changes by
$2\pi$ five times (panel b) such that the total phase variation is $10\pi$.

We now move to $\omega>{\bar{g}}$ and use Eq. (32). At the largest $\omega$,
we have $B(\omega)\approx 1$ and

$\displaystyle
A(\omega)=\left(\frac{\bar{g}}{|\omega|}\right)^{\gamma}Q_{\gamma,\infty}\cos{\frac{\pi\gamma}{2}}$
$\displaystyle C(\omega)\approx
i\left(\frac{\bar{g}}{|\omega|}\right)^{\gamma}\sin{\frac{\pi\gamma}{2}}{\text{sign}}\omega{\bar{Q}}_{\gamma,\infty}$
(43)

where

$\displaystyle
Q_{\gamma,\infty}=\int_{0}^{\infty}d\omega_{m}\frac{\Delta_{\infty}(\omega_{m})}{\omega_{m}},$
$\displaystyle{\bar{Q}}_{\gamma,\infty}=\int_{0}^{\infty}d\omega\frac{\Delta_{\infty}(\omega)}{\omega}.$
(44)

The two integrals in (44) are actually identical. To prove this, we recall
that $D_{\infty}(\omega)=\Delta_{\infty}(\omega)/\omega$ satisfies KK relation
$D^{{}^{\prime\prime}}(\omega)=-(2/\pi)P\int_{0}^{\infty}dxD^{{}^{\prime}}(x)x/(x^{2}-\omega^{2})$,
where $P$ stands for principle part, and that
$D_{\infty}(\omega_{m})=(2/\pi)\int_{0}^{\infty}dxD^{{}^{\prime}}(x)x/(x^{2}+\omega^{2}_{m})$.
Using these relations, we obtain
$\int_{0}^{\infty}d\omega_{m}D_{\infty}(\omega_{m})=\int_{0}^{\infty}d\omega
D^{{}^{\prime}}(\omega)$. Using further that $\int_{0}^{\infty}d\omega
D^{{}^{\prime\prime}}(\omega)=(2/\pi)\int_{0}^{\infty}dxxD^{{}^{\prime}}(x)P\int_{0}^{\infty}d\omega/(x^{2}-\omega^{2})=0$
we find ${\bar{Q}}_{\gamma,\infty}=\int_{0}^{\infty}d\omega
D_{\infty}(\omega)=\int_{0}^{\infty}d\omega
D^{{}^{\prime}}(\omega)=\int_{0}^{\infty}d\omega_{m}D_{\infty}(\omega_{m})=Q_{\gamma}$.
Hence, at the largest $\omega$,

$A(\omega)+C(\omega)=e^{i(\pi\gamma/2){\text{s}gn}\omega}\left(\frac{\bar{g}}{|\omega|}\right)^{\gamma}$
(45)

Substituting into (32) and using (38), we obtain

$\Delta_{\infty}(\omega)\propto\left(\frac{\bar{g}}{|\omega|}\right)^{\gamma}e^{(i\pi\gamma/2){\text{sign}}\omega}Q_{\infty\gamma}$
(46)

This result could also be obtained by a direct rotation of
$\Delta_{\infty}(\omega_{m})\propto 1/|\omega_{m}|^{\gamma}$ from Matsubara to
real axis.

Figure 11:  (a) Comparison between the exact $\Delta_{\infty}(\omega)$ along
the real frequency axis (blue and orange thick lines) and the universal
contribution to $\Delta_{\infty}(\omega)$ from Eq. (51) (black dotted lines)
for $\gamma=1.91$. The agreement is nearly perfect. (b) Variation of the phase
of the gap function, $\eta_{\infty}(\omega)$, between $\omega\sim{\bar{g}}$
and $\omega=\infty$. For convenience of presentation we confined
$\eta_{\infty}(\omega)$ to $(-\pi,\pi)$, up to small variations. In the inset
we plot the continuous $\eta_{\infty}(\omega)$. We see that the total phase
variation between $\omega\sim\bar{g}$ and $\omega=\infty$ is
$16\pi+\pi\gamma/2$.

For $\gamma\geq 1$ there is a single crossover between the two limiting forms
of $\Delta_{\infty}(\omega)$, Eqs. (6) and (46). However, for $\gamma\leq 2$,
the new intermediate regime emerges at $\omega>{\bar{g}}$, as we now
demonstrate. In this regime, $\Delta^{\prime}_{\infty}(\omega)$ and
$\Delta^{{}^{\prime\prime}}_{\infty}(\omega)$ again oscillate, but with the
period set by a power of $\omega$ rather than by $\log|\omega|$.

To show this, we recall that in Sec. III.2 we found an exponentially small
oscillating component of $\Delta_{\infty}(\omega_{m})$ on the Matsubara axis
at $\omega_{m}>{\bar{g}}$. We now convert this component onto the real axis.
Replacing $\omega_{m}$ by $-i\omega+\delta$, we obtain for $\omega>0$:

$\Delta_{\infty}(\omega)=\int_{k_{min}}^{\infty}dk\left(e^{-\pi
k(1+\gamma/2)}e^{-ikS_{k}(\omega)}+e^{-\pi
k(1-\gamma/2)}e^{ikS_{k}(\omega)}\right)$ (47)

where $k_{min}=O(1)$ and
$S_{k}(\omega)=\log{\left[(|\omega|/{\bar{g}})^{\gamma}(e/A_{\gamma}k)^{\gamma-1}\right]}$
Comparing this result with Eq. (17) we see that the exponentially small factor
$e^{-\pi|k|}$ splits into $e^{-\pi k(1+\gamma/2)}$ and $e^{-\pi
k(1-\gamma/2)}$. The first term is smaller and the second one is larger than
the original term. Keeping only the larger term and evaluating the integral
over $k$ in the same way as in Sec. III.2, we obtain the oscillating term on
the real axis:

$\Delta_{\infty}(\omega)\propto|\omega|^{\frac{1}{2}\frac{\gamma}{\gamma-1}}\exp\left(-\frac{2-\gamma}{A_{\gamma}}\left(\frac{|\omega|}{\bar{g}}\right)^{\frac{\gamma}{\gamma-1}}\left[\frac{\pi}{2}-i\frac{\gamma-1}{2-\gamma}\left(1-\frac{\pi^{2}}{2(\gamma-1)^{2}}\right)\right]\right)$
(48)

For $\gamma\leq 2$, the exponential part of the prefactor is small in
$2-\gamma$, and the power-law part increases with $\omega$. As the result,
this oscillating contribution exceeds the one from Eq. (46) for frequencies
between $O({\bar{g}})$ and $\omega_{cr}$, where

$\omega_{cr}\sim{\bar{g}}\left(\frac{|\log{(2-\gamma)}|}{2-\gamma}\right)^{1/2}\gg{\bar{g}}.$
(49)

As a result, at ${\bar{g}}<\omega<\omega_{cr}$, $\Delta^{{}^{\prime}}(\omega)$
and $\Delta^{{}^{\prime\prime}}(\omega)$ oscillate as functions of
$(|\omega|/{\bar{g}})^{\gamma/(\gamma-1)}$ with weakly, yet exponentially
decaying prefactor. To the leading order in $2-\gamma$, Eq. (48) gives

$\Delta_{\infty}(\omega)\propto|\omega|e^{-(1-\gamma/2)(|\omega|/{\bar{g}})^{2}}e^{i(|\omega|/{\bar{g}})^{2}/\pi}$
(50)

We used $A_{2}=\pi$. For completeness, we computed the subleading term under
$e^{i...}$. It changes Eq. (50) into

$\Delta_{\infty}(\omega)\propto|\omega|e^{-(1-\gamma/2)(|\omega|/{\bar{g}})^{2}}e^{i\left((|\omega|/{\bar{g}})^{2}+\log{((|\omega|/{\bar{g}})^{2})}\right)/\pi}$
(51)

In Fig. 11 (a), we compare the exact $\Delta_{\infty}(\omega)$ with Eq. (51)
for $\gamma=1.91$. We see that the agreement is quite good, and the range of
$\omega^{2}$ oscillations is quite wide for this $\gamma$.

In Fig. 11 (b), we show the variation of the phase $\eta(\omega)$ between
$\omega\sim{\bar{g}}$ and $\omega=\infty$. We see that the phase changes by
$2\pi m$, where $m$ is an integer, and the total phase variation is $2\pi
m+\pi\gamma/2$. The last piece just follows from (46), and the first one is
due to oscillations given by (50). The integer $m$ increases one-by-one as
$\gamma$ increases towards $2$. In Sec. V below we associate $2\pi m$ phase
variation with the emergence of $m$ dynamical vortices in the upper half-plane
of frequency.

Figure 12: The density of states $N_{\infty}(\omega)$, counted from the normal
state value $N_{F}$. To better display oscillations we used different scales
at $\omega>\bar{g}$ and $\omega<\bar{g}$.

In Fig. 12, we plot the DoS $N_{\infty}(\omega)$ for several $\gamma$. The DoS
is defined as $N(\omega)=(-1/\pi){\text{Im}}G_{l}(\omega)$, where

$G_{l}(\omega)=-i\pi\sqrt{\frac{\omega^{2}}{\omega^{2}-\Delta^{2}(\omega)}}$
(52)

is the retarded Green’s function, integrated over the dispersion. 333 One can
verify that the extension of (52) to the upper frequency half-plane yields the
correct
$G(\omega_{m})=-i\pi{\text{s}gn}\omega_{m}\sqrt{\omega^{2}_{m}/(\omega^{2}_{m}+\Delta^{2}(\omega))}$.
We also caution that $\int d\omega N(\omega)$ is not equal to the number of
particles. The latter is expressed via a local time-ordered Green’s function.
Because $\Delta_{\infty}(\omega)$ is vanishingly small,

$N_{\infty}(\omega)=N_{F}\left(1-{\text{R}e}\frac{\Delta^{2}_{\infty}(\omega)}{\omega^{2}}\right)$
(53)

where $N_{F}$ is the DoS in the normal state. We see that $N_{\infty}(\omega)$
oscillates around $N_{F}$ up to $\omega\sim\omega_{cr}$. At
$\omega<{\bar{g}}$, the period of oscillations is set by $\log{\omega}$, and
at ${\bar{g}}<\omega<\omega_{max}$, the period is set by
$\omega^{\gamma/(\gamma-1)}$.

###  IV.3 Frequency dependence of $\Delta_{n}(\omega)$. The case $n=0$.

We now consider the opposite limit $n=0$. We remind that on the Matsubara axis
$\Delta_{0}(\omega_{m})$ is sign-preserving. We show that on the real axis
both real and imaginary parts of $\Delta_{0}$ again oscillate in a finite
frequency range at $\omega>{\bar{g}}$, and the phase of complex
$\Delta_{0}(\omega)$ is winding by $2\pi m$ in this range.

At high frequencies $\omega\gg{\bar{g}}$, the form of $\Delta_{0}(\omega)$ can
be readily obtained from (37) in the same way as it was done in the previous
section. For the case $n=0$ we have

$\displaystyle
A(\omega)=\left(\frac{\bar{g}}{|\omega|}\right)^{\gamma}Q_{\gamma}\cos{\frac{\pi\gamma}{2}}$
$\displaystyle C(\omega)\approx
i\left(\frac{\bar{g}}{|\omega|}\right)^{\gamma}Q_{\gamma}\sin{\frac{\pi\gamma}{2}}{\text{sign}}\omega,$
(54)

where $Q_{\gamma}$ is given by (29). Substituting these forms along with
$B(\omega)\approx 1$ into (32), we obtain at $\omega\gg{\bar{g}}$,

$\Delta_{0}(\omega)\approx\left(\frac{\bar{g}}{|\omega|}\right)^{\gamma}e^{(i\pi\gamma/2){\text{sign}}\omega}Q_{\gamma}$
(55)

This is consistent with
$\Delta_{0}(\omega_{m})\propto({\bar{g}}/|\omega_{m}|)^{\gamma}$ on the
Matsubara axis.

At small $\omega$,
$\Delta_{0}(\omega)=(A(\omega)+C(\omega))/B(\omega)\approx\Delta_{0}(0)$ where
$\Delta_{0}(0)$ is the same as we found in the $\omega_{m}\to 0$ limit on the
Matsubara axis. We note, however, that the relative strength of
$A(\omega),B(\omega)$, and $C(\omega)$ changes between $\gamma<1$ and
$\gamma>1$. For $\gamma<1$, $A(\omega)$ and $B(\omega)$ tend to finite values
at $\omega\to 0$, while $C(\omega)\propto\omega^{1-\gamma}$ vanishes. Then
$\Delta_{0}(0)=A(0)/B(0)$. For $\gamma>1$, $B(\omega)$ and $C(\omega)$ scale
as $\omega^{1-\gamma}\sin(\pi\gamma/2)/(2-\gamma)$ at $\omega\to 0$, while
$A(\omega)\propto\omega^{1-\gamma}\sin(\pi\gamma/2)/(1-\gamma)$ is smaller, at
least for $\gamma\leq 2$, and, moreover, its sign changes to negative at
$\gamma>1$. The gap still tends to a finite positive value at $\omega=0$
because once $A(\omega)\ll C(\omega)$, $\Delta_{0}(0)=\lim_{\omega\to
0}C(\omega)/B(\omega)={\text{c}onst}$. However, we see that $\Delta_{0}(0)$
now scales with $C(\omega)$, which, we recall, is present in the gap equation
on the real axis because of the need to cancel out parasitic contributions
from the branch cut. This hints that for $\gamma>1$ the behavior of
$\Delta_{0}(\omega)$ on the real axis may be quite different from that of
$\Delta_{0}(\omega_{m})$.

We now show that this is indeed the case. We focus on frequencies
$\omega>{\bar{g}}$, where we detected the new behavior at $n=\infty$. Eq. (54)
for $A(\omega)$ is valid for all $\omega>{\bar{g}}$ as typical $\omega_{m}$ in
(35) are of order ${\bar{g}}$. Similarly, $B(\omega)\approx 1$ for all
$\omega>{\bar{g}}$. For $C(\omega)$, we have to be more careful and include
not only the contribution in (37) from $\omega-\Omega=O({\bar{g}})$, which
yields (38), but also contributions from $\Omega$ immediately below $\omega$
(i.e., $\omega-\Omega\ll\bar{g}$). These contributions are expressed in terms
of derivatives of $D_{0}(\omega)=\Delta_{0}(\omega)/\omega$. The leading term
here is the one with the first derivative:
$D_{0}(\omega-\Omega)-D_{0}(\omega)\approx-\Omega{\dot{D}}_{0}(\omega)$. For
$\gamma\leq 2$, the integral
$\int_{0}^{\omega}d\Omega/\Omega^{\gamma-1}=\omega^{2-\gamma}/(2-\gamma)$ is
determined by small $\Omega$, and the prefactor $1/(2-\gamma)$ compensates the
smallness of the overall factor $\sin(\pi\gamma/2)\approx(\pi/2)(2-\gamma)$ in
(37). This term then gives

$-i\frac{\pi}{2}\omega^{2-\gamma}\frac{{\dot{D}}_{0}(\omega)}{\sqrt{1-D^{2}_{0}(\omega)}}$
(56)

We will also need terms of order $2-\gamma$. For this, we keep the subleading
terms with $({\dot{D}}_{0})^{2}$ and ${\ddot{D}}_{0}$. Evaluating the
prefactors, we obtain differential gap equation in the form ($\omega>0$)

$-i\frac{\pi}{2}{\bar{g}}^{\gamma}\frac{\omega^{2-\gamma}}{\sqrt{1-D^{2}_{0}(\omega)}}\left[{\dot{D}}_{0}(\omega)-(2-\gamma)\omega\left(\frac{1}{2}{\ddot{D}}_{0}(\omega)+\frac{D_{0}(\omega)({\dot{D}}_{0}(\omega))^{2}}{1-D^{2}_{0}(\omega)}\right)\right]=D_{0}(\omega)\omega-
Q_{\gamma}\left(\frac{\bar{g}}{\omega}\right)^{\gamma}e^{i\pi\gamma/2}$ (57)

We follow Refs. Karakozov  _et al._ (1991); Combescot (1995) and introduce
$D_{0}(\omega)=1/\sin(\phi(\omega))$. Both $D_{0}(\omega)$ and $\phi(\omega)$
are complex functions of $\omega$. Substituting into (57) we obtain

${\dot{\phi}}-\frac{(2-\gamma)}{2}\omega\left(\ddot{\phi}+({\dot{\phi}})^{2}\tan{\phi}\right)=\frac{2}{\pi{\bar{g}}^{\gamma}}\left(\omega^{\gamma-1}-Q_{\gamma}\frac{\bar{g}^{\gamma}}{\omega^{2}}e^{i\pi\gamma/2}\sin{\phi}\right)$
(58)

This equation is similar to the one for $\gamma=2$ and a finite $T$, analyzed
by Combescot in Ref. Combescot (1995).

At the highest frequencies, the gap function must obey Eq. (46). This gap
function is reproduced if we choose

$\displaystyle\phi^{\prime}(\omega)=2\pi m+\frac{\pi}{2}(\gamma+1)$ (59)
$\displaystyle\phi^{{}^{\prime\prime}}(\omega)=\log{\frac{2{\bar{g}}}{Q_{\gamma}}}+(\gamma+1)\log{\frac{\omega}{{\bar{g}}}}$
(60)

where $m$ is integer. We see that at large enough $\omega$,
$\phi^{{}^{\prime\prime}}(\omega)\gg\phi^{{}^{\prime}}(\omega)$. On the other
hand, at $\omega\geq{\bar{g}}$,
$\phi^{\prime}(\omega)\approx(2/\pi\gamma)(\omega/{\bar{g}})^{\gamma}$, and
$\phi^{{}^{\prime\prime}}$ is small, of order $(2-\gamma)$. We use this to set
the boundary condition at $\omega={\bar{g}}$ as

$\phi^{\prime}=\frac{2}{\pi\gamma},~{}~{}\phi^{{}^{\prime\prime}}=a(2-\gamma),$
(61)

where $a=O(1)$. We will argue that the solution of (58) is largely independent
on $a$, as long as $a(2-\gamma)\ll 1$. To simplify the calculation, below we
neglect $\ddot{\phi}$ term in Eq. (58) and use the boundary condition (61) as
the initial condition for the first order differential equation. We show that
the solution of (58) without $\ddot{\phi}$ by itself satisfies the boundary
condition (59,60). We discuss the role of $\ddot{\phi}$ and the validity of
dropping it at the end of this section.

A simple analysis of the Eq. (58) without $\ddot{\phi}$ term shows that
$\phi^{{}^{\prime\prime}}(\omega)$ rapidly increases shortly before
$\phi^{\prime}(\omega)$ reaches $\pi/2$. To see this, we neglect momentarily
the $Q_{\gamma}$ term in (58), which gives rise to a small initial
$\phi^{{}^{\prime\prime}}$ and solve the remaining equation as a quadratic
equation on $\dot{\phi}$. We obtain

$\dot{\phi}=\frac{1-\sqrt{1-\frac{4}{\pi}(2-\gamma)\left(\frac{\omega}{\bar{g}}\right)^{\gamma}\tan\phi(\omega)}}{(2-\gamma)\omega\tan\phi(\omega)}.$
(62)

(the sign is chosen to satisfy the initial condition). An elementary analysis
shows that $\phi^{{}^{\prime\prime}}(\omega)$ emerges at $\omega=\omega_{a}$,
where when

$\left(\frac{\omega_{a}}{\bar{g}}\right)^{\gamma}\tan(\phi^{\prime}(\omega_{a}))=\frac{\pi}{4}\frac{1}{2-\gamma},$
(63)

This $\omega_{a}$ is smaller than the one at which $\phi^{\prime}(\omega)$
reaches $\pi/2$. This is essential as in the absence of
$\phi^{{}^{\prime\prime}}$, the behavior near $\phi^{\prime}(\omega)=\pi/2$
would be singular. Once $\phi^{{}^{\prime\prime}}$ is non-zero, the
singularity is cut. Keeping the $Q_{\gamma}$ term, we find that
$\phi^{{}^{\prime\prime}}$ is non-zero at all frequencies, where Eq.(58) is
valid, but still rapidly increases around $\omega_{a}$, specified by (63).

At larger $\omega$, $\phi^{{}^{\prime\prime}}(\omega)$ increases, and
eventually $e^{\phi^{{}^{\prime\prime}}}$ becomes larger than one. At such
frequencies, $\tan{\phi}\approx i$, and (58) simplifies to

$\displaystyle{\dot{\phi}}^{\prime}(\omega)\approx\frac{2}{\pi}\frac{\omega^{\gamma-1}}{{\bar{g}}^{\gamma}}$
$\displaystyle{\dot{\phi}}^{{}^{\prime\prime}}(\omega)\approx\frac{2(2-\gamma)}{\pi^{2}}\frac{\omega^{2\gamma-1}}{{\bar{g}}^{2\gamma}}$
(64)

Solving, we find
$\phi^{\prime}(\omega)\approx(2/\pi\gamma)(\omega/{\bar{g}})^{\gamma}$,
$\phi^{{}^{\prime\prime}}(\omega)\approx(2-\gamma)/(\pi^{2}\gamma)(\omega/{\bar{g}})^{2\gamma}$.
Note that these forms are universal and do not depend on the boundary
condition (i.e., on the factor $a$ in Eq. (61)). For $\gamma\approx 2$,
$\phi^{\prime}(\omega)\approx(\omega/{\bar{g}})^{2}/\pi$. To compare with
$n=\infty$ case, we computed the subleading term. It comes from the second
term in $B(\omega)$ in (38) and changes $\phi^{\prime}(\omega)$ to

$\phi^{\prime}(\omega)=\frac{1}{\pi}\left(\left(\frac{\omega}{\bar{g}}\right)^{2}+\log{\left(\frac{\omega}{\bar{g}}\right)^{2}}\right).$
(65)

Observe that the r.h.s. of (65) is the same function as we found for
$n=\infty$, Eq. (51).

Figure 13: Panel (a). The numerical solution of Eq. (58) for $\gamma=1.995$.
Frequency is in units of ${\bar{g}}$. The crossover at $\omega_{cr}\approx
15{\bar{g}}$ is clearly visible. Panels (b) and (c) – zoom into regions
$\omega<\omega_{cr}$ and $\omega>\omega_{cr}$. At $\omega<\omega_{cr}$,
$\phi^{\prime}$ grows as $\omega^{\gamma}$, while $\phi^{\prime\prime}$ first
increases step-like, and then behaves as $\omega^{2\gamma}$. At
$\omega>\omega_{cr}$, $\phi^{\prime}$ saturates at $2\pi m+(\gamma+1)\pi/2$,
where $m=9$ for $\gamma=1.995$, while $\phi^{{}^{\prime\prime}}$ increases as
$(\gamma+1)\log{\omega/\bar{g}}$.

Eq. (64) holds up to $\omega$ at which
$e^{\phi^{{}^{\prime\prime}}}\sim(\omega/{\bar{g}})^{3}$. At larger
frequencies, the second term in the r.h.s. of (58) cannot be neglected, and
the functional form of $\phi(\omega)$ changes. In Fig.13 (a) we show the
result of numerical solution of (58) We see that there is a single crossover,
at
$\omega=\omega_{cr}\sim{\bar{g}}(|\log{(2-\gamma)}|/(2-\gamma))^{1/(2\gamma)}$,
from $\phi(\omega)$ given by (64) to $\phi(\omega)$ given by (59,60). In
Fig.13 (b) and (c) we show separately the behavior of $\phi^{\prime}(\omega)$
and $\phi^{{}^{\prime\prime}}(\omega)$ at $\omega<\omega_{cr}$ and
$\omega>\omega_{cr}$. The fits to Eq. (64) Eqs. (59,60), respectively, are
almost perfect. Deviations from the high-frequency forms decrease as
$1/\omega^{\gamma}$ and oscillate as trigonometric functions of
$(2/\pi\gamma)(\omega/{\bar{g}})^{\gamma}$.

Figure 14:  (a)-(b) The functions $\phi^{\prime}(\omega)$ and
$\phi^{\prime\prime}(\omega)$ for different boundary conditions, set by a
parameter $a$ in eq.(61). The figure shows that the value $2\pi
m+(\gamma+1)\pi/2$, at which $\phi^{\prime}(\omega)$ saturates is independent
on $a$. The behavior of $\phi^{\prime}(\omega)$ in the universal regime, where
both $\phi^{\prime}$ and $\phi{{}^{\prime\prime}}$ are continuous functions of
$\omega$, also does not depend on $a$. (c) Variation of the integer $m$ with
$\gamma$. There is a discrete set of $\gamma_{i}$, at which $m$ changes by
$1$. The set becomes progressively more dense at $\gamma\to 2$. (d) The
behavior of $\phi^{\prime}(\omega)$ near one of these $\gamma_{i}\approx
1.924$. The value of $\phi^{\prime}$ at large $\omega$ jumps by $2\pi$ as
$\gamma$ passes through $\gamma_{i}$ and a new vortex moves into the upper
half-plane of frequency.

In Fig.14((a) and (b)) we show that the value of $\phi^{\prime}$ for
$\omega>\omega_{cr}$ is independent on the parameter $a$ in the boundary
condition (61), as long as $a=O(1)$. At the same time, the value of an integer
$m$ in (59) changes if we change $\gamma$, as is shown in Fig.14(c).
Specifically, $m$ jumps to the nearest integer at a discrete set of
$\gamma_{i}\leq 2$ (the smaller is $2-\gamma$, the larger is $m$). We
demonstrate this in Fig.14(d).

In Fig 15(a) we show real and imaginary parts of the gap function

$\Delta_{0}(\omega)=\frac{\omega}{\sin{\phi(\omega)}}=2\omega\frac{\sin{\phi^{\prime}}\cosh{\phi^{{}^{\prime\prime}}}-i\cos{\phi^{\prime}}\sinh{\phi^{{}^{\prime\prime}}}}{\cosh{2\phi^{{}^{\prime\prime}}}-\cos{2\phi^{\prime}}}$
(66)

We see that both $\Delta^{\prime}_{0}(\omega)$ and
$\Delta^{{}^{\prime\prime}}_{0}(\omega)$ oscillate between
$\omega=O({\bar{g}})$ and $\omega_{cr}$ and display sign-preserving
$1/\omega^{\gamma}$ behavior at $\omega>\omega_{cr}$. At $\omega\geq\bar{g}$,
$\Delta^{\prime}(\omega)$ varies roughly as
$\omega/\sin{((2/(\pi\gamma))(\omega/{\bar{g}})^{\gamma})}$ and
$\Delta^{{}^{\prime\prime}}(\omega)$ has almost $\delta$-functional spikes in
where $\sin(2/(\pi\gamma))(\omega/{\bar{g}})^{\gamma})$ is small. At larger
$\omega\leq\omega_{cr}$, both $\Delta^{\prime}$ and
$\Delta^{{}^{\prime\prime}}$ oscillate with progressively decreasing
magnitude.

Figure 15:  Upper panel: frequency dependence of the gap function
$\Delta_{0}(\omega)=\Delta_{0}^{\prime}(\omega)+i\Delta_{0}^{\prime\prime}(\omega)$,
from Eq. (66). We set $\gamma=1.95$. For $\bar{g}<\omega<\omega_{cr}$, both
$\Delta_{0}^{\prime}(\omega)$ and $\Delta_{0}^{\prime\prime}(\omega)$ display
oscillations with a decreasing magnitude. For $\omega\geq\omega_{cr}$,
$\Delta_{0}^{\prime}(\omega)$ becomes negative and does not oscillate. Lower
panel: Variation of the phase of the gap $\eta_{0}(\omega)$ as a function of
$\omega$. As before, we confined phase variation to $(-\pi,\pi)$, up to small
variations. The phase changes by $2\pi$ three times between
$\omega=O(\bar{g})$ and $\omega_{cr}$. There are no more $2\pi$ phase
variations, despite that the phase shows wiggling at large frequencies. In the
inset we plot the continuous $\eta_{0}(\omega)$. We see that the total phase
variation between $\omega_{m}=\bar{g}$ and $\omega_{m}=\infty$ is
$6\pi+\pi\gamma/2$.

In Fig. 15(b) we show the variation of the phase of the gap function
$\eta(\omega)$. The total variation of $\eta$ between $\omega=0$ and
$\omega=\infty$ is $2\pi m+\pi\gamma/2$. We emphasize that this is the result
for $\Delta_{0}(\omega)$, which is sign-preserving on the Matsubara axis. We
clearly see that there is a huge difference between the forms of the gap
function on the real and the Matsubara axis in between ${\bar{g}}$ and
$\omega_{cr}$.

We now use more precise analysis to determine $\omega_{cr}$. In Eq. (58) we
only included the first two terms in the expansion of $C(\omega)$ in powers of
$\dot{\phi}$. Meanwhile, the expansion in derivatives holds in powers of
${\dot{\phi}}\omega/{\bar{g}}$ without additional $(2-\gamma)$ in the
prefactor. This implies that higher-order terms are not negligible at
$\omega>{\bar{g}}$. We now use the fact that before $\phi(\omega)$ crosses
over to (54,38), it shows the universal behavior in the range where
$1<e^{\phi^{{}^{\prime\prime}}}<(\omega/{\bar{g}})^{3}$. In this regime, (i)
${\dot{\phi}}\approx{\dot{\phi}}^{\prime}\approx(2\pi)\omega^{\gamma-1}/{\bar{g}}^{\gamma}\gg{\bar{g}}$,
while higher-order derivatives are smaller, and (ii) $D(\omega)$ is small,
such that $\sqrt{1-D^{2}}\approx 1$. In this situation, one can sum up the
full Taylor series for $(2-\gamma)$ term in $C(\omega)$. We use that
$\tan\phi\approx i$ and the $n$-th detivative $D^{n}(\omega)\approx
i^{n}(\dot{\phi})^{n}D$. Integrating each term in Talor series over $\Omega$,
we obtain the modified l.h.s. of (58) in the form

${\dot{\phi}}-i(2-\gamma){\dot{\phi}}F\left(\frac{\omega}{\bar{g}}\dot{\phi}\right)$
(67)

where

$F(x)=\int_{0}^{x}\frac{1-\cos{y}}{y^{2}}=SI(x)-\frac{1-\cos{x}}{x}$ (68)

and $SI(x)$ is SinIntegral. Eq. (58) is reproduced if we approximate
$1-\cos{y}$ by $y^{2}/2$. Then $F(x)\approx x/2$. If we use the full
expression, we find that at large $x$, which we are interested in,
$F(x)\approx\pi/2$. The equation for $\phi$ then reduces to

$\dot{\phi}=\frac{2}{\pi\bar{g}^{\gamma}}\frac{\omega^{\gamma-1}-\frac{Q_{\gamma}\bar{g}^{\gamma}}{\omega^{2}}e^{i\pi\gamma/2}\sin\phi}{1-i(2-\gamma)\pi/2}$
(69)

We show the solution in Fig.16 and present the plots of the gap function
$\Delta_{0}(\omega)$, the phase $\eta_{0}(\omega)$, and the variation of $m$
with $\gamma$ in Fig.17.

Figure 16:  Panel (a). The solution of Eq. (69) for $\gamma=1.88$. Panels (b)
and (c) - zooms into the regions below and above the crossover frequency. The
behavior of $\phi^{\prime}(\omega)$ and $\phi^{{}^{\prime\prime}}(\omega)$ is
qualitatively similar to that in Fig. 14, but differs in detail. In
particular, both $\phi^{\prime}(\omega)$ and
$\phi^{{}^{\prime\prime}}(\omega)$ approach asymptotic forms
$12\pi+(\gamma+1)\pi/2$ and $(\gamma+1)\log{\omega}$, respectively, without
oscillations.

The behavior of $\phi(\omega)$ is qualitatively similar to the one in Fig.13,
i.e., there is a single crossover frequency, and the total variation of
$\eta(\omega)$ between $\omega=0$ and $\omega=\infty$ is $2\pi m+\pi\gamma/2$.
However the crossover scale $\omega_{cr}$ has different dependence on
$2-\gamma$ compared to Eq. (58), and also the dependence of $m$ on $\gamma$ is
different from that in Fig.14(b) (and there are no wiggles in $\phi(\omega)$
at $\omega>\omega_{cr}$).

Figure 17: Panels (a) and (b) – the same as in Fig. 15, but for
$\phi^{\prime}$ and $\phi^{{}^{\prime\prime}}$ from Eq.(69). At large
frequencies, $\Delta(\omega)\propto e^{i\pi\gamma/2}/\omega^{\gamma}$. The
total variation of the phase of $\Delta$ is $\delta\eta=2\pi m+\pi\gamma/2$.
Panel (c) – variation of $m$ with $\gamma$. As before, there is a discrete set
of $\gamma_{i}$, where $m$ jumps by 1. The set becomes progressively more
dense at $\gamma\to 2$.

To obtain the modified crossover frequency, $\omega_{cr}$, we note that in the
universal regime $1<e^{\phi^{\prime\prime}}<(\omega/{\bar{g}})^{3}$,
${\dot{\phi}}^{\prime}(\omega)$ is still given by (64), while

${\dot{\phi}}^{{}^{\prime\prime}}(\omega)\approx(2-\gamma)\frac{\omega^{\gamma-1}}{{\bar{g}}^{\gamma}}$
(70)

Accordingly,
$\phi^{\prime}(\omega)\approx(2\pi\gamma)(\omega/{\bar{g}})^{\gamma}$,
$\phi^{{}^{\prime\prime}}(\omega)=(2-\gamma)/\gamma(\omega/{\bar{g}})^{\gamma}$.
The crossover frequency is then determined by
$e^{\phi^{{}^{\prime\prime}}}\sim(\omega/{\bar{g}})^{3}$. Solving for
$\omega$, we obtain, for $\gamma\leq 2$,

$\omega_{cr}\sim\bar{g}\left(\frac{|\log{(2-\gamma)}|}{2-\gamma}\right)^{1/2}$
(71)

This is the same scale as $\omega_{cr}$ that we obtained for $n=\infty$, Eq.
(49).

We see therefore that for $n=\infty$ and $n=0$ the new behavior of the gap
function emerges in the same frequency range ${\bar{g}}<\omega<\omega_{cr}$.
Moreover, the period of oscillations in this range is the same function of
frequency, Eqs. (51) and (65).

In Fig.18 we plot the DoS

$N_{0}(\omega)=N_{F}{\text{R}e}\sqrt{\frac{\omega^{2}}{\omega^{2}-\Delta^{2}_{\infty}(\omega)}}$
(72)

for several $\gamma$.

Figure 18: The DoS $N(\omega)$ for the $n=0$ solution, for several $\gamma$.
Peaks in $N(\omega)$ sharpen as $\gamma\to 2$, while the weight between the
peaks is reduced.

We see that $N_{0}(\omega)$ vanishes below a certain threshold frequency, and
at ${\bar{g}}<\omega<\omega_{max}$ has a set of maxima and minima. The peaks
in $N_{0}(\omega)$ get progressively sharper as $\gamma$ approaches $2$, while
in between the peaks $N_{0}(\omega)$ gets smaller. The positions of the peaks
almost coincide with the maxima of $N_{\infty}(\omega)$.

Before we move to finite $n$, we pause to discuss the validity of neglecting
the $\ddot{\phi}$ term in Eq. (58). This is justified if $|\ddot{\phi}|$ is
smaller than $(\dot{\phi})^{2}|\tan\phi|$ for all frequencies. Of particular
relevance here are frequencies near $\omega_{a}$, specified by Eq. (63). If we
set $a=0$ in the initial condition (61) and use Eq. (62), we find that
$\phi^{{}^{\prime\prime}}$ initially increases as $(\omega-\omega_{0})^{3/2}$.
In this situation, $\ddot{\phi}$ diverges at $\omega=\omega_{0}+0$, and the
behavior of $\phi(\omega)$, which we found earlier in this Section, is valid
outside of the vicinity of $\omega_{0}$, while near $\omega_{0}$,
$\ddot{\phi}$ and higher derivatives cannot be neglected. Specifically,
evaluating $\ddot{\phi}$, $\dot{\phi}$, and $\tan\phi$ near $\omega_{0}$, we
find that for $\gamma\leq 2$, the condition
$|\ddot{\phi}|<(\dot{\phi})^{2}|\tan\phi|$ is satisfied when

$\phi^{\prime\prime}>\pi(2-\gamma)$ (73)

For $a=0+$, this holds outside a finite range near $\omega_{0}$. However, if
$a=O(1)$, $|\ddot{\phi}|$ never becomes much larger than
$(\dot{\phi})^{2}|\tan\phi|$, hence neglecting $|\ddot{\phi}|$ does not change
the results qualitatively even near $\omega_{0}$.

###  IV.4 Frequency dependence of $\Delta_{n}(\omega)$. Finite $n$

The behavior of $\Delta_{n}(\omega)$ at $\omega<{\bar{g}}$ on the real axis
parallels the one on the Matsubara axis. Namely, the phase varies by $n\pi$
between $\omega=0$ and ${\bar{g}}$ due to n vortices on the Matsubara axis. At
$\gamma\to 2$, oscillations shift to smaller $\omega$. At $\omega>{\bar{g}}$,
it is natural to expect oscillations with the period set by (51) at
${\bar{g}}<\omega<\omega_{cr}$ and $e^{i\pi\gamma/2}/\omega^{\gamma}$ behavior
at $\omega>\omega_{cr}$. The total phase variation between $\omega=0$ and
$\omega=\infty$ is

$n\pi+2\pi m+\frac{\pi\gamma}{2}$ (74)

The DoS at a finite $n$ vanishes below a certain threshold frequency. Above
the threshold, it oscillates around $N_{F}$ with a period set by
$\log{\omega}$, and at larger ${\bar{g}}<\omega<\omega_{cr}$ has a set of
peaks and dips at about the same frequencies as $N_{\infty}(\omega)$ and
$N_{0}(\omega)$.

##  V Vortices near the real axis

We now use the Cauchy relation

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

and extend the gap functions $\Delta_{n}(\omega)$ into the upper frequency
half-plane, to $z=\omega^{\prime}+i\omega^{{}^{\prime\prime}}$
($\omega^{{}^{\prime\prime}}>0$). Earlier, we demonstrated that
$\Delta_{n}(z)$ has $n$ vortices on the Matsubara axis. Here, we analyze the
behavior of $\Delta_{n}(z)$ between the Matsubara and the real axis. We show
that as $\gamma$ increases from one to two, new vortices appear one-by-one in
the upper frequency half-plane, near the real axis. These vortices are located
at $|z|>{\bar{g}}$, and their number, $m$, is determined by $\gamma$ and is
the same for all $n$. When $\gamma\to 2$, $m$ tends to infinity. The emergence
of vortices obviously correlates with the oscillations of $\Delta_{n}(\omega)$
on the real axis, and as such is another consequence of the change of sign of
the real part of the interaction on the real axis, $V^{\prime}(\Omega)$, which
becomes repulsive at $\gamma=1$. The increase of the number of vortices as
$\gamma$ approaches $2$ in turn correlates with the decrease of
$V^{{}^{\prime\prime}}(\Omega)$.

That such vortices must be present can be understood by comparing the behavior
of $\Delta_{n}(z)$ for $|z|>{\bar{g}}$ near the Matsubara axis and the real
axis. Along the Matsubara axis, $\Delta_{n}(\omega_{m})$ is real and does not
change sign at $|\omega_{m}|>{\bar{g}}$. By continuity, $\Delta^{\prime}(z)$
near the Matsubra axis should remain sign-preserving, hence $\eta(z)$ does not
wind. On the other hand, on the real axis, the phase winds by $2\pi m$, as we
found in the previous Section. This number is topologically protected against
perturbations and can only disappear upon rotation from real to imaginary $z$
if there are vortices at some complex $z$. Indeed, let’s compute

$\delta\eta_{\Gamma}=\text{Im}\int_{\Gamma}dz\left[\partial\log\Delta(z)/\partial
z\right]$ (76)

along the path $\Gamma$, which starts at the large negative real $z=-R$, goes
along the real axis up to $+R$, and then closes along the large arc in the
upper half plane (see Fig.19)

Figure 19: The contour $\Gamma$ used in Eq.(76). The contour consists of a
real axis and an semi-circle(arc) at $z=|z|e^{i\psi}$, $|z|\to\infty$ and
$0<\psi<\pi$. Along the arc, $\Delta_{n}(z)\propto
e^{i(\pi/2-\psi)\gamma}/|z|^{\gamma}$, so the corresponding contribution to
$\delta\eta_{L}$ is $-\pi\gamma$. Along the real axis,
$\delta\eta_{\Gamma}=2\pi(n+2m)+\pi\gamma$. The total $\delta\eta_{\Gamma}$
along $\Gamma$ is then $2\pi(n+2m)$. The same $\delta\eta_{\Gamma}$ must be
obtained by counting contributions from the poles inside $\Gamma$. Each vortex
is a pole with residue $2\pi$, hence there should be $n+2m$ vortices. This is
consistent with our analysis of $\Delta(z)$: there are $n$ vortices on the
Matsubara axis, and $2m$ vortices in the upper half plane near the real axis.

The arc is chosen such that along it $\Delta_{n}(z)\propto
e^{i(\pi/2-\psi)}/|z|^{\gamma}$ . The total phase variation along the arc is
$-\pi\gamma$, and the phase variation along the real axis is
$2\pi(n+2m)+\pi\gamma$. The contour integral then gives

$\delta\eta_{\Gamma}=2\pi(n+2m)$ (77)

This $\delta\eta_{\Gamma}$ should be equal to the contribution from inside the
contour. Because $\partial\log\Delta(z)/\partial z$ has a simple pole at each
point where $\Delta(z)$ vanishes, there must be $n+2m$ vortices inside the
contour (we recall that $\Delta(z)$ by itself has no poles in the upper half
plane). There are $n$ vortices on the Matsubara axis, the other $2m$ should be
located in between the Matsubara and the real axis. By symmetry, there must
then be $m$ vortices in the first quadrant, and another $m$ is in the second
one.

To determine the location of the vortices, it is instructive to again consider
separately the cases $n=\infty$, $n=0$, and a finite $n$.

###  V.1 $n=\infty$

We first quickly verify that there are no vortices away from Matsubara axis
for $|z|<{\bar{g}}$. We express $\Delta(z)$ in terms of $|z|$ and $\psi$,
defined via $z=|z|e^{\psi}$. Replacing $\omega_{m}$ by $-iz$ in Eq. (7) we
obtain the expressions similar to (41) and (42):

$\displaystyle\Delta^{\prime}_{\infty}(z)\propto\cos{\left(\frac{\gamma}{2}\left(\frac{\pi}{2}-\psi\right)\right)}\cos{D_{1}}\cosh{D_{2}}+\sin{\left(\frac{\gamma}{2}\left(\frac{\pi}{2}-\psi\right)\right)}\sin{D_{1}}\sinh{D_{2}}$
$\displaystyle\Delta^{{}^{\prime\prime}}_{\infty}(z)\propto\cos{\left(\frac{\gamma}{2}\left(\frac{\pi}{2}-\psi\right)\right)}\sin{D_{1}}\sinh{D_{2}}-\sin{\left(\frac{\gamma}{2}\left(\frac{\pi}{2}-\psi\right)\right)}\cos{D_{1}}\cosh{D_{2}}$
(78)

but $D_{1}$ and $D_{2}$ are now given by

$\displaystyle D_{1}=\beta\log{(|z|/{\bar{g}})^{\gamma}}+\phi,$ (79)
$\displaystyle D_{2}=\beta\left((\pi-2\psi)\frac{\gamma}{2}\right)$ (80)

The vortices are the points where
$\Delta^{\prime}_{\infty}(z)=\Delta^{{}^{\prime\prime}}_{\infty}(z)=0$. For
$\Delta_{\infty}(z)$ given by (78) this is satisfied if $D_{1}=\pi/2+m\pi$ and
$D_{2}=0$. The second condition is satisfied only if $\psi=\pi/2$, i.e., on
the Matsubara axis. The first condition coincides with
$\Delta_{\infty}(\omega_{m})=0$. In Fig.20 we show that the phase
$\eta_{\infty}(z)$ of
$\Delta_{\infty}(z)=|\Delta_{\infty}(z)|e^{i\eta_{\infty}(z)}$ evolves in a
”rectangular” way upon rotation from the Matsubara to the real axis. The white
curves in this Figure are determined by $\Delta^{{}^{\prime\prime}}(z)=0$ and
$\Delta^{\prime}(z)<0$. The phase $\eta(z)$ changes by $2\pi$ upon crossing
each of these curves. We see that the phase winds five times by $2\pi$ at
Re$z>0$, consistent with the presence of 10 vortices on the Matsubara axis.

Figure 20: The phase $\eta_{\infty}(z)$ in the complex plane at
$|z|<O(\bar{g})$, for $\gamma=1.96$. We set the lower boundary for $|z|$ at
$\omega_{0}\sim 10^{-18}\bar{g}$, as in Fig.10. Along the white lines,
$\Delta^{{}^{\prime\prime}}(z)=0$ and $\Delta^{\prime}(z)<0$. The phase of the
gap then changes by $2\pi$ upon crossing each of these lines. Observe that
each white line starts in between a pair of vortices on the Matsubara axis,
and ends on the real axis. As a result, if there are $k$ vortices on the
Matsubara axis, the phase variation $\delta\eta$ on the positive real axis is
exactly $k\pi$. In our case, $k=10$.  Figure 21: The case $n=\infty$. We plot
the phase of the gap function $\eta_{\infty}(z)$ in the upper half plane
($z=\omega^{\prime}+i\omega^{\prime\prime}$) for different $\gamma$ near
$\gamma=2$. The locations of the vortices are marked by red dots. At the
vortex core, $\Delta_{\infty}(z)=0$ and $\eta_{\infty}$ is undefined. This
number of vortices $m$ rapidly increases as $\gamma\to 2$.

We now consider the range $|z|>{\bar{g}}$. We recall that on the real axis,
$\Delta_{\infty}(\omega)$ oscillates as a function of
$(|\omega|/{\bar{g}})^{\gamma/(\gamma-1)}$ between ${\bar{g}}$ and a larger
$\omega_{cr}$ given by Eq. (49). We detected the oscillating term on the
Matsubara axis $(\psi=\pi/2$) , where it is exponentially small, and converted
it to the real axis ($\psi=0$), where it becomes the largest piece in
$\Delta_{\infty}(\omega)$ at ${\bar{g}}<\omega<\omega_{cr}$. We now analyze
this term at arbitrary $\psi$. Replacing $\omega_{m}$ by
$|z|e^{-i(\pi/2-\psi)}$ in (18) and using (6), we obtain the oscillating
component of $\Delta_{\infty}(z)$ in the form

$\Delta_{\infty}(z)=\int_{k_{min}}^{\infty}dk\bigg{(}e^{-(\pi+\theta)k}e^{-ikS_{k}(|z|)}+e^{-(\pi-\theta)k}e^{ikS_{k}(|z|)}\bigg{)}.$
(81)

Here, $\theta=(\pi/2-\psi)\gamma$ and
$S_{k}(|z|)=\log{\left[(|z|/{\bar{g}})^{\gamma}(e/A_{\gamma}k)^{\gamma-1}\right]}$,
where $A_{\gamma}$ is defined in (15). For $\gamma\to 2$,
$\theta\to\pi-2\psi$, and evaluating the integral over $k$ in the same way as
in Sec. III.2, we obtain

$\Delta_{\infty}(z)\propto|z|^{\frac{1}{2}\frac{\gamma}{\gamma-1}}\exp\left(i\frac{\gamma-1}{A_{\gamma}}\left(\frac{|z|}{\bar{g}}\right)^{\frac{\gamma}{\gamma-1}}\left[1-\frac{1}{2}\left(\frac{\pi-\theta}{\gamma-1}\right)^{2}+i\frac{\pi-\theta}{\gamma-1}\right]\right)+(\theta\rightarrow-\theta)$
(82)

We plot $\Delta_{\infty}(z)$ from Eq. (82) in Fig. 21 for several $\gamma$
near $\gamma=2$. We see that there is an array of vortices near the real axis,
at $\rvert z\rvert>{\bar{g}}$. We now recall that Eq. (82) is valid in a range
${\bar{g}}<|z|<\omega_{cr}(\theta)$. At larger $|z|$, the oscillating term
becomes smaller than the regular piece $e^{i\theta}/|z|^{\gamma}$ and there
are no vortices. The boundary frequency for vortices is
$\omega_{cr}(\theta)\sim{\bar{g}}(|\log{(\pi-\theta)}|/(\pi-\theta))^{1/2}$
for $\gamma\leq 2$. Because $\theta\approx\pi-2\psi$ in this limit, the range
of $|z|$, where Eqn. (82) is valid, remains wide for small $\psi$. The
vortices move closer to the real axis at larger $|z|$ and for
$|z|\geq\omega_{cr}(\theta)$ escape into the lower half-plane of frequency. We
see from Fig. 21 that the number of vortices increases rapidly as $\gamma\to
2$.

###  V.2 Case $n=0$

We now show that the same vortex structure appears for $\Delta_{0}(z)$, which,
we recall, is sign-preserving along the Matsubata axis.

Figure 22: The case $n=0$. Left panel: Plot of $\log{|\Delta_{0}(z)|}$ in the
upper half plane, for $\gamma=1.75$. The gap $\Delta_{0}(z)$ for a generic $z$
in the upper half plane is obtained by analytic continuation from
$\Delta^{{}^{\prime\prime}}(\omega)$ on the real axis, Eq.(69). The two bright
points $z_{1}$ and $z_{2}$ are the locations of the vortices (points where
$|\Delta_{0}(z)|=0$). Right panel: Plot of $\eta_{0}(z)$ in the same region.

We obtain $\Delta_{0}(z)$ by Cauchy relation, Eq. (75), using as an input
$\Delta^{{}^{\prime\prime}}(\omega)$ from Eq. (66). In Fig.22 we show the
amplitude of the gap $|\Delta_{0}(z)|$ in the first quadrant of complex
$z=\omega^{\prime}+i\omega^{{}^{\prime\prime}}$ for $|z|>\bar{g}$ for
$\gamma=1.75$. We see that there are two vortices at complex $z$. We verified
that this is consistent with $m=2$ in Eq. (59) for this $\gamma$. In Fig.23 we
show that the number of vortices increases when $\gamma$ increases towards
$2$.

Figure 23: The case $n=0$. The phase of the gap function $\eta_{0}(z)$ in the
upper half plane ($z=\omega^{\prime}+i\omega^{\prime\prime}$) for different
$\gamma$ near $\gamma=2$. The locations of the vortices are marked by red
dots. At the vortex core $\Delta_{0}(z)=0$ and $\eta_{0}$ is undefined. The
number of vortices is set by $m$ in (59). This number rapidly increases as
$\gamma\to 2$.

In Fig.24

Figure 24: The gap amplitude $|\Delta_{0}(\omega)|$ as a function of $\omega$
for a set of $\gamma_{i}$. For each of these $\gamma_{i}$, the gap amplitude
vanishes at a particular $\omega$, indicating that a vortex crosses the real
axis on its way from the lower to the upper frequency half-plane.

we show how vortices emerge in the upper frequency half-plane one-by-one at a
set of $\gamma_{i}$. A given vortex is located in the lower frequency half-
plane for $\gamma<\gamma_{i}$ and moves into the upper half-plane at
$\gamma>\gamma_{i}$. Right at $\gamma=\gamma_{i}$ it appears on the real axis.
At this $\gamma$, both $\Delta^{\prime}_{0}(\omega)$ and
$\Delta^{\prime\prime}_{0}(\omega)$ vanish at some frequency, hence
$|\Delta_{0}(\omega)|=0$. We see from Fig.24 that this indeed happens for
$\gamma_{i}$ specified in this Figure. Once $\gamma$ becomes larger than
$\gamma_{i}$, $m$ increases by one, and $\delta\eta_{0}$ increases by $2\pi$.
We extended $\omega$ to $\omega\pm i\delta$ and verified that the vortex
indeed moves from the lower to the upper frequency half-plane as $\gamma$
increases through $\gamma_{i}$.

The emergence of the line of dynamical vortices (dynamical in the sense that
they are at complex $z$) can also be analyzed within a semi-phenomenological
analytical model, in which we combine oscillating behavior of
$\Delta_{0}(\omega)$ at ${\bar{g}}<\omega<\omega_{cr}$ and non-oscillating
$1/\omega^{\gamma}$ form at $\omega>\omega_{cr}$ by approximating
$\Delta_{0}(\omega)$ as the sum of the two terms:

$\Delta_{0}(\omega)=2\omega\frac{\sin{\phi^{\prime}}\cosh{\phi^{\prime\prime}}-i\cos{\phi^{\prime}}\sinh{\phi^{\prime\prime}}}{\cosh{2\phi^{\prime\prime}}-\cos{2\phi^{\prime}}}+C{\bar{g}}\left(\frac{\bar{g}}{\omega}\right)^{\gamma}e^{i\pi\gamma/2},$
(83)

where $\phi^{\prime}$ and $\phi^{\prime}$ are some functions of $\omega$. We
treat $C$ as a phenomenological parameter, which controls the width of the
range where $\Delta_{0}(\omega)$ oscillates. We assume that $C$ scales as some
power of $2-\gamma$ and vanishes at $\gamma=2$. The smaller $C$ is, the larger
is the width of this range where $\Delta_{0}(\omega)$ oscillates. The precise
forms of $\phi^{\prime}(\omega)$ and $\phi^{{}^{\prime\prime}}(\omega)$ are
not essential for this consideration, except that $\phi^{\prime}(\omega)$ must
be odd in $\omega$ and $\phi^{\prime\prime}(\omega)$ must be even, and KK
relations for $\Delta_{0}(\omega)$ must be satisfied with enough accuracy
within the range where $\Delta_{0}(\omega)$ oscillates. We use
$\phi^{\prime}(\omega)=\omega$ and
$\phi^{\prime\prime}(\omega)=(2-\gamma)|\omega|^{1.5}$ to generate the plots
in Fig.25.

Figure 25: Left panel: The plots of
$\log|\Delta^{\prime}_{0}(\omega)/\bar{g}|$ and
$\log|\Delta^{{}^{\prime\prime}}_{0}(\omega)/\bar{g}|$, from Eq. (83). The
form of $\Delta_{0}(\omega)$ is controlled by the parameter $C$. At small $C$
(which corresponds to larger $\gamma$), the gap function oscillates, and the
range of oscillations increases with decreasing $C$. Middle panel. The phase
$\eta_{0}(\omega)$, constrained to $(-\pi,\pi)$, for different $C$. The
winding number increases as $C$ decreases. Right panel. The phase
$\eta_{0}(z)$ in the upper half-plane of frequency. At large $C$, there are no
vortices in the upper half-plane. As $C$ decreases , vortices appear one-by-
one .

At large $C$, $\Delta^{\prime}_{0}(\omega)$ and
$\Delta^{{}^{\prime\prime}}_{0}(\omega)$ gradually decrease without changing
sign (the top panel). This models the behavior of the gap function for
$\gamma\sim 1$. As $C$ decreases, $\Delta^{\prime}_{0}(\omega)$ and
$\Delta^{{}^{\prime\prime}}_{0}(\omega)$ evolve and start undergoing sign
changes. At first, this does not give rise to phase winding (second panel from
the top). Then, at some critical $C=C_{1}\approx 110$, there appears a point
infinitesimally close to the real axis, where $\Delta^{\prime}_{0}(\omega)$
and $\Delta^{{}^{\prime\prime}}_{0}(\omega)$ change sign at the same
frequency. Once $C$ becomes smaller than this value, a vortex appears in the
upper frequency half-plane, and the phase $\eta_{0}(\omega)$ rapidly changes
by $2\pi$ in a narrow frequency range. At a smaller $C=C_{2}\approx 25$,
another vortex emerges on the real axis and moves to the upper frequency half-
plane at smaller $C$, causing another rapid change of $\eta_{0}(\omega)$ by
$2\pi$ (the third panel). The third vortex emerges at $C=C_{3}=1.5$ and so on
(the lower panels).

Comparing the results for $n=0$ and $n=\infty$ we see that in both cases a
line of vortices emerges at complex $z$ in the first quadrant, and equal
number of vortices forms in the second quadrant. The vortices develop for
$\gamma>1$. The number of vortices increases with $\gamma$ and tends to
infinity at $\gamma\to 2$. In this limit, the line of vortices extends to
$z=\infty$. Comparing Figs. 21 and 23, we see that the phase winding along the
real axis and the number of vortices in the upper frequency half-plane for any
given $\gamma$ are very likely the same for $n=0$ and $n=\infty$ A small
difference in the values of $\gamma_{i}$, where the new vortices appear in the
two cases, is most probably related to the fact that the $n=0$ solution is an
approximate one.

###  V.3 Finite $n$

Given the equivalence between the number of vortices at complex $z$ for $n=0$
and $n=\infty$ and near-equivalence of their position, we expect that the
geometry of the vortices near the real axis is determined solely by $\gamma$
and is the same for all values of $n$.

The emergence of the array of vortices is another indication that the pairing
at $\gamma>1$, where $V^{\prime}(\omega)$ is repulsive, is quite special. The
argument for this goes as follows. Take the set of $m$ points, where
$\Delta_{n}(z)=0$, and make it infinite by adding points $z_{i}$ at larger
$z>\omega_{cr}$, for which $\Delta_{n}(z_{i})\propto C_{n}/|z_{i}|^{\gamma}$.
As long as $\Delta_{n}(z)$ vanishes at $|z|=\infty$ everywhere in the upper
half-plane, one can continue analytically from such set and obtain a unique
$\Delta_{n}(z)$ for all $z$. We see that $\Delta_{n}(z)$ is non-zero only
because we extended the set of vortex points by adding additional $z_{i}$ with
$|z_{|}>\omega_{cr}$. The implication is that $\Delta_{n}(z)$ for all $z$ is
actually determined by the set of complex frequencies, above the crossover
scale $\omega_{cr}$. At $\gamma\to 2$, $\omega_{cr}$ tends to infinity. In
this limit, the set of the gap functions on the Matsubara axis is determined
solely by an essential singularity at $z=\infty$. This obviously points out
that the pairing at $\gamma>1$ is highly unconventional and becomes more so as
$\gamma\to 2$.

##  VI Conclusions

In this paper we continued with our analysis of the interplay between the
pairing and the non-Fermi liquid behavior in a metal for a set of quantum-
critical (QC) systems with an effective dynamical electron-electron
interaction $V(\Omega_{m})\propto 1/|\Omega_{m}|^{\gamma}$, mediated by a
critical massless boson (the $\gamma$-model). In previous papers we considered
the cases $0<\gamma<1$ and $\gamma\approx 1$. We argued that the pairing by a
gapless boson is fundamentally different from BCS/Eliashberg pairing by a
massive boson as for QC pairing there exists an infinite discrete set of
topologically distinct solutions for the gap function $\Delta_{n}(\omega_{m})$
at $T=0$ ($n=0,1,2...$), each with its own condensation energy $E_{c,n}$. Here
we extended the analysis to larger $1<\gamma<2$. We argued that the discrete
spectrum of $E_{c,n}$ get progressively denser as $\gamma$ increases towards
$2$. This increases the strength of ”longitudinal” gap fluctuations, which
tend to reduce the actual superconducting $T_{c}$ compared to the onset
temperature for the pairing. We also reported two new features on the real
axis, which again become critical at $\gamma\to 2$. First, the density of
states evolves towards a set of discrete $\delta-$functions. Second, on a real
axis, $\Delta_{n}(\omega)$ with all $n$, including $n=0$, oscillate between
${\bar{g}}$ and $\omega_{cr}>{\bar{g}}$, and the phase of the gap
$\eta_{n}(\omega)$ varies in this range by an integer number of $\pi$. We
associated this variation with the presence of dynamical vortices in the upper
half-plane of frequency, at a complex
$z=\omega^{\prime}+i\omega^{{}^{\prime\prime}}$. The vortices appear one-by-
one at a discrete set of $\gamma_{i}>1$, and the number of vortices tends to
infinity at $\gamma\to 2$. We related the emergence of oscillations on the
real axis and vortices at complex $z$ with the fact that for real $\omega$,
the interaction
$V^{\prime}(\Omega)\propto\cos(\pi\gamma/2)=/|\Omega|^{\gamma}$ becomes
repulsive for $\gamma>1$, and the imaginary part
$V^{\prime}(\Omega)\propto\sin(\pi\gamma/2)/|\Omega|^{\gamma}$ gets
progressively smaller as $\gamma\to 2$. We speculated that there is the
emergence of an infinite number of vortices and a continuum spectrum of
$E_{c,n}$ at $\gamma\to 2$ are related phenomena.

The case $\gamma=2$ requires separate consideration and will be discussed in
the next paper in the series. There we also discuss in detail longitudinal gap
fluctuations.

###### Acknowledgements.

We thank I. Aleiner, B. Altshuler, E. Berg, R. Combescot, 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, D. Maslov, F. Marsiglio, M. Metlitski, W. Metzner, A. Millis, D.
Mozyrsky, C. Pepin, 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
Y. Wang, for useful discussions. The work by YW and AVC was supported by the
NSF DMR-1834856.

##  Appendix A The exact form of $\Delta_{\infty}(\omega_{m})$

###  A.1 The Linearized gap equation

We choose the unit for frequency to be $\bar{g}$ and re-write the linearized
Eliashberg equation at zero temperature in dimensionless variable
${\bar{\omega}}=\omega/{\bar{g}}$ as

$\Delta({\bar{\omega}})=\frac{1}{2}\int
d{\bar{\omega}}^{\prime}\frac{\Delta({\bar{\omega}}^{\prime})-\Delta({\bar{\omega}})\frac{{\bar{\omega}}^{\prime}}{{\bar{\omega}}}}{\rvert{\bar{\omega}}^{\prime}\rvert}\frac{1}{\rvert{\bar{\omega}}-{\bar{\omega}}^{\prime}\rvert^{\gamma}},$
(84)

To obtain $\Delta({\bar{\omega}})$, we will explore the same strategy as in
Paper I, where we found the exact solution for $\gamma<1$. For this, we
introduce
$S({\bar{\omega}})=\text{sign}({\bar{\omega}})\left(\rvert{\bar{\omega}}\rvert-{1\over\gamma-1}\rvert{\bar{\omega}}\rvert^{1-\gamma}+i0^{+}\right)$
and

$\tilde{\Phi}({\bar{\omega}})=S({\bar{\omega}})\frac{\Delta({\bar{\omega}})}{{\bar{\omega}}}.$
(85)

The choice of the sign of $i0^{+}$ term in $S({\bar{\omega}})$ is arbitrary.
The result for $\Delta({\bar{\omega}})$ is the same for either sign.

The gap equation in terms of $\tilde{\Phi}({\bar{\omega}})$ takes the form

$\tilde{\Phi}({\bar{\omega}})=\frac{1}{2}\int
d{\bar{\omega}}^{\prime}\left(\frac{\tilde{\Phi}({\bar{\omega}}^{\prime})}{S({\bar{\omega}}^{\prime})}-\frac{\tilde{\Phi}({\bar{\omega}})}{S({\bar{\omega}})}\right)\frac{\text{sign}({\bar{\omega}}^{\prime})}{\rvert{\bar{\omega}}-{\bar{\omega}}^{\prime}\rvert^{\gamma}}-\frac{{\bar{\omega}}-S({\bar{\omega}})}{S({\bar{\omega}})}\tilde{\Phi}({\bar{\omega}}).$
(86)

Like for $\gamma<1$, we introduce complete and orthogonal basis, specified by

$\displaystyle\Phi_{\beta}(\Omega)$ $\displaystyle=$
$\displaystyle\frac{\rvert\Omega\rvert^{-2i\beta+\delta_{\Omega}}}{\rvert\Omega\rvert^{\gamma/2}},$
(87)

where $\beta\in[-\infty,\infty]$ and
$\delta_{\Omega}=\delta\mbox{sign}(1-|\Omega|)$. The functions
$\Phi_{\beta}(\Omega)$ satisfy the orthogonality and completeness relations

$\displaystyle\int_{-\infty}^{\infty}\Phi_{\beta}^{*}(\omega)\Phi_{\beta}(\Omega)\frac{d\beta}{2\pi}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\rvert\Omega\rvert^{1-\gamma}\delta(\rvert\Omega\rvert-\rvert\omega\rvert).$
(88)
$\displaystyle\int_{-\infty}^{\infty}\Phi_{\beta}^{*}(\omega)\Phi_{\beta^{\prime}}(\omega)\frac{d\omega}{|\omega|^{1-\gamma}}$
$\displaystyle=$ $\displaystyle 2\pi\delta(\beta-\beta^{\prime}).$ (89)

The function $\tilde{\Phi}({\bar{\omega}})$ is expressed via
$\Phi_{\beta}(\Omega)$ as

$\tilde{\Phi}({\bar{\omega}})=\int_{-\infty}^{\infty}\frac{d\beta}{2\pi}a_{\beta}\Phi_{\beta}({\bar{\omega}}),$
(90)

where

$a_{\beta}=\int\frac{d{\bar{\omega}}}{\rvert{\bar{\omega}}\rvert^{1-\gamma}}\Phi_{\beta}^{\ast}({\bar{\omega}})\tilde{\Phi}({\bar{\omega}}).$
(91)

Multiplying both sides of Eq. 86 by $\Phi_{\beta}^{\ast}({\bar{\omega}})$ and
integrating over $d{\bar{\omega}}/|{\bar{\omega}}|^{1-\gamma}$, we obtain the
integral equation on $a_{\beta}$ in the form

$a_{\beta}=\frac{1}{2}\int_{-\infty}^{\infty}\frac{d\beta^{\prime}}{2\pi}A_{\beta,\beta^{\prime}}a_{\beta^{\prime}}-\int_{-\infty}^{\infty}\frac{d\beta^{\prime}}{2\pi}B_{\beta,\beta^{\prime}}a_{\beta^{\prime}}\\\
$ (92)

where

$\displaystyle
A_{\beta\beta^{\prime}}=\int\frac{d{\bar{\omega}}}{|{\bar{\omega}}|^{1-\gamma}}\int
d{\bar{\omega}}^{\prime}\left(\frac{\Phi_{\beta^{\prime}}({\bar{\omega}}^{\prime})}{S({\bar{\omega}}^{\prime})}-\frac{\Phi_{\beta^{\prime}}({\bar{\omega}})}{S({\bar{\omega}})}\right)\frac{\text{sign}({\bar{\omega}}^{\prime})}{\rvert{\bar{\omega}}-{\bar{\omega}}^{\prime}\rvert^{\gamma}}\Phi_{\beta}^{\ast}({\bar{\omega}})$
(93) $\displaystyle
B_{\beta,\beta^{\prime}}=\int\frac{d{\bar{\omega}}}{|{\bar{\omega}}|^{1-\gamma}}\frac{{\bar{\omega}}-S({\bar{\omega}})}{S({\bar{\omega}})}\Phi_{\beta^{\prime}}({\bar{\omega}})\Phi_{\beta}^{\ast}({\bar{\omega}})$
(94)

One can verify that
$A_{\beta\beta^{\prime}}=(\gamma-1)B_{\beta\beta^{\prime}}F_{\beta}$, where

$\displaystyle F_{\beta}$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}dx\frac{1}{\rvert
x-1\rvert^{\gamma}}\left(|x|^{2i\beta+\gamma/2-1}-\text{sign}(x)\right)$ (95)
$\displaystyle=$
$\displaystyle\frac{2}{\gamma-1}\left(1-\epsilon_{\beta}\right),$

and

$\epsilon_{\beta}=\frac{1-\gamma}{2}\frac{\Gamma\left(\gamma/2-2i\beta\right)\Gamma\left(\gamma/2+2i\beta\right)}{\Gamma(\gamma)}\left(1+\frac{\cosh(2\pi\beta)}{\cos(\pi\gamma/2)}\right).$
(96)

Evaluating the integral for $B_{\beta\beta^{\prime}}$, we obtain
($\text{P}\int$ means Principal Value)

$\displaystyle
B_{\beta\beta^{\prime}}={2\over\gamma-1}\text{P}\int_{0}^{\infty}\frac{d{\bar{\omega}}}{{\bar{\omega}}}\frac{{\bar{\omega}}^{2i(\beta-\beta^{\prime})+\delta_{{\bar{\omega}}}}}{{(\gamma-1)}{\bar{\omega}}^{\gamma}-i0^{+}}$
(97)
$\displaystyle=\frac{2i\pi}{\gamma}\frac{e^{-2\pi(\beta-\beta^{\prime})/\gamma-2i(\beta-\beta^{\prime})/\gamma\log(\gamma-1)}}{\sinh(2\pi(\beta-\beta^{\prime})/\gamma-i0^{+})},$
(98)

The convergence factor $\delta_{{\bar{\omega}}}$ in the first line in (97) is
only relevant at small ${\bar{\omega}}$, where it is positive. Substituting
$A_{\beta\beta^{\prime}}$ and $B_{\beta\beta^{\prime}}$ into the integral
equation on $a_{\beta}$, we obtain

$a_{\beta}=-\frac{2i\pi\epsilon_{\beta}}{\gamma}\int\frac{d\beta^{\prime}}{2\pi}\frac{e^{-2\pi(\beta-\beta^{\prime})/\gamma-2i(\beta-\beta^{\prime})/\gamma\log(\gamma-1)}}{\sinh(2\pi(\beta-\beta^{\prime})/\gamma-i0^{+})}a_{\beta^{\prime}}.$
(99)

It is convenient to change the variables to $\beta=\gamma k/2$,
$a_{\beta}=e^{-\pi k-ik\log(\gamma-1)}\epsilon_{\gamma k/2}\tilde{b}_{k}$. The
gap equation then takes the form

$\tilde{b}_{k}=\frac{i}{2}\int\frac{\epsilon_{k^{\prime}\gamma/2}}{\sinh(\pi(k^{\prime}-k)+i0^{+})}\tilde{b}_{k^{\prime}}dk^{\prime},$
(100)

Once we know $\tilde{b}_{k}$, one can compute $\tilde{\Phi}({\omega})$ and the
gap function $\Delta({\omega})$ using Eqs. (90) and (85).

###  A.2 Solution of the gap equation

Let’s define the function $B_{z}$ of a complex argument $z$ as

$B_{z}=\frac{i}{2}\int\frac{\epsilon_{k^{\prime}\gamma/2}}{\sinh(\pi(k^{\prime}-z))}b_{k^{\prime}}dk^{\prime}.$
(101)

For real $k$, $\tilde{b}_{k}=B_{k-i0^{+}}$ The function $B_{z}$ satisfies the
periodicity condition

$B_{z+in}=(-1)^{n}B_{z},$ (102)

where $n$ integer, and has branch cuts at $z=x+in$, where $x$ is real. In
particular, for $n=0$,

$B_{k+i0^{+}}=\left(1-\epsilon_{k\gamma/2}\right)B_{k-i0^{+}}.$ (103)

We now take the logarithm of both sides of (103). For $\gamma>1$, the function
$\epsilon_{k\gamma/2}$ increases monotonically with $|k|$, and
$\epsilon_{k\gamma/2}=1$ at $k=k_{\gamma}=\pm\beta$. Using this, we obtain

$\log B_{k+in+i0^{+}}-\log B_{k+in-i0^{+}}=\log\rvert
1-\epsilon_{k\gamma/2}\rvert+i\pi\chi_{n}\Theta(k-\beta)+i\pi\xi_{n}\Theta(-k-\beta).$
(104)

Here $\chi_{n},\xi_{n}=\pm 1$ are two parameters, which reflect the ambiguity
of evaluating $\log{(1-\epsilon_{k\gamma/2})}$ for $\rvert k\rvert>\beta$ due
to a discontinuity of a logarithm across its branch cut. Because of the
periodicity condition (102), the l.h.s. of (104) is independent to $n$, hence
$\chi_{n}\equiv\chi$ and $\xi_{n}\equiv\xi$ are just two numbers, each is
either $+1$ or $-1$. Using the Sokhotski-Plemelj theorem, we can then re-
express (104) as

$\displaystyle\log B_{k}$ $\displaystyle=$
$\displaystyle\frac{1}{2i}\int_{-\infty}^{\infty}\log\rvert
1-\epsilon_{k\gamma/2}\rvert\coth\left(\pi(k^{\prime}-k)\right)dk^{\prime}+\frac{\pi\chi}{2}\int_{\beta}^{\infty}\coth\left(\pi(k^{\prime}-k)\right)dk^{\prime}$
(105)
$\displaystyle+\frac{\pi\xi}{2}\int_{-\infty}^{-\beta}\coth\left(\pi(k^{\prime}-k)\right)dk^{\prime}+G_{k},$

where $G_{k}$ is an analytic function, which satisfies the periodicity
relation

$\displaystyle G_{k+in}$ $\displaystyle=$ $\displaystyle
G_{k}+(2l+1)in\pi,~{}~{}l,n\in\mathbb{Z}.$ (106)

The last condition uniquely specifies $G_{k}=c+(2l+1)\pi k$, where $c$ is a
constant. The latter is irrelevant as one can easily verify that it only
contributes to an irrelevant overall factor in
$\Delta_{\infty}({\bar{\omega}})$. Using (106) and absorbing the divergent
piece into $c$, we obtain from (105):

$\displaystyle\tilde{b}_{k}=B_{k-i0^{+}}$ $\displaystyle=$
$\displaystyle\exp\Bigg{[}\frac{1}{2i}\int_{-\infty}^{\infty}\log\rvert
1-\epsilon_{k\gamma/2}\rvert\coth\left[\pi(k^{\prime}-k+i0^{+})\right]dk^{\prime}$
(107)
$\displaystyle-\frac{\chi}{2}\log\frac{\sinh[\pi(k_{\gamma}-k+i0^{+})]}{\sinh(\pi
k_{\gamma})}+\frac{\xi}{2}\log\frac{\sinh[\pi(k_{\gamma}+k-i0^{+})]}{\sinh(\pi
k_{\gamma})}$ $\displaystyle+(2l+1-\frac{\chi+\xi}{2})\pi k\Bigg{]}.$

This $\tilde{b}_{k}$ depends on three parameters: $\chi,\xi=\pm 1$ and
$l\in\mathbb{Z}$. Below we select $\chi,\xi$, and $l$ based on the two
conditions: (1) the integrand in the r.h.s. of the gap equation (100) must be
convergent; (2) $b_{k}$ must give rise to a gap function
$\Delta({\bar{\omega}})$, which vanishes at ${\bar{\omega}}\to\infty$. The
first condition puts the following constraints:

  * •

If $\chi=\xi=1$, then $l=0$.

  * •

If $\chi=\xi=-1$, then $l=-1$.

  * •

If $\chi=-\xi=1$, then $l=0,-1$.

  * •

If $\chi=-\xi=-1$, there is no solution for $l$.

We verified that the second condition is satisfied only if $\chi=-\xi=1$ and
$l=0,-1$. The dependence on $l$ in (106), and the difference between the two
choices for $l$ is an irrelevant overall factor for $\Delta({\bar{\omega}})$.
We then concludes that the two conditions uniquely specify ${\tilde{b}}_{k}$.

Substituting these $\chi,\xi$, and $2l+1$ into (106) and re-expressing the
result back in terms of $b_{k}=e^{-\pi k-ik\log(\gamma-1)}{\tilde{b}}_{k}$, we
obtain

$\displaystyle b_{k}$ $\displaystyle=$
$\displaystyle\exp\left[-ik\log(\gamma-1)\right]\exp\Bigg{[}\frac{1}{2i}\int_{-\infty}^{\infty}\log\rvert
1-\epsilon_{k^{\prime}\gamma/2}\rvert\coth[\pi(k^{\prime}-k+i0^{+})]dk^{\prime}$
(108)
$\displaystyle-\frac{1}{2}\log\frac{\sinh[\pi(k_{\gamma}-k+i0^{+})]}{\sinh[\pi
k_{\gamma}]}-\frac{1}{2}\log\frac{\sinh[\pi(k_{\gamma}+k-i0^{+})]}{\sinh[\pi
k_{\gamma}]}\Bigg{]}.$

###  A.3 Computation of $\Delta_{\infty}({\omega})$

We now substitute this $b_{k}$ into Eq. (90) and compute
$\tilde{\Phi}({\omega})$ and the gap function $\Delta({\omega})$ using Eq.
(85). It is convenient to introduce
$y\equiv\rvert{\bar{\omega}}\rvert^{\gamma}$. In terms of $y$, we obtain

$\displaystyle\Delta_{\infty}(y)=y^{1/2}\int_{-\infty}^{\infty}dkb_{k}e^{-ik\log
y}.$ (109)

The $y^{1/2}$ in the prefactor can be eliminated by shifting the integration
contour away from the real axis, to $k-i/2$. This is a safe procedure because
$b_{k}$ is analytic within the interval $-1<\text{Im}(k)<0$. Shifting the
integration, we obtain the final expression

$\Delta_{\infty}(y)=\int_{-\infty}^{\infty}dkb_{k}e^{-ik\log y},$ (110)

where

$b_{k}=\frac{e^{-iI_{k}-ik\log(\gamma-1)}}{\sqrt{\cosh(\pi(k-\beta))\cosh(\pi(k+\beta))}}$
(111)

and

$I_{k}=\frac{1}{2}\int_{-\infty}^{\infty}dk^{\prime}\log{|\epsilon_{k^{\prime}}-1|}\tanh{\pi(k^{\prime}-k)}.$
(112)

###  A.4 Series expansion for $\Delta_{\infty}(y)$

The function $I_{k}$ can be expressed as an infinite product of the Gamma-
functions (see Paper I for details). For $b_{k}$, this yields, up to an
overall factor,

$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+\frac{2m}{\gamma}-ik\right)}{\Gamma\left(\frac{1}{2}-i(k+i\beta_{m})\right)\Gamma\left(1+\frac{2m}{\gamma}+ik\right)}$
(113)

Here $\beta_{m}>0$ are the solutions of $\epsilon_{i\beta_{m}}=1$ for
imaginary argument. There is an infinite set of such $\beta_{m}$, specified by
integer $m=0,1,2..$ and located at $1/2+2m/\gamma<b_{m}<1/2+2(m+1)/\gamma$
(see Fig. 3 (b) in the main text).

The integral in (110) can be evaluated by closing the integration contour
along an infinite arc in the complex plane of frequency. For $y<1$, i.e.,
$\omega_{m}<{\bar{g}}$, the arc must be in the upper half-plane, and for
$y>1$, i.e., $\omega_{m}>{\bar{g}}$, in the lower half-plane. Viewed as a
function of complex $k$, $b_{k}$ has poles from individual $\Gamma$-functions
in the upper frequency half-plane at $k=\pm\beta+i(n+1/2)$, where $n=0,1,2..$
and at $k=i\beta_{m}+i(n+1/2)$, and in the lower half-plane, at $k=-i(n+1)$
and $k=-i(1+2m/\gamma+n)$, where $n=0,1,2,...$ and $m=1,2,...$.

####  A.4.1 $y<1$

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

$\Delta_{\infty}(y)=\text{Re}~{}\sum_{n=0}^{\infty}e^{i(\beta\log{y}+\phi)}C^{<}_{n}~{}y^{n+1/2}+\sum_{n,m=0}^{\infty}D^{<}_{n,m}y^{n+b_{m}/\gamma+1/2}$
(114)

We cited this result in Eq. (13) in the main text. Here $\phi$ is some
particular, $\gamma-$dependent number, and $C^{<}_{n}$ and $D^{<}_{n,m}$ are
$\gamma$-dependent coefficients. The leading term in (114) at small $y$ is

$\Delta_{\infty}(y)=y^{1/2}C^{<}_{n}\cos{\left(\beta\log{y}+\phi\right)}$
(115)

The first subleading term scales as $y^{3/2}$ for $\gamma\leq 1.81$, where
$\beta_{0}/\gamma>1$, and as $y^{1/2+\beta_{0}/\gamma}$ for $\gamma$ closer to
$2$, where $\beta_{0}/\gamma<1$. We verified this explicitly by subtracting
the leading term, given by (115), from the exact solution and identifying the
leading term in the leftover (see Fig. 26).

Figure 26: The expansion of the exact $\Delta_{\infty}(\omega_{m})$ in powers
of $\omega_{m}$. The leading term at small $\omega_{m}$ is
$\Delta^{{\text{l}ead}}_{\infty}|\omega_{m}|^{\gamma/2}\cos{(\beta\gamma\log{|\omega_{m}|/{\bar{g}}}+\phi)}$.
The subleading terms form series of local and non-local terms. Local series
hold in powers of $|\omega_{m}|^{n\gamma}$ ($n=1,2...$), and all terms in
these series oscillate, like the leading term. The non-local terms do not
oscillate and form series in $|\omega_{m}|^{(n^{\prime}+\beta_{m})\gamma}$,
where $n^{\prime}=(0,1,2..)$ and $\beta_{m}$ ($m=0,1,2..$) are some
$\gamma-$dependent numbers. We show
$|\Delta^{\infty}-\Delta^{{\text{l}ead}}_{\infty}|$ for $\gamma=1.01$ and
$\gamma=1.99$. An analytical analysis shows that in the first case,
$\beta_{0}>1$, such that the leading correction comes from the local series
and oscillates. For $\gamma=1.99$, $\beta_{0}\simeq 0.9<1$. In this case the
leading correction comes from non-local series and does not oscillate. The
numerical evaluation of $|\Delta^{\infty}-\Delta^{{\text{l}ead}}_{\infty}|$
shown here, confirms these results.

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

$\displaystyle\Delta_{\infty,L}(y)=\text{
Re}~{}\sum_{n=0}^{\infty}e^{i(\beta\log{y}+\phi)}C^{<}_{n}~{}y^{n+1/2}$ (116)
$\displaystyle\Delta_{\infty,NL}(y)=\sum_{n,m=0}^{\infty}D^{<}_{n,m}y^{n+b_{m}/\gamma+1/2}$
(117)

The coefficients $C^{<}_{n}$ in (116) can be obtained analytically, as we
already found in Papers I and III for $\gamma<1$. The computations for
$\gamma>1$ are similar, and we only present the result. We obtained

$C^{<}_{n}=C^{<}_{0}\displaystyle\prod_{m=1}^{n}\frac{1}{{\bar{I}}_{m}},$
(118)
$\displaystyle{\bar{I}}_{m}=\frac{1}{2}\left[\frac{\Gamma((m+1/2)\gamma+i{\beta}\gamma)\Gamma((1/2-m)\gamma-i{\beta}\gamma)}{\Gamma(\gamma)}-\frac{\Gamma(\gamma(1/2+i{\beta}))\Gamma(\gamma(1/2-i\beta))}{\Gamma(\gamma)}\right]+$
$\displaystyle\frac{\Gamma(1-\gamma)}{2}\left(\frac{\Gamma((m+1/2)\gamma+i{\beta}\gamma)}{\Gamma(1-(1/2-m)\gamma+i{\beta}\gamma)}+\frac{\Gamma((1/2-m)\gamma-i{\beta}\gamma)}{\Gamma(1-(m+1/2)\gamma-i{\beta}\gamma)}\right)$
$\displaystyle-\frac{\Gamma(1-\gamma)}{2}\left(\frac{\Gamma(\gamma(1/2+i{\beta}))}{\Gamma(1-\gamma(1/2-i{\beta}))}+\frac{\Gamma(\gamma(1/2-i{\beta}))}{\Gamma(1-\gamma(1/2+i\beta))}\right)$
(119)

Using the expansion for $\Gamma$-functions, one can verify that at large $m$,
${\bar{I}}_{m}\sim m^{\gamma-1}$, i.e., $C^{<}_{n}\sim 1/n^{n(\gamma-1)}$.
Interestingly, the sum in (116) then converges for any $y$ and can be obtained
numerically by summing up the proper number of terms in (116). We show the
results for $\Delta_{\infty,L}(y)$ for different $\gamma$ in Fig. (27),
extending also into the range $y>1$.

Figure 27: Panels (a)-(g). The gap function $\Delta_{\infty,L}(y)$ as a
function of $y=(|\omega_{m}|/{\bar{g}})^{\gamma}$ for various $\gamma$. At
$y<1$, $\Delta_{\infty,L}(y)$ oscillates for all $\gamma$ with the period set
by $\log y$. As $\gamma$ increases towards $2$, the new oscillating regime
emerges if one extends $\Delta_{\infty,L}(y)$ to $y>1$. In this new regime the
period of oscillations is set by $y$. (h) The analytical form of
$\Delta_{\infty,L}(y)$ at $\gamma\to 2$. (i) The approximately linear
dependence of the upper boundary for new oscillations, $y^{*}$, on $2-\gamma$.

At $\gamma\leq 2$, one can expand the $\Gamma$ functions in (119) in
$2-\gamma$. For ${\bar{I}}_{m}$ we then obtain, to first order in $2-\gamma$,

${\bar{I}}_{m}=\frac{im}{\beta}\left(1+iQ(2-\gamma)(m+1)\right)$ (120)

where $Q\approx 0.7$. This holds for $Q(2-\gamma)m<1$. For larger $m$,
$Q(2-\gamma)m$ becomes $\sin(Q(2-\gamma)m)$. Substituting ${\bar{I}}_{m}$ from
(120) into (118), evaluating the product to first order in $(2-\gamma)$, and
substituting the result into (116), we obtain, to the same accuracy,

$\Delta_{\infty,L}(y)\propto\sqrt{y}e^{-2Q(2-\gamma)\beta y}\cos{(\beta(\log
y-y)-Q(2-\gamma)(y\beta)^{2}/2)+\phi)}$ (121)

The expression becomes particularly simple for $\gamma=2$, where

$\Delta_{\infty,L}(y)\propto\sqrt{y}\cos{(\beta(\log y-y)+\phi)}$ (122)

####  A.4.2 $y>1$

For $y>1$, relevant poles are in the lower half-plane, at $k=-i(n+1)$ and
$k=-i(n+1+2m/\gamma)$. This yields the expansion for $\Delta_{\infty}(y)$ in
powers of $1/y$ in the form

$\Delta_{\infty}(y)=\sum_{n=0}^{\infty}C^{>}_{n}~{}\left(\frac{1}{y}\right)^{n+1}+\sum_{n,m=0}^{\infty}D^{>}_{n,m}\left(\frac{1}{y}\right)^{n+1+2(m+1)/\gamma}$
(123)

We cited this result in Eq. (14) in the main text. The leading term in the
series is $1/y$, i.e., at $\omega_{m}\gg{\bar{g}}$,
$\Delta_{\infty}(\omega)\propto 1/|\omega|^{\gamma}$. The exact
$\Delta_{\infty}(\omega)$ clearly shows this behavior (see Fig. 4)

Eq. (123) is formally the same as the series expansion result for $\gamma<1$,
however the coefficients $C^{>}_{n}$ and $D^{>}_{n,m}$ depend on $\gamma$. We
now argue that this dependence is qualitatively different for $\gamma<1$ and
$\gamma>1$. Namely, we argue below that for $\gamma>1$, there is a universal
piece in $\Delta_{\infty}(y)$ at large $y$, which does not fit into the power-
law series in (123). We obtain this piece by analyzing the form of $b_{k}$ in
Eq. (110) and evaluating the integral over $k$ directly, without closing the
itegation contour in the complex plane of $k$. While we didn’t compute
$C^{>}_{n}$ and $D^{>}_{n,m}$ explicitly, the presence of a piece that does
not fit into power-law series implies that the series expansion in (123) does
not converge at $\gamma>1$ (e.g., the coefficients $C^{>}_{n}$ grow fast
enough with $n$ and eventually overcome the smallness of $1/y^{n}$) and would
yields incorrect results, starting from some critical $n$, which depends on
$y$. In Paper V, where we specifically focus on $\gamma=2$, we show explicitly
that this is the case and determine a critical $n$.

Below we focus on the universal, non-power-law term in $\Delta_{\infty}(y)$.
We show that this term is present at $\gamma>1$, because the complex phase of
the integrand in Eq. (110) has an extreme at $k=k_{*}\sim y^{1/(\gamma-1)}$.
The universal contribution then comes from the expansion of the integrand
around $k=k^{*}$. There is no such extreme for $\gamma<1$.

To identify the universal term, we consider large $y$ and analyze the
contribution to the integral in (110) from $k$ above some $k_{min}={\cal
O}(1)$. We label the corresponding term as $\Delta_{\infty}^{u}(y)$. In
explicit form,

$\displaystyle\Delta_{\infty;u}(y)$
$\displaystyle=\int_{k_{\text{min}}}^{\infty}dk{\cos{(I_{k}+k\log
y(\gamma-1))}\over(\cosh\left[\pi(k-\beta)\right]\cosh\left[\pi(k+\beta)\right])^{1/2}}.$
(124)

For large $k$, the integral, which determines $I_{k}$ in Eq. (112), is
determined by $k^{\prime}\sim k$. This integral contains ${\epsilon}_{\gamma
k^{\prime}/2}$. For $\gamma>1$, ${\epsilon}_{\gamma k^{\prime}/2}$ is an
increasing function of $k^{\prime}$ and its leading term at large $k^{\prime}$
is ${\epsilon}_{\gamma
k^{\prime}/2}\simeq(\gamma-1)(A_{\gamma}k^{\prime})^{\gamma-1}$, where
$A_{\gamma}=\gamma\left(\frac{\pi}{2\Gamma(\gamma)\cos(\pi(2-\gamma)/2)}\right)^{1/(\gamma-1)}$.
Substituting this form into Eq. (112), we obtain

$\displaystyle
I(k)+\log{[y(\gamma-1)]}\simeq(\gamma-1)k\ln\frac{ye}{A_{\gamma}k},$ (125)

Substituting this form into (124), we find that the argument of the integrand
has a maximum at

$k_{*}=y^{\frac{1}{\gamma-1}}A_{\gamma}^{-1}.$ (126)

Expanding $I(k)+k\log\left[y(\gamma-1)\right]$ to quadratic order in
$u=k/k_{*}-1$, as

$I(k)+k\log\left[y(\gamma-1)\right]\approx(\gamma-1)k_{*}\left(1-\frac{1}{2}u^{2}\right),$
(127)

substituting into (124), and approximating the denominator in (124) by its
large $k$ value, we obtain

$\Delta_{\infty}^{u}(y)\sim\sqrt{k_{*}}e^{-\pi
k^{*}}\text{Re}\left[e^{-i(\gamma-1)k_{*}\left(1-\frac{\pi^{2}}{(\gamma-1)^{2}}\right)}\int_{-[k^{*}(\gamma-1)/2]^{1/2}[1-i\pi/(\gamma-1)]}^{\infty}d{\tilde{u}}e^{i{\tilde{u}}^{2}}\right].$
(128)

where ${\tilde{u}}=u((\gamma-1)k^{*}/2)^{1/2}$. The universal part of this
expression is obtained by taking the lower limit to $-\infty$, in which case
$\int_{-\infty}^{\infty}e^{i{\tilde{u}}^{2}}=\sqrt{\pi}e^{i\pi/4}$.
Substituting into (128) and expressing $k^{*}$ in terms of
$|\omega_{m}|/{\bar{g}}$, we obtain for the universal term

$\displaystyle\Delta_{\infty}^{u}(\omega_{m})\sim$
$\displaystyle\left({\rvert\omega_{m}\rvert\over{\bar{g}}}\right)^{{\gamma\over
2(\gamma-1)}}\exp\left[-\frac{\pi}{A_{\gamma}}\left({\rvert\omega_{m}\rvert\over{\bar{g}}}\right)^{{\gamma\over\gamma-1}}\right]$
$\displaystyle\cos\left[\frac{\gamma-1}{A_{\gamma}}\left(1-\frac{\pi^{2}}{(\gamma-1)^{2}}\right)\left({\rvert\omega_{m}\rvert\over{\bar{g}}}\right)^{{\gamma\over\gamma-1}}-\frac{\pi}{4}\right].$
(129)

Working along the same lines, but extending the analysis to include the
subleading terms at large $k$, we obtain $\Delta_{\infty}^{u}(\omega_{m})$
with a more accurate argument of the $\cos$:

$\Delta_{\infty}^{u}(\omega_{m})\sim\left({\rvert\omega_{m}\rvert\over{\bar{g}}}\right)^{{\gamma\over
2(\gamma-1)}}\exp\left({-\frac{\pi}{A_{\gamma}}\left({\rvert\omega_{m}\rvert\over{\bar{g}}}\right)^{{\gamma\over\gamma-1}}}\right)\cos{\Psi_{m}(\omega_{m}/{\bar{g}})}$
(130)

where

$\Psi_{m}(x)=\frac{\gamma-1}{A_{\gamma}}\left[1-{1\over
2}\left({\pi\over\gamma-1}\right)^{2}\right]x^{{\gamma\over\gamma-1}}+\frac{\pi^{2}}{2(\gamma-1)A_{\gamma}}x^{\frac{\gamma(2-\gamma)}{\gamma-1}}+\frac{x^{{\gamma(2-\gamma)\over\gamma-1}}-1}{2-\gamma}A_{\gamma}^{1-\gamma}-\frac{\pi}{4}.$
(131)

At $\gamma\to 2$, $\Psi(x)$ acquires a simple form
$\Psi(x)=(x^{2}(1-\pi^{2}/2)+\log{x^{2}})/\pi+\pi/4$.

###  A.5 Analytical continuation

The gap function given by (110) can be analytically continued away from the
Matsubara axis by a simple rotation of the argument: $i\omega_{m}\to
z=|z|e^{i\psi}$. Under this transformation, $\log y$ in (110) transforms into
$\log{y_{z}}-i\theta$, where $y_{z}=(\rvert z\rvert/{\bar{g}})^{\gamma}$ and
$\theta=(\pi/2-\psi)\gamma$. The gap function transforms to

$\Delta_{\infty}(z)=\int_{-\infty}^{\infty}dk\frac{e^{-\theta
k}e^{-iI_{k}-ik\log y_{z}}}{\sqrt{\cosh(\pi(k-\beta))\cosh(\pi(k+\beta))}}.$
(132)

Applying this transformation, we extend (131) to

$\Delta_{\infty}^{u}(z)\sim\left({\rvert\omega_{m}\rvert\over{\bar{g}}}\right)^{{\gamma\over
2(\gamma-1)}}\left(Q(\theta,|z|)+Q(-\theta,|z|)\right)$ (133)

where

$Q(\theta,|z|)=\exp\left({-\frac{\pi-\theta}{A_{\gamma}}\left({\rvert
z\rvert\over{\bar{g}}}\right)^{{\gamma\over\gamma-1}}}\right)e^{i\Psi(|z|/{\bar{g}},\theta)}$
(134)

and

$\displaystyle\Psi(x,\theta)=$
$\displaystyle\frac{\gamma-1}{A_{\gamma}}\left(1-{1\over