Interplay between superconductivity and non-Fermi liquid at a quantum-critical point in a metal.
I: The $\gamma$ model and its phase diagram at $T=0$. The case $0<\gamma<1$.

The interplay between superconductivity and pairing near a quantum-critical point (QCP) in a metal is a fascinating subject, which attracted substantial attention in the correlated electron community after the discovery of superconductivity (SC) in the cuprates, heavy-fermion and organic materials, and, more recently, Fe-pnictides and Fe-chalcogenides. (see, e.g., Refs. Monthoux et al. (2007); Scalapino (2012); Norman (2014); Maiti and Chubukov (2014); Keimer et al. (2015); Shibauchi et al. (2014); Fernandes and Chubukov (2016); Fradkin et al. (2010); Yang et al. (2006); Fratino et al. (2016); Sachdev (2018); Coleman (2015)). Itinerant QC models, analyzed analytically in recent years include, e.g., models of fermions in spatial dimensions $D\leq 3$ (Refs. Mross et al. (2010); Metlitski et al. (2015); Mahajan et al. (2013); *raghu2; *raghu3; *raghu4; *raghu5; Wang et al. (2018, 2017a); Fitzpatrick et al. (2015b)), two-dimensional (2D) models near a spin-density-wave (SDW) and charge-density-wave (CDW) instabilities (Refs.  Millis (1992); Altshuler et al. (1995a); Sachdev et al. (1995); Abanov et al. (2001a, 2003); *acs2; *finger_2001; *acn; Sachdev et al. (2009); *Subir2; *Sachdev2019; Vojta and Sachdev (1999); Wang and Chubukov (2013a); Efetov et al. (2013); *efetov2; *efetov3; Tsvelik (2017); Metlitski and Sachdev (2010a); Dalidovich and Lee (2013); *sslee_2018; Bauer and Sachdev (2015); Castellani et al. (1995); *ital2; *ital3; Wang and Chubukov (2013a); *wang23; *wang_22; Khodas et al. (2010); *tsvelik; *tsvelik_1; *tsvelik_2; *tsvelik_3; Metzner et al. (2003); *DellAnna2006; Vilardi et al. (2018); *metzner_1; *metzner_2; Sordi et al. (2012); Chubukov and Wölfle (2014)), $2k_{F}$ density-wave instability (Refs.Altshuler et al. (1995b); *2kf2; *2kf3), pair-density-wave instability  Kim et al. (2008); Lawler and Fradkin (2007); *Fradkin2009; *Fradkin2014, $q=0$ instabilities towards a Pomeranchuk order instabilityBonesteel et al. (1996); *Oganesyan2001; *q=01; *q=02; Metzner et al. (2003); *DellAnna2006; Rech et al. (2006); Lee (2009); Metlitski and Sachdev (2010b); Lederer et al. (2015); Maslov and Chubukov (2012); *maslov2; *maslov3; *maslov4; You et al. (2016); Wang et al. (2001); *triplet2; *triplet3; Lederer et al. (2017); Klein et al. (2019); *avi_2 and towards circulating currents Bok et al. (2016), 2D fermions at a half-filled Landau level Lee (1989); *Monien1993; *Nayak1994; *Altshuler1994; *Kim_1994, generic QC models with different critical exponents Raghu et al. (2015); Moon and Chubukov (2010); Khveshchenko and Shively (2006); She and Zaanen (2009); *She2010; Wang et al. (2016); Lee et al. (2018); Wu et al. (2019a, b); *Abanov_19; Yuzbashyan et al. (a), Sachdev-Ye-Kitaev (SYK) and SYK-Yukawa models Patel and Sachdev (2019); Esterlis and Schmalian (2019); Wang (2020); Hauck et al. (2019); Chowdhury and Berg (2020), strong coupling limit of electron-phonon superconductivity, Combescot (1995); Bergmann and Rainer (1973); *Bergmann2; *ad; Marsiglio et al. (1988); *Marsiglio_91; Karakozov et al. (1991); Chubukov et al. (2020a), and even color superconductivity of quarks, mediated by gluon exchangeSon (1999); *son2. These problems have also been studied using various numerical techniques Yang et al. (2011); Gerlach et al. (2017); *berg_2; *berg_3; *berg_4; Haule and Kotliar (2007); *kotliar2; Sordi et al. (2012); Georges et al. (2013); *georges2.

From theory perspective, the key interest in the pairing near a QCP is due to the fact that an effective dynamic electron-electron interaction, $V(q,\Omega)$, mediated by a critical collective boson, which condenses at a QCP, provides strong attraction in one or more pairing channels and, at the same time, gives rise to non-Fermi liquid (NFL) behavior in the normal state. The two tendencies compete with each other: fermionic incoherence, associated with NFL behavior, destroys Cooper logarithm and reduces the tendency to pairing, while the opening of a SC gap eliminates the scattering at low energies and reduces the tendency to NFL. To find the winner (SC or NFL), one needs to analyze the set of integral equations for the fermionic self-energy, $\Sigma({\bf k},\omega)$, and the gap function, $\Delta({\bf k},\omega)$, for fermions with momentum/frequency $({\bf k},\omega)$ and $(-{\bf k},-\omega)$.

We consider the subset of models in which collective bosons are slow modes compared to dressed fermions, for one reason or the other. In this situation, which bears parallels with Eliashberg theory for electron-phonon interaction Eliashberg (1960), the self-energy and the pairing vertex can be approximated by their values at the Fermi surface (FS) and computed within one-loop approximation. The self-energy on the FS, $\Sigma({\bf k},\omega)$ is invariant under rotations from the point group of the underlying lattice. The rotational symmetry of the gap function $\Delta({\bf k}_{F},\omega)$ and the relation between the phases of $\Delta({\bf k}_{F},\omega)$ on different FS’s in multi-band systems are model specific. Near a ferromagnetic QCP, the pairing interaction mediated by a soft boson is attractive in the p-wave channel. Near an antiferromagnetic QCP, the strongest attraction is in $d-$wave channel for the case of a single FS and fermionic density near half-filling. For a nearly compensated metal with hole and electron pockets, as in Fe-based superconductors, the two attractive channels near an antiferromagnetic QCP are $s^{+-}$ and $d-$wave. Near a $q=0$ nematic QCP, the pairing interaction, mediated by soft nematic fluctuations, is attractive in all channels: $s-$wave, $p-$wave, $d-$wave, etc. In each particular case, one has to project the pairing interaction into the proper irreducible channel, find the strongest one, and solve for the pairing vertex within a given pairing symmetry. In principle, even after projection, one has to solve an infinite set of coupled equations in momentum space as in a lattice system each irreducible representation contains an infinite set of eigenfunctions. However, near a QCP the pairing is often confined to a narrow range on the FS around special ”hot spots”. In this situation, the momentum integration near the FS can be carried out exactly, and the set of coupled equations for the self-energy and the gap function reduce to two one-dimensional (1D) equations for $\Sigma(\omega)$ and $\Delta(\omega)$, each with frequency-dependent effective ’local” interaction $V(\Omega)$ (Ref. Abanov et al. (2001a). The same holds for the cases of s-wave pairing by a soft optical phonon, when momentum-dependencies of $\Sigma$ and $\Delta$ are not crucial and can be neglected, of non-s-wave pairing, when one eigenfunction gives the dominant contribution to the gap (e.g., $\cos k_{x}-\cos k_{y}$ for $d-$wave pairing in the cuprates), and for the case when the fermionic density of states is peaked at particular ${\bf k}_{F}$ due to, e.g. van-Hove singularities ( see e.g., Maiti and Chubukov, 2014 and references therein).

Away from a QCP, the effective $V(\Omega)$ tends to a finite value at $\Omega=0$. In this situation, the fermionic self-energy has a FL form at the smallest frequencies, the equation for $\Delta(\omega)$ is similar to that in a conventional Eliashberg theory for phonon-mediated superconductivity, and the only qualitative distinction for electronically-mediated pairing is that $V(\Omega)$ by itself changes below $T_{c}$ due to feedback from fermionic pairing on collective modes.

At a QCP, the situation becomes qualitatively different because the effective interaction $V(\Omega)$, mediated by a critical massless boson, becomes a singular function of frequency: on Matsubara axis $V(\Omega_{m})\propto 1/|\Omega_{m}|^{\gamma}$ (Fig. 1). The exponent $\gamma>0$ depends on the model, ranging from small $\gamma=0(\epsilon)$ in models in $D=3-\epsilon$ to $\gamma\leq 1$ in 2D models at SDW, CDW, and nematic QCP and in Yukawa-SYK model Chubukov et al. (2020b); Maiti and Chubukov (2014); Chubukov et al. (2020a); Esterlis and Schmalian (2019); Wang (2020). The case $\gamma=2$ corresponds to fermions, interacting with a critical Einstein phonon. The set of models with $V(\Omega_{m})\propto 1/|\Omega_{m}|^{\gamma}$ has been nicknamed the $\gamma$-model, and we will use this notation. We present an (incomplete) set of models in Tables 1- 3.

In this and subsequent papers, we present comprehensive analysis of the competition between NFL and SC within the $\gamma-$model. We define the dimensionless $V(\Omega_{m})$ as $V(\Omega_{m})=({\bar{g}}/|\Omega_{m}|)^{\gamma}$, where ${\bar{g}}$ is the effective fermion-boson coupling. This ${\bar{g}}$ is the only parameter in the model with the dimension of energy, and we will see that it determines both the magnitude of the gap function $\Delta(\omega_{m})$ and the upper limit of NFL behavior in the normal state, i.e., the scale below which $\Sigma(\omega_{m})>\omega_{m}$. We show that the competition between NFL and SC holds for all values of $\gamma$, but the physics and the computational analysis are different for the models with $\gamma<1$, $\gamma=1$, $1<\gamma<2$, $\gamma=2$, and $\gamma>2$, which we consider separately. We show that for all $\gamma$, a NFL self-energy in the normal state does not prevent the formation of bound pairs of fermions at a non-zero onset temperature $T_{p}$. However, we argue that $T_{p}$ is generally higher than the actual superconducting $T_{c}$, and at $T_{c}<T<T_{p}$ bound pairs remain incoherent. In this $T$ range various observables, e.g., the spectral function and the density of states, display pseudogap behavior. This holds for all $\gamma$, but the width of the pseudogap phase increases with $\gamma$. We relate the pseudogap behavior to strong gap fluctuations, but argue that these fluctuations can be viewed as long-wavelength phase fluctuations only in some range ${\bar{T}}_{p}<T<T_{p}$. At smaller $T_{c}<T<{\bar{T}}_{p}$, pseudogap behavior is predominantly due to the existence of low-energy ’longitudinal” gap fluctuations, which change the functional form of $\Delta(\omega)$ and cause phase slips. We argue that longitudinal gap fluctuations develop a zero mode at $\gamma=2$ in which case SC order gets destroyed already at $T=0$. Specifically, the presence of a zero mode implies that there exists an infinite number of solutions for $\Delta(\omega)$, all with the same condensation energy. In this situation, the ground state is a mixture of different $\Delta(\omega)$, each with its own phase. We will argue that $T_{c}$ gradually vanishes at $\gamma\to 2$ and remains zero at larger $\gamma$, while $T_{p}$ and ${\bar{T}}_{p}$ stay finite. We show the phase diagram in Fig. 2. Our results present the scenario for the pseudogap, which holds even if for a given solution for $\Delta(\omega)$ SC stiffness is larger than $T_{c}$, i.e., conventional phase fluctuations are weak. This scenario is complementary to the one in which the smallness of the superfluid stiffness is caused by the closeness to a Mott transition (see, e.g., Ref. Simard et al. (2019) and references therein).

We consider itinerant fermions at the onset of a long-range order in either spin or charge channel. At the critical point, the propagator of a soft boson becomes massless and mediates singular interaction between fermions. We follow earlier works Abanov et al. (2001a, 2003); Moon and Chubukov (2010); Metlitski and Sachdev (2010b); Mross et al. (2010); Monthoux et al. (2007); Efetov et al. (2013); Metlitski et al. (2015); Raghu et al. (2015); Haslinger and Chubukov (2003); Wang et al. (2016); Lee et al. (2018) and 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. This is the case for a conventional phonon-mediated superconductivity, where for $q\sim k_{F}$ a typical fermionic frequency is of order $E_{F}$, while typical bosonic frequency is of order Debye frequency $\omega_{D}$. The ratio $\delta_{E}=\omega_{D}/E_{F}$ is the small parameter for Eliashberg theory of phonon-mediated superconductivity. This theory allows one to obtain a set of coupled integral equations for frequency dependent fermionic self-energy and the pairing vertex. By analogy, the theory of electronic superconductivity, mediated by soft collective bosonic excitations in spin or charge channel, is also often called Eliashberg theory. We will use this convention.

Within the Eliashberg approximation, one can explicitly integrate over the momentum component perpendicular to the Fermi surface (for a given pairing symmetry) and reduce the pairing problem to a set of coupled integral equations for frequency dependent self-energy $\Sigma(\omega_{m})$ and the pairing vertex $\Phi(\omega_{m})$ with effective frequency-dependent dimensionless interaction $V(\Omega)=({\bar{g}}/|\Omega|)^{\gamma}$. This interaction gives rise to NFL form of the self-energy in the normal state and, simultaneously, gives rise to the pairing.

At $T=0$, the coupled Eliashberg equations for the pairing vertex and the fermionic self-energy are, in Matsubara formalism,

	$\displaystyle\Phi(\omega_{m})$	$\displaystyle=$	$\displaystyle\frac{{\bar{g}}^{\gamma}}{2}\int d\omega^{\prime}_{m}\frac{\Phi(\omega^{\prime}_{m})}{\sqrt{{\tilde{\Sigma}}^{2}(\omega^{\prime}_{m})+\Phi^{2}(\omega^{\prime}_{m})}}~{}\frac{1}{|\omega_{m}-\omega^{\prime}_{m}|^{\gamma}},$
	$\displaystyle{\tilde{\Sigma}}(\omega_{m})$	$\displaystyle=$	$\displaystyle\omega_{m}+\frac{{\bar{g}}^{\gamma}}{2}\int d\omega^{\prime}_{m}\frac{{\tilde{\Sigma}}(\omega^{\prime}_{m})}{\sqrt{{\tilde{\Sigma}}^{2}(\omega^{\prime}_{m})+\Phi^{2}(\omega^{\prime}_{m})}}~{}\frac{1}{|\omega_{m}-\omega^{\prime}_{m}|^{\gamma}}$		(1)

where ${\tilde{\Sigma}}(\omega_{m})=\omega_{m}+\Sigma(\omega_{m})$. In these equations, both $\Sigma(\omega_{m})$ and $\Phi(\omega_{m})$ are real functions. Observe that we define $\Sigma(\omega_{m})$ with the overall plus sign and without the overall factor of $i$. ¹¹1In obtaining the Eliashberg equations we assumed that states away from the Fermi surface are irrelevant for the pairing and NFL behavior, and extended the integration over momenta to infinite limits. In this situation, fermionic self-energy $\Sigma(\omega_{m})$ is real. If the momentum integration over $k-k_{F}$ holds in finite and non-symmetric limits, $\Sigma(\omega_{m})$ has both real and imaginary components. The same holds in SYK-type models at a non-zero chemical potential (see Ref. Gu_2020 and references therein). In the normal state ($\Phi\equiv 0$),

$${\tilde{\Sigma}}(\omega_{m})=\omega_{m}+\omega_{0}^{\gamma}|\omega_{m}|^{1-\gamma}\operatorname{sgn}{\omega_{m}},$$

(2)

where $\omega_{0}={\bar{g}}/(1-\gamma)^{1/\gamma}$. At small $\gamma$, $\omega_{0}={\bar{g}}e$.

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

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

(3)

The equation for $\Delta(\omega_{m})$ is readily obtained from (II):

$$\Delta(\omega_{m})=\frac{{\bar{g}}^{\gamma}}{2}\int d\omega^{\prime}_{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}}.$$

(4)

This equation contains a single function $\Delta(\omega_{m})$, but for the cost that $\Delta(\omega_{m})$ appears also in the r.h.s., which makes Eq. (4) less convenient for the analysis than Eqs. (II).

Eqs. (II)-(4) describe color superconductivity Son (1999); *son2 and pairing in 3D ($\gamma=0_{+}$, $V(\Omega_{m})\propto\log{|\omega_{m}|}$), spin- and charge-mediated pairing in $D=3-\epsilon$ dimensions Mross et al. (2010); Metlitski et al. (2015); Raghu et al. (2015) and superconductivity in graphene Khveshchenko (2009) ($\gamma=O(\epsilon)\ll 1$), a 2D pairing  Altshuler et al. (1995b) with interaction peaked at $2k_{F}$ ($\gamma=1/4$), pairing at a 2D nematic/Ising-ferromagnetic QCP Bonesteel et al. (1996); Lederer et al. (2015); Wang et al. (2001); *triplet2; *triplet3 ($\gamma=1/3$), pairing at a 2D $(\pi,\pi)$ SDW QCP Millis (1992); Abanov et al. (2001a, 2003); Wang and Chubukov (2013a) and an incommensurate CDW QCP Castellani et al. (1995); *ital2; *ital3; Chowdhury and Sachdev (2014b); *wang_22; *wang23 ($\gamma=1/2$), dispersionless fermions randomly interacting with an Einstein phonon Esterlis and Schmalian (2019); Wang (2020); Hauck et al. (2019) and a spin-liquid model for the cuprates Tsvelik (2017) ($\gamma=0.7$) a 2D pairing mediated by an undamped propagating boson ($\gamma=1$), pairing in several Fe-based superconductors Haule and Kotliar (2007) ($\gamma=1.2$) and even the strong coupling limit of phonon-mediated superconductivity for either dispersion-full Combescot (1995); Bergmann and Rainer (1973); *Bergmann2; *ad; Marsiglio et al. (1988); *Marsiglio_91; Karakozov et al. (1991) or dispersion-less Esterlis and Schmalian (2019); Wang (2020) fermions ($\gamma=2$). The pairing models with parameter-dependent $\gamma$ have been analyzed as well (Refs. Sachdev et al., 2009; Moon and Chubukov, 2010). The case $\gamma=0$ describes a BCS superconductor. We list some of the model in Tables 1 - 3.

A justification of Eliashberg theory for electronically mediated superconductivity (e.g., the reasoning to neglect vertex corrections) is case specific. For Ising-nematic fluctuations (the case $\gamma=1/3$), vertex corrections are small in ${\bar{g}}/E_{F}$, which is a small parameter of the theory Metzner et al. (2003); *DellAnna2006; Maslov and Chubukov (2010); Klein et al. (2020); *avi_last_1. In other cases, e,g., in SYK models, a small parameter for Eliashberg approximation is $1/N$, where $N$ is the number of of fermionic flavors Patel and Sachdev (2019); Esterlis and Schmalian (2019); Wang (2020); Hauck et al. (2019); Chowdhury and Berg (2020). For several 2D models, the corrections to Eliashberg approximation for the self-energy in the normal state are logarithmically singular and in the absence of the pairing would change the system behavior at the smallest frequencies Abanov et al. (2003); Metlitski and Sachdev (2010a); Dalidovich and Lee (2013); *sslee_2018. One-loop logarithmic corrections just change $\gamma$ but keep the model intact (Ref. Abanov et al. (2003)), but higher-order corrections go beyond the $\gamma$ model (Refs. Metlitski and Sachdev (2010a); Dalidovich and Lee (2013); *sslee_2018). Here, we assume that the onset temperature for the pairing, $T_{p}$, is larger, at least numerically, than the scale at which corrections to Eliashberg approximation become relevant, and stick with the Eliashberg theory. Recent Quantum Monte Carlo (QMC) calculations for superconductivity, mediated by antiferromagnetic spin fluctuations Wang et al. (2017b) and Ising-nematic fluctuations Klein et al. (2020), have found the onset of superconductivity at a temperature, almost identical to the one in the $\gamma$ model with $\gamma=1/2$ and $1/3$, respectively, without vertex corrections. Normal state QMC calculations for the models with antiferromagnetic and Ising-nematic fluctuations also found Klein et al. (2020); *avi_last_1 very good agreement with the behavior of the $\gamma$ models with $\gamma=1/2$ and $1/3$. A good agreement has been also found between QMC calculations of superconducting $T_{c}$ for electron-phonon interaction and the $\gamma$-model with $\gamma=2$ (Ref. Chubukov et al. (2020a)), for couplings smaller than the Fermi energy.

In this paper, we consider the set of $\gamma$-models with $0<\gamma<1$. The analysis for $\gamma\geq 1$ requires separate consideration because of divergencies in the r.h.s. of both equations in (II) (but not in (4)) and will be presented in subsequent papers.

Notice that for any $\gamma>0$ the pairing interaction $V(\Omega)$ decays at large frequencies as $1/|\Omega|^{\gamma}$. A simple experimentation shows that the integrals in (II) and (4) are then convergent in the ultraviolet. The only exceptions are the BCS case $\gamma=0$, when $V(\Omega)$ is just a constant, and the case $V(\Omega_{m})\propto\log{|\omega_{m}|}$, i.e., $\gamma=0_{+}$ (color superconductivity/pairing in 3D, Refs. Son (1999); *son2; Metlitski et al. (2015)). For these two cases one needs to set the upper cutoff of frequency integration at some $\Lambda$ to cut ultra-violet divergence. We discuss the limit $\gamma\to 0$ in Appendix A.

The full set of Eliashberg equations for electron-mediated pairing contains also the equation describing the feedback from the pairing on the bosonic propagator. This feedback is small by $\delta_{E}$ in the case of electron-phonon interaction, but is generally not small when the pairing is mediated by a collective mode because the dispersion of a collective mode may change qualitatively below $T_{p}$. The most known example of this kind is the transformation of Landau overdamped spin collective mode in the normal state to a propagating mode (often called a resonance mode) below the onset of $d-$wave pairing mediated by antiferromagnetic spin fluctuations (see, e.g., Refs.Millis and Monien (1996); Eschrig (2006); Abanov et al. (2002); Abanov and Chubukov (2000)). To avoid additional complications, we do not include this feedback explicitly into our consideration. In general, the feedback from the pairing makes bosons less incoherent and can be modeled by assuming that $\gamma$ moves towards a larger value as $T$ decreases ( e.g., from $\gamma=1/2$ to $\gamma=1$ for the case of antiferromagnetic fluctuations and from $\gamma=1/3$ to $\gamma=2/3$ for the case of Ising-nematic fluctuations).

The coupled equations for $\Phi$ and ${\tilde{\Sigma}}$, Eq. (II), describe the interplay between the two competing tendencies – one towards superconductivity, specified by $\Phi$, and the other towards incoherent non-Fermi liquid behavior, specified by ${\tilde{\Sigma}}$. The competition between the two tendencies is encoded in the fact that ${\tilde{\Sigma}}$ appears in the denominator of the equation for $\Phi$ and $\Phi$ appears in the denominator of the equation for ${\tilde{\Sigma}}$. In more physical terms a self-energy ${\tilde{\Sigma}}$ is an obstacle to Cooper pairing, while if $\Phi$ is non-zero, it reduces the strength of the self-energy and moves the system back into a FL regime.

In Eq. (II) the couplings in the particle-particle and particle-hole channels have the same magnitude ${\bar{g}}$. To study the interplay, it is convenient to have a parameter, which would increase either the tendency towards NFL or towards pairing. With this in mind, we multiply the coupling in the particle-particle channel by a factor $1/N$, i.e., set it to be ${\bar{g}}^{\gamma}/N$ instead of $\bar{g}^{\gamma}$, and keep the coupling in the particle-hole channel intact. For $N<1$, the tendency towards pairing is enhanced, for $N>1$ the tendency towards NFL effectively gets larger. We will treat $N$ as a free parameter, but use the extension as just a way to see the relevant physics more clearly, as our ultimate goal is to understand system behavior in the physical case of $N=1$. We note in passing that the extension to integer $N>1$ can be formalized by extending the original model to matrix $SU(N)$ model Raghu et al. (2015).

The modified equations for $\Phi(\omega_{m})$ and ${\tilde{\Sigma}}(\omega_{m})$ are

	$\displaystyle\Phi(\omega_{m})=\frac{{\bar{g}}^{\gamma}}{2N}\int d\omega^{\prime}_{m}\frac{\Phi(\omega^{\prime}_{m})}{\sqrt{{\tilde{\Sigma}}^{2}(\omega^{\prime}_{m})+\Phi^{2}(\omega^{\prime}_{m})}}~{}\frac{1}{|\omega_{m}-\omega^{\prime}_{m}|^{\gamma}},$		(5)
	$\displaystyle{\tilde{\Sigma}}(\omega_{m})=\omega_{m}+\frac{{\bar{g}}^{\gamma}}{2}\int d\omega^{\prime}_{m}\frac{{\tilde{\Sigma}}(\omega^{\prime}_{m})}{\sqrt{{\tilde{\Sigma}}^{2}(\omega^{\prime}_{m})+\Phi^{2}(\omega^{\prime}_{m})}}~{}\frac{1}{|\omega_{m}-\omega^{\prime}_{m}|^{\gamma}},$		(6)

and the equation for $\Delta(\omega_{m})$ becomes

$$\Delta(\omega_{m})=\frac{{\bar{g}}^{\gamma}}{2N}\int d\omega^{\prime}_{m}\frac{\Delta(\omega^{\prime}_{m})-N\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}}.$$

(7)

Below we will occasionally refer to the equation for $\Phi(\omega_{m})$ as the gap equation, notwithstanding that the gap equation is given by Eqs. (7). Indeed, once we know $\Phi(\omega_{m})$ and ${\tilde{\Sigma}}(\omega_{m})$, we also know $\Delta(\omega_{m})=\Phi(\omega_{m})\omega_{m}/{\tilde{\Sigma}}(\omega_{m})$.

We will be searching for normalized solutions for $\Phi(\omega_{m})$ and $\Delta(\omega_{m})$. Physically, the normalizability of a solution follows from the requirement that the Free energy of a superconductor should be free from divergencies. The Free energy of a superconductor, $F_{sc}$, can be obtained by either applying Hubbard-Stratonovich formalism  Yuzbashyan et al. (a) or by using a generic Luttinger-Ward-Eliashberg expression Luttinger and Ward (1960); Eliashberg (1960), and explicitly integrating over momentum, approximating the density of states by its value at the Fermi level Haslinger and Chubukov (2003); Wu et al. (2019b). The result is

$$\frac{F_{sc}}{N_{0}}=-\int d\omega_{m}\frac{\omega^{2}_{m}}{\sqrt{\omega^{2}_{m}+\Delta^{2}(\omega_{m})}}-\frac{{\bar{g}}^{\gamma}}{4}\int d\omega_{m}d\omega^{\prime}_{m}\frac{\omega_{m}\omega^{\prime}_{m}+\frac{1}{N}\Delta(\omega_{m})\Delta(\omega^{\prime}_{m})}{\sqrt{\omega^{2}_{m}+\Delta^{2}(\omega_{m})}\sqrt{(\omega^{\prime}_{m})^{2}+\Delta^{2}(\omega^{\prime}_{m})}}\frac{1}{|\omega_{m}-\omega^{\prime}_{m}|^{\gamma}}.$$

(8)

The gap equation (7) is obtained from the condition $\delta F_{sc}/\delta\Delta(\omega_{m})=0$.

The condensation energy $E_{c}$ is the difference between $F_{sc}$, with $\Delta(\omega_{m})$ satisfying the Eliashberg equation (7), and $F_{n}$ (the free energy for $\Delta=0$). Using Eq. (8) and the gap equation (7) we obtain Wada (1964); Bardeen and Stephen (1964); Yuzbashyan et al. (a)

	$\displaystyle\frac{E_{c}}{N_{0}}=-\int d\omega_{m}\frac{|\omega_{m}|}{\sqrt{1+D^{2}(\omega_{m})}}\left(\sqrt{1+D^{2}(\omega_{m})}-1\right)^{2}-\frac{{\bar{g}}^{\gamma}}{4}\int\frac{d\omega_{m}d\omega^{\prime}_{m}}{\sqrt{1+D^{2}(\omega_{m})}\sqrt{1+D^{2}(\omega^{\prime}_{m})}}$
	$\displaystyle\times\left(\sqrt{1+D^{2}(\omega_{m})}-\sqrt{1+D^{2}(\omega^{\prime}_{m})}\right)^{2}\left(\frac{1}{|\omega_{m}-\omega^{\prime}_{m}|^{\gamma}}-\frac{1}{|\omega_{m}+\omega^{\prime}_{m}|^{\gamma}}\right)$		(9)

where $D(\omega_{m})=\Delta(\omega_{m})/\omega_{m}$. Both terms in the r.h.s. of (9) are negative, hence, the condensation energy is negative for any solution of the gap equation $\Delta(\omega_{m})$. Note that the factor $N$ is not present in (9).

To understand whether or not the ground state at $T=0$ is a NFL state or a superconducting state, we first consider infinitesimally small $\Phi(\omega_{m})$. The corresponding equation is obtained by neglecting $\Phi$ in the denominator of Eq. (5) and using (2) for ${\tilde{\Sigma}}$:

$$N\Phi(\omega_{m})=\frac{1-\gamma}{2}\int d\omega^{\prime}_{m}\frac{\Phi(\omega^{\prime}_{m})}{|\omega^{\prime}_{m}|^{1-\gamma}|\omega_{m}-\omega^{\prime}_{m}|^{\gamma}}\frac{1}{1+\left(\frac{|\omega^{\prime}_{m}|}{\omega_{0}}\right)^{\gamma}}$$

(10)

This is an equation for an eigenfunction of a linear operator, in which $N$ plays a role of the eigenvalue. Observe that the fermion-boson coupling ${\bar{g}}$ appears only in the last term in the denominator, via $\omega_{0}\propto{\bar{g}}$. Without this term, the r.h.s. of (10) is marginal by power counting (the total exponent in the denominator is $1-\gamma+\gamma=1$). Then, once we rescale frequency to ${\bar{\omega}}_{m}=\omega_{m}/\omega_{0}$, the equation for $\Phi({\bar{\omega}}_{m})$ becomes fully universal:

$$\Phi({\bar{\omega}}_{m})=\frac{1-\gamma}{2N}\int d{\bar{\omega}}^{\prime}_{m}\frac{\Phi({\bar{\omega}}^{\prime}_{m})}{|{\bar{\omega}}^{\prime}_{m}|^{1-\gamma}|{\bar{\omega}}_{m}-{\bar{\omega}}^{\prime}_{m}|^{\gamma}}\frac{1}{1+|{\bar{\omega}}^{\prime}_{m}|^{\gamma}}$$

(11)

The same is true for the linearized equation for $\Delta(\omega_{m})$:

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

(12)

Using ${\bar{g}}^{\gamma}=\omega^{\gamma}_{0}(1-\gamma)$ and introducing again ${\bar{\omega}}_{m}=\omega_{m}/\omega_{0}$, we re-express (12) as

$$\Delta({\bar{\omega}}_{m})=\frac{1-\gamma}{2N}\int d{\bar{\omega}}^{\prime}_{m}\frac{\Delta({\bar{\omega}}^{\prime}_{m})-N\Delta({\bar{\omega}}_{m})\frac{{\bar{\omega}}^{\prime}_{m}}{{\bar{\omega}}_{m}}}{|{\bar{\omega}}^{\prime}_{m}|}~{}\frac{1}{|{\bar{\omega}}_{m}-{\bar{\omega}}^{\prime}_{m}|^{\gamma}}.$$

(13)

For infinitesimal $\Phi$ and $\Delta$, $\Delta(\omega_{m})=\Phi(\omega_{m})/\left(1+(\omega_{0}/|\omega_{m}|)^{\gamma}\right)$, or, equivalently, $\Delta({\bar{\omega}}_{m})=\Phi({\bar{\omega}}_{m})/(1+|{\bar{\omega}}_{m}|^{-\gamma})$.

A conventional wisdom, borrowed from BCS theory, would imply that Eqs. (11) and (13) must have solutions at a single critical $N_{cr}$, separating NFL and superconducting states, assuming that such critical $N_{cr}$ exists. By the same logic, the full non-linear gap equation has no solution at $N>N_{cr}$, while at $N<N_{cr}$, the solution $\Phi(\omega_{m})$ (or $\Delta(\omega_{m})$) has a finite magnitude, which increases with $N_{cr}-N$.

We show below that the situation in the $\gamma$-model is different. Specifically, we show that critical $N_{cr}$ exists, but the linearized gap equation has a solution not only at $N=N_{cr}$, but for all $N<N_{cr}$. Furthermore, $N=N_{cr}$ turns out to be a multi-critical point of the $\gamma$-model (for $0<\gamma<1$), below which there exists an infinite number of solutions of the full non-linear gap equation, $\Delta_{n}(\omega_{m})$. The magnitude of the $n$th solution initially increases as $e^{-A_{n}/\sqrt{N_{cr}-N}}$, where the factor $A_{n}$ depends on the number $n$ of the solution.

In the next four sections we discuss the linearized gap equation, Eqs. (11) and (13). We return to the non-linear gap equation in Sec. VI.

To verify that a potential solution of the linearized gap equation is normalizable, we will need to analyze the Free energy (8) order $\Delta^{2}$. Expanding in (8) we obtain

	$\displaystyle F_{sc}$	$\displaystyle=$	$\displaystyle F_{norm}+\frac{N_{0}}{2N}J(\Delta,N)$		(14)
	$\displaystyle J(\Delta,N)$	$\displaystyle=$	$\displaystyle N\int d{\bar{\omega}}_{m}\frac{\Delta^{2}(\bar{\omega}_{m})(1+|\bar{\omega}_{m}|^{\gamma})}{|\bar{\omega}_{m}|^{1+\gamma}}-\frac{1-\gamma}{2}\int d\omega_{m}d\omega^{\prime}_{m}\frac{\Delta(\omega_{m})\Delta(\omega^{\prime}_{m})}{|\bar{\omega}_{m}||\bar{\omega}^{\prime}_{m}|}\frac{1}{|\bar{\omega}_{m}-\bar{\omega}^{\prime}_{m}|^{\gamma}}.$

The linearized gap equation (13) is obtained by varying this $F_{sc}$ over $\Delta(\omega_{m})$. A solution is normalizable if the variation $\delta J(\Delta,N)/\delta N$ is finite, i.e., if

$$\int d\bar{\omega}_{m}\frac{\Delta^{2}(\bar{\omega}_{m})(1+|\bar{\omega}_{m}|^{\gamma})}{|\bar{\omega}_{m}|^{1+\gamma}}$$

(15)

is non-divergent. In terms of $\Phi(\bar{\omega}_{m})$, the same integral is

$$\int d\bar{\omega}_{m}\frac{\Phi^{2}(\bar{\omega}_{m})}{|\bar{\omega}_{m}|^{1-\gamma}(1+|\bar{\omega}_{m}|^{\gamma})}$$

(16)