Interplay between superconductivity and non-Fermi liquid at a quantum-critical point in a metal.
III: The $\gamma$ model and its phase diagram across $\gamma=1$ .

Yi-Ming Wu 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 our analysis of the interplay between the pairing and the non-Fermi liquid behavior in a metal for a set of quantum-critical models with an effective dynamical electron-electron interaction $V(\Omega_{m})\propto 1/|\Omega_{m}|^{\gamma}$ (the $\gamma$ -model). We analyze both the original model and its extension, in which we introduce an extra parameter $N$ to account for non-equal interactions in the particle-hole and particle-particle channel. In two previous papers,Abanov and Chubukov (2020); Wu et al. (2020) we considered the case $0<\gamma<1$ and argued that (i) at $T=0$ , there exists an infinite discrete set of topologically different gap functions, $\Delta_{n}(\omega_{m})$ , all with the same spatial symmetry, and (ii) each $\Delta_{n}$ evolves with temperature and terminates at a particular $T_{p,n}$ . In this paper, we analyze how the system behavior changes between $\gamma<1$ and $\gamma>1$ , both at $T=0$ and a finite $T$ . The limit $\gamma\to 1$ is singular due to infra-red divergence of $\int d\omega_{m}V(\Omega_{m})$ , and the system behavior is highly sensitive to how this limit is taken. We show that for $N=1$ , the divergencies in the gap equation cancel out, and $\Delta_{n}(\omega_{m})$ gradually evolve through $\gamma=1$ both at $T=0$ and a finite $T$ . For $N\neq 1$ , divergent terms do not cancel, and a qualitatively new behavior emerges for $\gamma>1$ . Namely, the form of $\Delta_{n}(\omega_{m})$ changes qualitatively, and the spectrum of condensation energies, $E_{c,n}$ becomes continuous at $T=0$ . We introduce different extension of the model, which is free from singularities for $\gamma>1$ .

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) in a metal with four-fermion interaction, mediated by a critical soft boson. We consider a class of models, for which soft bosons are slow modes compared to dressed electrons. In this situation, the low-energy physics at a QCP is governed by an effective dynamical interaction $V(\Omega_{m})={\bar{g}}^{\gamma}/|\Omega_{m}|^{\gamma}$ integrated along the Fermi surface (the $\gamma$ -model). This interaction is singular, and gives rise to two opposite tendencies: NFL behavior in the normal state, with fermionic self-energy $\Sigma(\omega_{m})\propto\omega^{1-\gamma}_{m}$ , and an attraction in at least one pairing channel. The two tendencies compete with each other as a NFL self-energy reduces the magnitude of the pairing kernel, while the feedback from the pairing reduces fermionic self-energy.

In the first paper in the series, 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 an extensive list of references to earlier publications on this subject. In this and in the subsequent paper, Ref. Wu et al. (2020), hereafter referred to as Paper I and Paper II, respectively, we analyzed the behavior of the $\gamma$ -model for $0<\gamma<1$ at $T=0$ (Paper I) and at a finite $T$ (Paper II). We found that the system does become unstable towards pairing. However, in qualitative distinction with BCS/Eliashberg theory of superconductivity, in which there is a single solution of the gap equation, $\Delta(\omega_{m})$ , here we found an infinite discrete set of solutions $\Delta_{n}(\omega_{m})$ . All solutions have the same spatial symmetry, but are topologically distinct as $\Delta_{n}(\omega_{m})$ changes sign $n$ times as a function of Matsubara frequency (each such point is a center of a dynamical vortex). The gap functions $\Delta_{n}(\omega_{m})$ with finite $n$ tend to finite values at zero frequency, but the magnitude of $\Delta_{n}(0)$ decreases with $n$ and at large enough $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 the solution of the linearized gap equation. We found the exact form of $\Delta_{\infty}(\omega_{m})$ . It oscillates as a function of $\log(|\omega_{m}|/{\bar{g}})$ down to the lowest frequencies and up to $\omega_{max}$ , which is generally of order ${\bar{g}}$ , except for the smallest $\gamma$ , where $\omega_{max}\sim{\bar{g}}(1/\gamma)^{1/\gamma}$ . At $\omega>\omega_{max}$ , $\Delta_{\infty}(\omega_{m})$ decays as $1/|\omega_{m}|^{\gamma}$ . A function $\Delta_{n}(\omega_{m})$ with a finite $n$ saturates at $\Delta_{n}(0)$ below $\omega_{m}\sim\Delta_{n}(0)$ , and at larger $\omega$ retains the functional form of $\Delta_{\infty}(\omega_{m})$ . At a finite $T$ , each $\Delta_{n}(\omega_{m})$ evolves with $T$ and terminates at its own $T_{p,n}\sim\Delta_{n}(0)$ .

In this paper we extend the analysis to larger $\gamma$ . We will be particularly interested in the evolution of the system behavior between $\gamma\leq 1$ and $\gamma\geq 1$ . At $T=0$ , the frequency integral $\int d\Omega_{m}V(\Omega_{m})\propto\int d\Omega_{m}/|\Omega_{m}|^{\gamma}$ diverges at small $\Omega_{m}$ for $\gamma\geq 1$ , and from a general perspective one could expect that this divergence introduces qualitative changes in the system behavior. Indeed, the pairing vertex and the fermionic self-energy at $T=0$ do become singular for $\gamma\geq 1$ . We show, however, that singular terms cancel out in the equation for the gap function $\Delta(\omega_{m})$ . As a consequence, for both $\gamma\leq 1$ and $\gamma\geq 1$ the full non-linear gap equation has an infinite number of solutions $\Delta_{n}(\omega_{m})$ at $T=0$ and each $\Delta_{n}$ terminates at its own $T_{p,n}$ . All functions $\Delta_{n}(\omega_{m})$ evolve smoothly through $\gamma=1$ . The corresponding condensation energies, $E_{c,n}$ form a discrete set in which $E_{c,0}$ is the largest.

We next analyze a more general model with different interaction strength in particle-particle and particle-hole channels. A natural way to account for this is to multiply the interaction in the pairing channel by a factor $1/N$ leaving the interaction in the particle-hole channel intact Abanov and Chubukov (2020); Wu et al. (2020); Wang et al. (2016); Chubukov et al. (2020a); Abanov et al. (2019); *Wu_19_1 Another way to split the strength of the two interactions is to extend the original $\gamma-$ model to matrix $SU(N)$ modelRaghu et al. (2015); Wang et al. (2017, 2018).

The factor $N$ plays the role of the eigenvalue in the linearized matrix gap equation, and understanding the system behavior for $N\neq 1$ is also essential for the interpretation of the flow of the eigenvalues and the eigenfunctions in the numerical analysis of the gap equation even for $N\to 1$ . We show that for $\gamma\to 1$ , the system behavior becomes very sensitive to small deviations from $N=1$ , because the $T=0$ divergencies in the self-energy and the pairing vertex do not cancel for any $N\neq 1$ . As a consequence, the limits $\gamma\to 1$ and $N\to 1$ do not commute, and the structure of the gap function strongly depends on the ratio $(N-1)/(1-\gamma)$ . We show that the behavior at $\gamma\to 1$ and $N<1$ is qualitatively different from that for $N\equiv 1$ . Namely, for $N<1$ the set of condensation energies becomes a continuous one at $\gamma=1$ : $E_{c,n}$ for all finite $n$ become the same as $E_{c,0}$ and $E_{c,\infty}$ form a continuous one-parameter gapless spectrum, similar to how a continuous phonon spectrum emerges in a continuum limit. This opens up a channel of massless ’longitudinal” gap fluctuations, in a truly qualitative distinction from BCS-type physics. We show that this behavior holds at $T=0$ for $\gamma>1$ , and that the structure of $\Delta_{n}(\omega_{m})$ at a finite $T$ also becomes qualitatively different from that for $\gamma<1$ and gives rise to highly unconventional form of the density of eigenvalues for $N<1$ , as we show both analytically and numerically.

We also discuss another extension of the theory, which does not introduce singular contributions that were responsible for qualitatively different behavior at $N=1$ and $N\neq 1$ . This extension also allows one to vary the relative strength of the interactions in the particle-hole and particle-particle channels (although in a less obvious way), and tune between NFL and SC states for $\gamma>1$ , similar to how it was done before for $\gamma<1$ in Refs. Wang et al. (2016); Chubukov et al. (2020a); Abanov et al. (2019); *Wu_19_1; Abanov and Chubukov (2020); Wu et al. (2020); Raghu et al. (2015); Wang et al. (2017, 2018).

The structure of the paper is the following. In Sec. II we briefly review the $\gamma$ model and present the equations for the pairing vertex, the self-energy, and the gap function, which we will use later in the paper. In Sec. III we show that the singularities, imposed by the divergence of $\int d\Omega_{m}V(\omega_{m})$ for $\gamma\geq 1$ , cancel out in the gap equation for $N=1$ . We argue that the full non-linear gap equation at $T=0$ has an infinite set of solutions $\Delta_{n}(\omega_{m})$ both for $\gamma\leq 1$ and $\gamma\geq 1$ , and show that the solutions vary smoothly through $\gamma=1$ . In Sec. IV we extend the $\gamma$ model to $N\neq 1$ and show that for a generic $N$ the system behavior changes qualitatively between $\gamma<1$ and $\gamma\geq 1$ . We discuss the double limit $\gamma\to 1$ , $N\to 1$ at $T=0$ , show how the set of the condensation energies, $E_{c,n}$ , becomes continuous at $\gamma>1$ , and discuss the new structure of $\Delta_{n}(\omega_{m})$ at $\gamma>1$ and a finite $T$ . In Sec. V we discuss another extension of the $\gamma$ model, which does not introduce the divergencies.

In Paper IV, the next in the series, we consider the case $N=1$ , $1<\gamma<2$ in more detail, and argue that as $\gamma$ increases, the dynamical vortices emerge one by one and form an array in the upper frequency half-plane. The number of vortices tends to infinity for $\gamma\to 2$ .

II $\gamma$ -model, Eliashberg equations

The $\gamma$ -model was introduced in Paper I and in earlier publications as a low-energy model for the interaction between soft bosons and electrons Abanov and Chubukov (2020); Wu et al. (2020); Abanov et al. (2001, 2003); Moon and Chubukov (2010); Metlitski and Sachdev (2010); Mross et al. (2010); Monthoux et al. (2007); Efetov et al. (2013); Metlitski et al. (2015); Raghu et al. (2015); Wang et al. (2016); Lee et al. (2018); Abanov et al. (2019); Wu et al. (2019a); Chubukov et al. (2020a), and we refer the reader to these works for the justification of the model and its relation to various quantum-critical systems. The model describes low-energy fermions with an effective dynamical interaction $V(\Omega_{m})={\bar{g}}^{\gamma}/|\Omega_{m}|^{\gamma}$ , averaged over momenta on the Fermi surface with a proper weight. The case $\gamma\approx 1$ corresponds to, e.g., pairing by a weakly damped soft optical phonon with static susceptibility peaked some finite momentum $Q_{0}$ .Chubukov et al. (2020b) The coupled equations for the fermionic self-energy $\Sigma(\omega_{m})$ and the pairing vertex $\Phi(\omega_{m})$ in the most attractive pairing channel are similar to Eliashberg equations for the case of a dispersionless phonon, and we will use the term “Eliashberg equations” for our case.

At a finite $T$ the coupled Eliashberg equations for $\Phi(\omega_{m})$ and $\Sigma(\omega_{m})$ are, in Matsubara formalism

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

where ${\tilde{\Sigma}}(\omega_{m})=\omega_{m}+\Sigma(\omega_{m})$ . In these notations, $\Sigma(\omega_{m})$ is a real function, odd in frequency.

The SC 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}}

(2)

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

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

(3)

This equation contains a single function $\Delta(\omega_{m})$ , but at the cost that $\Delta(\omega_{m})$ appears also in the r.h.s. Both $\Phi(\omega_{m})$ and $\Delta(\omega_{m})$ are defined up to an overall $U(1)$ phase factor, which we set to zero for definiteness. Eqs. (II) and (3) exclude the self-action term with $m^{\prime}=m$ . This term cancels out by Anderson theorem Abrikosov et al. (1965), because scattering with zero frequency transfer mimics the effect of scattering by non-magnetic impurities.

Below we will analyze the full non-linear equations and the linearized equations, for infinitesimally small $\Phi(\omega_{m})$ and $\Delta(\omega_{m})$ . The latter determine, e.g., critical temperatures $T_{p,n}$ . The linearized gap equation is

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

(4)

The linearized equation for the pairing vertex $\Phi(\omega_{m})$ is

\Phi(\omega_{m})={\bar{g}}^{\gamma}\pi T\sum_{m^{\prime}\neq m}\frac{\Phi(\omega_{m^{\prime}})}{|\omega_{m^{\prime}}+\Sigma_{\text{norm}}(\omega_{m^{\prime}})|}~{}\frac{1}{|\omega_{m}-\omega_{m^{\prime}}|^{\gamma}}

(5)

where $\Sigma_{\text{norm}}(\omega_{m})$ is the self-energy of the normal state,

\Sigma_{\text{norm}}(\omega_{m})={\bar{g}}^{\gamma}(2\pi T)^{1-\gamma}\sum_{m^{\prime}=1}^{m}\frac{1}{|m^{\prime}|^{\gamma}}={\bar{g}}^{\gamma}(2\pi T)^{1-\gamma}H_{m,\gamma}

(6)

and $H_{m,\gamma}$ is the Harmonic number. This expression holds for $\omega_{m}\neq\pm\pi T$ . For the two lowest Matsubara frequencies, $\Sigma_{\text{norm}}(\pm\pi T)=0$ . We emphasize that $\Sigma_{\text{norm}}(\omega_{m})$ in (6) is not the full normal state self-energy, as the summation in (5) excludes the term $m^{\prime}=m$ .

At $T=0$ ,

	$\displaystyle\Delta(\omega_{m})$	$\displaystyle=$	$\displaystyle\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}},$
	$\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}}$		(7)

The linearized equations are

	$\displaystyle\Delta(\omega_{m})$	$\displaystyle=$	$\displaystyle\frac{{\bar{g}}^{\gamma}}{2}\int d\omega^{\prime}_{m}\frac{\Delta(\omega^{\prime}_{m})-\Delta(\omega_{m})\frac{\omega^{\prime}_{m}}{\omega_{m}}}{\|\omega^{\prime}_{m}\|}~{}\frac{1}{\|\omega_{m}-\omega^{\prime}_{m}\|^{\gamma}},$
	$\displaystyle\Phi(\omega_{m})$	$\displaystyle=$	$\displaystyle\frac{{\bar{g}}^{\gamma}}{2}\int d\omega^{\prime}_{m}\frac{\Phi(\omega^{\prime}_{m})}{\|\omega^{\prime}_{m}+\Sigma_{\text{norm}}(\omega^{\prime}_{m})\|}~{}\frac{1}{\|\omega_{m}-\omega^{\prime}_{m}\|^{\gamma}}$		(8)

where

\Sigma_{\text{norm}}(\omega_{m})=\omega^{\gamma}_{0}|\omega_{m}|^{1-\gamma}{\text{sgn}}(\omega_{m})

(9)

and

\omega_{0}={\bar{g}}/(1-\gamma)^{1/\gamma}.

(10)

At small $\gamma$ , $\omega_{0}={\bar{g}}e$ . At $\gamma\to 1$ , $\omega^{\gamma}_{0}$ diverges as $1/(1-\gamma)$ .

III Transformation from $\gamma\leq 1$ to $\gamma\geq 1$

Refer to caption — Figure 1: The solution of $\epsilon_{i\beta}=1$ , with $\epsilon_{i\beta}$ given by Eq. (12), as a function of the exponent $\gamma$ .

III.1 A generic $\gamma<1$

We found in Papers I and II that

•

The non-linear gap equation has an infinite discrete set of solutions $\Delta_{n}(\omega_{m})$ , $n=0,1,2...$ . All $\Delta_{n}(\omega_{m})$ with finite $n$ tend to finite $\Delta_{n}(0)$ at zero frequency and decay as $1/|\omega_{m}|^{\gamma}$ at large frequencies. The function $\Delta_{n}(\omega_{m})$ changes sign $n$ times. At large $n$ , $\Delta_{n}(0)\propto e^{-An}$ , where $A=O(1)$ is a function of $\gamma$ .
•

The end point of the set, $\Delta_{\infty}(\omega_{m})$ , is the solution of the linearized gap equation. At small $\omega_{m}\ll\omega_{0}$ ,

$\Delta_{\infty}(\omega_{m})=C|\omega_{m}|^{\gamma/2}\cos\left({\beta\log{\left(\frac{|\omega_{m}|}{\omega_{0}}\right)^{\gamma}}+\phi}\right)$ (11)

where $\beta$ depends on $\gamma$ (see Eq. (12) below and Fig.1), and $\phi$ is some $\gamma$ -dependent number. Eq. (11) is readily obtained if one neglects $\omega_{m}$ compared to $\Sigma_{\text{norm}}(\omega_{m})\propto\omega^{1-\gamma}_{m}$ in Eq. (8) for the pairing vertex. The corrections to log-oscillating form hold in powers of $z=(|\omega_{m}|/\omega_{0})^{\gamma}=\omega_{m}/\Sigma_{\text{norm}}(\omega_{m})$ .
•

We found the exact form of $\Delta_{\infty}(z)$ for all $z$ . It oscillates up to $z=O(1)$ and decays as $1/z$ at larger $z$ . A gap function $\Delta_{n}(z)$ with a finite $n$ also decays as $1/z$ for $z>1$ , oscillates $n$ times at smaller $z$ , and saturates at the lowest frequencies at a finite $\Delta_{n}(0)$ .
•

At a finite $T$ , each $\Delta_{n}(\omega_{m})$ develops at the onset temperature $T_{p,n}\sim\Delta_{n}(0)$ . At large $n$ , $T_{p,n}\propto e^{-An}$ . The magnitude of $\Delta_{n}(\omega_{m})$ increases with decreasing $T$ , and at $T=0$ it coincides with the $n$ -th solution of the non-linear gap equation.

III.2 $\gamma\approx 1$

III.2.1 Linearized gap equation, $T=0$

We now analyze what happens when $\gamma$ increases and approaches $1$ . We begin with the linearized gap equation for $\Delta_{\infty}(\omega_{m})$ . At low frequencies, the solution is Eq. (11). The pre-logarithmic factor $\beta$ there is the root of $\epsilon_{i\beta}=1$ , where

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

(12)

The solution exists for all $\gamma<1$ , and $\beta$ approaches a finite value $0.792$ when $\gamma\to 1$ (Fig. 1). However, other quantities do become singular at $\gamma\to 1$ . We see from (9) and (10) that the normal state self-energy diverges because $\omega_{0}={\bar{g}}/(1-\gamma)^{1/\gamma}\to\infty$ . Accordingly, $z=\omega_{m}/\Sigma_{\text{norm}}(\omega_{m})$ remains small at frequencies of order ${\bar{g}}$ and becomes $O(1)$ only at a much larger $\omega_{m}\sim\omega_{0}$ . We show this in Fig.2.

Taken at a face value, this would imply that at $\gamma=1$ , the corrections from the expansion in $\omega_{m}/\Sigma_{\text{norm}}(\omega_{m})$ become totally irrelevant, and log-oscillations of $\Delta_{\infty}(z)$ extend to all frequencies. This would have a profound effect on the behavior of all other $\Delta_{n}(z)$ and on $T_{p,n}$ , as it is set by a frequency at which log-oscillations end.

We show that this is not the case, and $\Delta_{\infty}(\omega_{m})$ evolves smoothly through $\gamma=1$ . Namely $\Delta_{\infty}(\omega_{m})$ displays log-oscillations only up to $\omega_{m}=O({\bar{g}})$ , even at $\gamma\to 1$ , and decays as $1/z$ at larger frequencies. We show that this happens because the expansion in $\omega_{m}/\Sigma_{\text{norm}}(\omega_{m})$ in the limit $\gamma\rightarrow 1$ actually holds in

y=\frac{z}{1-\gamma}=\left(\frac{|\omega_{m}|}{\bar{g}}\right)^{\gamma}

(13)

so that the singularity in $z$ is canceled in this limit.

To demonstrate this, we analyze the structure of the corrections to the log-oscillating form of $\Delta_{\infty}(z)$ . As we discussed in Paper I, there are two types of corrections from the expansion in $\omega_{m}/\Sigma_{\text{norm}}(\omega_{m})$ : local corrections, which come from fermions with frequencies of order $\omega_{m}$ , and non-local corrections, which come from fermions with frequencies of order ${\bar{g}}$ : $\Delta_{\infty}(\omega_{m})=\Delta_{\infty,L}(\omega_{m})+\Delta_{\infty,NL}(\omega_{m})$ The expansion in powers of $\omega_{m}/\Sigma_{\text{norm}}(\omega_{m})$ comes from the local corrections, and we analyze now the structure of these corrections for $\gamma\to 1$ .

The series of local corrections can be obtained analytically in the order-by-order expansion. For a generic $\gamma<1$ , this expansion holds in powers of $z$ with prefactors of order one. Specifically,

\Delta_{\infty,L}(z)\propto\frac{\sqrt{z}}{1+z}{\text{R}e}\left[e^{(i\beta\log{z}+i\phi)}\sum_{m=0}^{\infty}C_{m}~{}z^{m}\right]

(14)

Here $C_{m}$ , subject to $C_{0}=1$ , are complex coefficients given by

C_{m>0}=I_{m}\displaystyle\prod_{{m^{\prime}}=1}^{m}\frac{1}{I_{m^{\prime}}-1},

(15)

where

	$\displaystyle I_{m^{\prime}}=\frac{(1-\gamma)}{2}\frac{\Gamma((m^{\prime}+1/2)\gamma+i{\beta}\gamma)\Gamma((1/2-m^{\prime})\gamma-i{\beta}\gamma)}{\Gamma(\gamma)}+$
	$\displaystyle\frac{\Gamma(2-\gamma)}{2}\left(\frac{\Gamma((m^{\prime}+1/2)\gamma+i{\beta}\gamma)}{\Gamma(1-(1/2-m^{\prime})\gamma+i{\beta}\gamma)}+\frac{\Gamma((1/2-m^{\prime})\gamma-i{\beta}\gamma)}{\Gamma(1-(m^{\prime}+1/2)\gamma-i{\beta}\gamma)}\right),$		(16)

and $\Gamma(...)$ are Gamma-functions. The phase $\phi$ is a free parameter in $\Delta_{\infty,L}(z)$ . Its value is set by the requirement that the total $\Delta_{\infty}(z)=\Delta_{\infty,L}(z)+\Delta_{\infty,NL}(z)$ decay as $1/z$ at large $z$ .

For $\gamma\approx 1$ , all $I_{m^{\prime}}$ tend to 1, and the coefficients $C_{m}$ become singular. Expanding $I_{m^{\prime}}$ in (16) near $\gamma=1$ , we obtain $I_{m^{\prime}}=1+(1-\gamma){\bar{I}}_{m^{\prime}}$ , where

	$\displaystyle{\bar{I}}_{m^{\prime}}=\frac{((-1)^{m^{\prime}}-1)\pi}{2\cosh{(\pi\beta)}}+$
	$\displaystyle\frac{1}{2}\left(\Psi{(1/2+i\beta)}+\Psi{(1/2-i\beta)}-\Psi{(1/2+m^{\prime}+i\beta)}-\Psi{(1/2-m^{\prime}-i\beta)}\right)$
	$\displaystyle=\frac{((-1)^{m^{\prime}}-1)\pi}{2\cosh{(\pi\beta)}}-\sum_{p=0}^{m^{\prime}-1}\frac{1}{1/2+i\beta+p}$		(17)

where $\Psi(...)$ is a di-Gamma function. Substituting into (15), we find that the coefficients $C_{m}$ scale as $C_{m}~{}\propto 1/(1-\gamma)^{m}$ . Substituting these $C_{m}$ into (14) we find that the expansion actually holds in $z/(1-\gamma)=(|\omega_{m}|/{\bar{g}})^{\gamma}$ , which is non-critical at $\gamma=1$ . The corrections to log-oscillating behavior then become relevant at a finite characteristic frequency $\omega_{m}\sim{\bar{g}}$ . The same behavior can be detected by plotting the exact solution for $\Delta_{\infty}(\omega_{m})$ for $\gamma\to 1$ . We present the plot in Fig. 3. We see that $\Delta_{\infty}(\omega_{m})$ indeed oscillates up to $\omega_{m}\sim{\bar{g}}$ and then decays as $1/|\omega_{m}|$ . We emphasize again that the largest scale for the oscillations is a finite ${\bar{g}}$ , despite that the expansion in frequencies in the exact solution formally holds in powers of $z\propto(1-\gamma)$ . We discuss this issue in Appendix A.

We also note that at $m\gg 1$ , $z^{m}C_{m}\propto(-1)^{m}((|\omega_{m}|/{\bar{g}})/\log m)^{m}$ . Because of $\log m$ in the denominator, the series in (14) converge absolutely, i.e., one can obtain $\Delta_{\infty,L}(\omega_{m})$ for any $|\omega_{m}|/{\bar{g}}$ by summing up enough terms in the perturbation series, although in practice it can be done only up to some $\omega/{\bar{g}}\geq 1$ . We plot the result of the summation of $1000$ terms in Fig.4 along with the exact $\Delta_{\infty}(\omega_{m})$ for $\gamma=0.9999$ . We see that over the whole frequency range where $\Delta_{\infty}(\omega_{m})$ oscillates, it practically coincides with $\Delta_{\infty,L}(\omega_{m})$ . To reproduce the $1/|\omega_{m}|$ behavior at larger frequencies we would need to include the non-local part $\Delta_{\infty,NL}(\omega_{m})$ .

III.2.2 Nonlinear gap equation, $T=0$

We now look at the evolution of $\Delta_{n}(\omega_{m})$ with some finite $n$ . At $T=0$ , $\Delta_{n}(\omega_{m})$ tends to a finite value at $\omega_{m}\to 0$ , and we first check whether $\Delta_{n}(0)$ remain continuous through $\gamma=1$ .

The gap function $\Delta_{n}(\omega_{m})$ is the solution of the non-linear gap equation (8). For $\gamma<1$ , one can safely move the term with $\Delta(\omega_{m})$ to the l.h.s. of the gap equation and re-express it as

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

Each integral is non-singular in the infra-red limit, provided that $\Delta_{n}(0)$ is finite. Then relevant $\omega^{\prime}_{m}$ are finite, and at small $\omega$ , one can expand in the integrands as

\frac{1}{|\omega^{\prime}_{m}\mp\omega_{m}|^{\gamma}}\approx\frac{1}{|\omega^{\prime}_{m}|^{\gamma}}\left(1\pm\gamma\frac{\omega_{m}}{\omega^{\prime}_{m}}\right)

(19)

Substituting the expansion into (18) and taking the limit $\omega_{m}\to 0$ , we obtain the condition on $\Delta_{n}(0)$ :

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

(20)

At $\gamma\to 1$ , each integral diverges as $1/(1-\gamma)$ , but the divergent terms cancel each other. As the result, $\Delta_{n}(0)$ remain finite at $\gamma\to 1$ . To see this more explicitly, consider the solution with $n=0$ . A sign-preserving $\Delta_{0}(\omega_{m})$ remains roughly equal to $\Delta_{0}(0)$ up to $\omega_{m}\sim\Delta_{0}(0)$ , at which both integrals in (20) already converge. Approximating then $\Delta_{0}(\omega^{\prime}_{m})$ by $\Delta_{0}(0)$ , we obtain from (20)

1={\bar{g}}^{\gamma}(1-\gamma)\int_{0}^{\infty}\frac{d\omega^{\prime}_{m}}{|\omega^{\prime}_{m}|^{\gamma}\sqrt{(\Delta_{0}(0))^{2}+(\omega^{\prime}_{m})^{2}}}=(1-\gamma)\left(\frac{\bar{g}}{\Delta_{0}(0)}\right)^{\gamma}\frac{\Gamma\left(\frac{1}{2}-\frac{\gamma}{2}\right)\Gamma\left(\frac{\gamma}{2}\right)}{2\sqrt{\pi}}.

(21)

This yields

\Delta_{0}(0)=\bar{g}\left[\frac{(1-\gamma)\Gamma\left(\frac{1}{2}-\frac{\gamma}{2}\right)\Gamma\left(\frac{\gamma}{2}\right)}{2\sqrt{\pi}}\right]^{1/\gamma}\approx\bar{g}\left(1+(1-\gamma)\log 2\right)

(22)

We see that $\Delta_{0}(0)={\bar{g}}$ at $\gamma\to 1$ from below. At smaller $\gamma$ , $\Delta_{0}(0)$ increases in the same way as $T_{c,0}$ , and the ratio $2\Delta_{0}(0)/T_{c,0}$ remains of order one. This is consistent with the more detailed study of $2\Delta_{0}(0)/T_{c,0}$ ratio in Ref.Wu et al. (2019b)

For $\gamma>1$ , $\int dx/|x|^{\gamma}$ diverges, and one cannot separate the two terms in the r.h.s. of (8). However, we can now use the identity

\int_{-\infty}^{\infty}d\omega^{\prime}_{m}\frac{1-\frac{\omega^{\prime}_{m}}{\omega_{m}}}{|\omega_{m}-\omega^{\prime}_{m}|^{\gamma}}=0

(23)

This identity holds for $\gamma>1$ , but not for $\gamma\leq 1$ . Using (23), we re-express the equation on $\Delta_{n}(\omega_{m})$ as

	$\displaystyle\Delta_{n}(\omega_{m})=\frac{{\bar{g}}^{\gamma}}{2}\int_{-\infty}^{\infty}\frac{d\omega^{\prime}_{m}~{}\left(\Delta_{n}(\omega^{\prime}_{m})-\Delta_{n}(\omega_{m})\right)}{\sqrt{\Delta^{2}_{n}(\omega^{\prime}_{m})+(\omega^{\prime}_{m})^{2}}\|\omega_{m}-\omega^{\prime}_{m}\|^{\gamma}}$
	$\displaystyle-\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[\sqrt{\Delta^{2}_{n}(\omega^{\prime}_{m})+(\omega^{\prime}_{m})^{2}}-\Delta_{n}(\omega^{\prime}_{m})\right]$		(24)

Each integral in (24) is now regular. In the limit $\omega_{m}\to 0$ , one can again use (19) and obtain

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

Like we did for $\gamma<1$ , we set $n=0$ and approximate $\Delta_{0}(\omega_{m})$ by $\Delta_{0}(0)$ . Substituting into (25), we obtain

	$\displaystyle\Delta_{0}(0)$	$\displaystyle=$	$\displaystyle{\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)}{\sqrt{\Delta^{2}_{0}(0)+(\omega^{\prime}_{m})^{2}}\|\omega^{\prime}_{m}\|^{\gamma}}$		(26)
			$\displaystyle=\Delta_{0}(0)\left(\frac{\bar{g}}{\Delta_{0}(0)}\right)^{\gamma}\left[(1-\gamma)\frac{\Gamma\left(\frac{1}{2}-\frac{\gamma}{2}\right)\Gamma\left(\frac{\gamma}{2}\right)}{2\sqrt{\pi}}\right]$		(26)

This gives exactly the same $\Delta_{0}(0)$ as (22). This proves that $\Delta_{0}(0)$ evolves continuously through $\gamma=1$ .

The verification that the same holds for $\Delta_{n}$ with a finite $n>0$ requires more efforts as one has to solve the actual non-linear gap equation for $\gamma<1$ and $\gamma>1$ and check whether the solutions match at $\gamma=1$ . This is technically quite challenging, but from physics perspective one should indeed expect $\Delta_{n}(\omega_{m})$ to vary continuously through $\gamma=1$ .

III.2.3 Linearized gap equation, finite $T$

We next analyze how the onset temperatures for the pairing, $T_{p,n}$ change around $\gamma=1$ . For a generic $\gamma<1$ , we found in paper II that at large $n$ , $T_{p,n}\propto e^{-\pi n/(\gamma\beta)}$ . We now show that this relation holds also for $\gamma\geq 1$ , but the derivation requires more efforts than for $\gamma<1$ .

The computations are more transparent when done for the pairing vertex $\Phi(\omega_{m})$ , expressed via the normal state ${\tilde{\Sigma}}_{\text{norm}}(\omega_{m})$ . The gap function $\Delta(\omega_{m})=\Phi(\omega_{m})\omega_{m}/{\tilde{\Sigma}}_{\text{norm}}(\omega_{m})$ . We have from (5)

	$\displaystyle\Phi(\omega_{m})$	$\displaystyle=$	$\displaystyle{\bar{g}}^{\gamma}\pi T\sum_{m^{\prime}\neq m}\frac{\Phi(\omega_{m^{\prime}})}{\|{\tilde{\Sigma}}_{\text{norm}}(\omega_{m^{\prime}})\|}~{}\frac{1}{\|\omega_{m}-\omega_{m^{\prime}}\|^{\gamma}},$
	$\displaystyle{\tilde{\Sigma}}_{\text{norm}}(\omega_{m})$	$\displaystyle=$	$\displaystyle\omega_{m}+{\bar{g}}^{\gamma}\pi T\sum_{m^{\prime}\neq m}~{}\frac{{\text{sign}}(\omega_{m^{\prime}})}{\|\omega_{m}-\omega_{m^{\prime}}\|^{\gamma}}$		(27)

Evaluating ${\tilde{\Sigma}}_{\text{norm}}(\omega_{m})\equiv{\tilde{\Sigma}}_{\text{norm}}(m)$ , we obtain

{\tilde{\Sigma}}_{\text{norm}}(m)=\pi T\left(2m+1+KA(m){\text{sign}}(2m+1)\right)

(28)

where

$\displaystyle A(m)$	$\displaystyle=$	$\displaystyle 2\sum_{1}^{m}\frac{1}{n^{\gamma}},~{}~{}m>0$
$\displaystyle A(m)$	$\displaystyle=$	$\displaystyle A(-m-1),~{}~{}m<-1$
$\displaystyle A(0)$	$\displaystyle=$	$\displaystyle A(-1)=0$	(29)

and

K=\left(\frac{{\bar{g}}}{2\pi T}\right)^{\gamma}

(30)

The expression for $A(m)$ is the same for $\gamma<1$ and $\gamma>1$ . The distinction is in that for $\gamma<1$ , $A(m)\propto m^{1-\gamma}$ , and for $\gamma>1$ , $A(m)$ tends to finite value at $m\to\infty$ : $A(m\to\infty)=2\zeta(\gamma)$ , where $\zeta(\gamma)$ is a Zeta-function. Substituting the self-energy into the equation for $\Phi(\omega_{m})=\Phi(m)$ and eliminating the term with $m=0$ , we obtain

	$\displaystyle\Phi(m>0)$	$\displaystyle=$	$\displaystyle\sum_{n=1,n\neq m}^{\infty}\frac{\Phi(n)}{A(n)+\frac{2n+1}{K}}\frac{1}{\|n-m\|^{\gamma}}+\sum_{n=1}^{\infty}\frac{\Phi(n)}{A(n)+\frac{2n+1}{K}}\frac{1}{(n+m+1)^{\gamma}}-$		(31)
			$\displaystyle\frac{K}{K-1}\sum_{n=1}^{\infty}\frac{\Phi(n)}{A(n)+\frac{2n+1}{K}}\left(\frac{1}{n^{\gamma}}+\frac{1}{(n+1)^{\gamma}}\right)\left(\frac{1}{m^{\gamma}}+\frac{1}{(m+1)^{\gamma}}\right)$		(31)

At small $T$ , when $K\gg 1$ , we obtain from (31):

	$\displaystyle\Phi(m>0)$	$\displaystyle=$	$\displaystyle\sum_{n=1,n\neq m}^{\infty}\frac{\Phi(n)}{A(n)}\frac{1}{\|n-m\|^{\gamma}}+\sum_{n=1}^{\infty}\frac{\Phi(n)}{A(n)}\frac{1}{(n+m+1)^{\gamma}}-$		(32)
			$\displaystyle\sum_{n=1}^{\infty}\frac{\Phi(n)}{A(n)}\left(\frac{1}{n^{\gamma}}+\frac{1}{(n+1)^{\gamma}}\right)\left(\frac{1}{m^{\gamma}}+\frac{1}{(m+1)^{\gamma}}\right)$		(32)

For $m\gg 1$ , Eq. (32) reduces to

\Phi(m>0)=\sum_{n=1,n\neq m}^{\infty}\frac{\Phi(n)}{A(n)}\frac{1}{|n-m|^{\gamma}}+\sum_{n=1}^{\infty}\frac{\Phi(n)}{A(n)}\frac{1}{(n+m)^{\gamma}}-\frac{2}{m^{\gamma}}\sum_{n=1}^{\infty}\frac{\Phi(n)}{A(n)}\left(\frac{1}{n^{\gamma}}+\frac{1}{(n+1)^{\gamma}}\right)

(33)

One can easily verify that relevant $n$ in the sums are of order $m$ , are also large. It is tempting to replace the sum by the integral, with the lower limit of order $T$ . However, this can be done only for $\gamma<1$ , when the integral does not diverge. Keeping $\gamma<1$ , replacing the summation by integration, and restoring Matsubara frequencies $\omega_{m}$ instead of Matsubara numbers, we obtain

\Phi(\omega_{m})=\frac{1-\gamma}{2}\int_{-\infty}^{\infty}d\omega^{\prime}_{m}\frac{\Phi(\omega^{\prime}_{m})}{|\omega^{\prime}_{m}|^{1-\gamma}|\omega_{m}-\omega^{\prime}_{m}|^{\gamma}}-\frac{2(1-\gamma)(2\pi T)^{\gamma}}{|\omega_{m}|^{\gamma}}\int^{\infty}_{O(T)}\frac{\Phi(\omega^{\prime}_{m})}{\omega^{\prime}_{m}}

(34)

At $\omega_{m}\gg T$ , the last term is irrelevant, and $\Phi(\omega_{m})$ has the same form as at $T=0$ : $\Phi(\omega_{m})\propto|\omega_{m}|^{-\gamma/2}\cos{(\beta\log{(|\omega_{m}|/{\bar{g}})^{\gamma}}+\phi)}$ . The phase $\phi$ is set by matching this form and $\Phi(\omega_{m})\propto|\omega_{m}|^{-\gamma}$ at $\omega_{m}\sim{\bar{g}}$ . At $\omega_{m}\sim T$ , i.e., at Matsubara numbers $m=O(1)$ , the last term cannot be neglected. However, it vanishes for certain $T$ , then log-oscillating $\Phi(\omega_{m})$ is the solution of the full Eq. (34). Substituting log-oscillating form into the last term we find that it vanishes when

\beta\gamma\log{\frac{T}{\bar{g}}}=n\pi+{\text{const}},~{}~{}n=0,1,2,...

(35)

This yields the set of $T_{p,n}\propto e^{-\pi n/(\beta\gamma)}$ . Because we assumed that $K\gg 1$ , i.e., $2\pi T\ll{\bar{g}}$ , Eq. (35) is, strictly speaking, valid for $n\gg 1$ .

For $\gamma>1$ , one cannot convert the summation in (33) into integration as the integral will be divergent. Instead, we use the fact that $A(\infty)=2\xi(\gamma)$ is now finite and do the following trick:

(i) rewrite the normal state ${\tilde{\Sigma}}_{\text{norm}}(m)$ as

{\tilde{\Sigma}}_{\text{norm}}(m)=\pi T\left(2m+1+KA(\infty)-K(A(\infty)-A(m))\right){\text{sign}}(2m+1)

(36)

(ii) introduce ${\bar{\tilde{\Sigma}}}_{\text{norm}}(m)$ via

{\bar{\tilde{\Sigma}}}_{\text{norm}}(m)={\tilde{\Sigma}}_{\text{norm}}(m)\left(1-\frac{\pi TKA(\infty)}{{\tilde{\Sigma}}_{\text{norm}}(m)}\right)=-\pi TK\left(A(\infty)-A(m)\right)

(37)

(ii) introduce simultaneously

{\bar{\Phi}}(m)=\Phi(m)\left(1-\frac{\pi TKA(\infty)}{{\tilde{\Sigma}}_{\text{norm}}(m)}\right)

(38)

Because $\Phi(m)/{\tilde{\Sigma}}_{\text{norm}}(m)={\bar{\Phi}}(m)/{\bar{\tilde{\Sigma}}}_{\text{norm}}(m)$ , Eq. (33) becomes

$\displaystyle{\bar{\Phi}}(m>0)$	$\displaystyle=$	$\displaystyle\sum_{n=1,n\neq m}^{\infty}\left(\frac{{\bar{\Phi}}(m)}{A(\infty)-A(m)}-\frac{{\bar{\Phi}}(n)}{A(\infty)-A(n)}\right)\frac{1}{\|n-m\|^{\gamma}}$	(39)
		$\displaystyle+\sum_{n=1}^{\infty}\left(\frac{{\bar{\Phi}}(m)}{A(\infty)-A(m)}-\frac{{\bar{\Phi}}(n)}{A(\infty)-A(n)}\right)\frac{1}{(n+m)^{\gamma}}-$
		$\displaystyle\frac{2}{m^{\gamma}}\left(\sum_{n=1}^{\infty}\frac{{\bar{\Phi}}(n)}{A(\infty)-A(n)}\left(\frac{1}{n^{\gamma}}+\frac{1}{(n+1)^{\gamma}}\right)+\frac{{\bar{\Phi}}(m)}{A(\infty)-A(m)}\right)$

Converting the summation over $n$ into integration, we see that the integral is now free from divergencies. Using that at large $m$ , $A(\infty)-A(m)=m^{1-\gamma}/(\gamma-1)$ and replacing $m$ by $\omega_{m}$ , we obtain

	$\displaystyle{\bar{\Phi}}(\omega_{m})$	$\displaystyle=$	$\displaystyle\frac{\gamma-1}{2}\int_{-\infty}^{\infty}d\omega^{\prime}_{m}\left({\bar{\Phi}}(\omega_{m})\|\omega_{m}\|^{\gamma-1}-{\bar{\Phi}}(\omega^{\prime}_{m})\|\omega^{\prime}_{m}\|^{\gamma-1}\right)\frac{1}{\|\omega_{m}-\omega^{\prime}_{m}\|^{\gamma}}-$		(40)
			$\displaystyle 2\left(\frac{2\pi T}{\|\omega_{m}\|}\right)^{\gamma}\left(\int_{O(T)}^{\infty}\frac{{\bar{\Phi}}(\omega^{\prime}_{m})}{\omega^{\prime}_{m}}+\pi T\frac{{\bar{\Phi}}(\omega_{m})}{\omega_{m}}\right)$		(40)

At $\omega_{m}\gg 2\pi T$ , the last term can be neglected, and we obtain

{\bar{\Phi}}(\omega_{m})=\frac{\gamma-1}{2}\int_{-\infty}^{\infty}d\omega^{\prime}_{m}\left({\bar{\Phi}}(\omega_{m})|\omega_{m}|^{\gamma-1}-{\bar{\Phi}}(\omega^{\prime}_{m})|\omega^{\prime}_{m}|^{\gamma-1}\right)\frac{1}{|\omega_{m}-\omega^{\prime}_{m}|^{\gamma}}

(41)

The solution of this equation is the same log-oscillating function $\Phi(\omega_{m})\propto|\omega_{m}|^{-\gamma/2}\cos{(\beta\log{(|\omega_{m}|/{\bar{g}})^{\gamma}}+\phi)}$ as for $\gamma<1$ , and $\beta$ is again determined by $\epsilon_{i\beta}=1$ , where $\epsilon_{i\beta}$ is given by Eq. (12). Like for $\gamma<1$ , the set of $T_{p,n}$ , where Eq. (40) is valid, is determined by the condition that the last term in (40) vanishes. Substituting log-oscillating form of $\Phi(\omega_{m})$ , we find the same condition on $T_{p,n}$ as in Eq. (40): $T_{p,n}\propto e^{-\pi n/(\gamma\beta)}$ . In Fig.5 we show numerical result for $T_{p,n}$ for $\gamma=1.5$ as a function of $n$ . We see that its dependence on $n$ is exponential, like for $\gamma<1$ .

The computation of the prefactor for $T_{p,n}$ requires more efforts, and we didn’t find it analytically. In Fig. 6 we show numerical results for the onset temperatures $T_{p,n}$ for $\gamma$ around one and $n=0,1,2$ . We see that all $T_{p,n}$ evolve smoothly through $\gamma=1$ .

III.2.4 Nonlinear gap equation, finite $T$

We did not attempt to solve the non-linear gap equation at a finite $T<T_{p,n}$ . Given that $\Delta_{n}(\omega_{m})$ with different $n$ are topologically distinct, and that there is a set of $\Delta_{n}(\omega_{m})$ at $T=0$ , we conjecture that the amplitude of $\Delta_{n}(\omega_{m})$ , which emerges at $T_{p,n}$ , increases as $T$ decreases, and at $T=0$ it coincides with the $n$ th solution of the non-linear gap equation. We illustrate this in Fig. 7.

IV Extension to $N\neq 1$

We now extend the model and introduce a parameter $N$ , which controls the relative strength of the interactions in the particle-hole and particle-particle channels. Like we said in the Introduction, we treat $N$ as a continuous variable. With this extension,

	$\displaystyle\Phi(\omega_{m})$	$\displaystyle=$	$\displaystyle\frac{{\bar{g}}^{\gamma}}{N}\pi T\sum_{m^{\prime}\neq 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}+{\bar{g}}^{\gamma}\pi T\sum_{m^{\prime}\neq 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}}$		(42)

and

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

(43)

Note that here, like in earlier papers, we extend Eliashberg equations to $N\neq 1$ after cancelling out the divergent contribution from thermal fluctuations (the $m^{\prime}=m$ term in the sum over Matsubara frequencies). An alternative approach, suggested in Ref. Wang et al. (2018), is to extend to $N\neq 1$ without first subtracting the $m^{\prime}=m$ terms in Eq. (43). In this case, one has to regularize the divergencies in the r.h.s of these equations and also in the gap equation. In general, the contribution from thermal fluctuations has to be computed differently from other terms in the frequency sum because one cannot factorize the momentum integration based on the separation between fast electrons and slow bosons. We refer a reader to RefsWang et al. (2018); Abanov et al. (2019); *Wu_19_1; Chubukov et al. (2020a), where this issue has been addressed in detail.

We now consider how the solutions of the gap equation, $\Delta_{n}(\omega_{m})$ , evolve near $\gamma=1$ . For this we consider separately the cases $\gamma<1$ and $\gamma>1$ .

IV.1 A generic $\gamma<1$

We first briefly summarize the results for a generic $\gamma<1$ (Parts I and II) and then move to $\gamma\to 1$ .

IV.1.1 Linearized gap equation, $T=0$ .

The linearized gap equation at $T=0$ is

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

(44)

or, equivalently,

D_{\infty}(\omega_{m})\omega_{m}\left(1+\lambda\left(\frac{\bar{g}}{|\omega_{m}|}\right)^{\gamma}\right)=\frac{{\bar{g}}^{\gamma}}{2N}\int d\omega^{\prime}_{m}\frac{D_{\infty}(\omega^{\prime}_{m})-D_{\infty}(\omega_{m})}{|\omega_{m}-\omega^{\prime}_{m}|^{\gamma}}{\text{sign}}\omega^{\prime}_{m},

(45)

where $D(\omega_{m})=\Delta(\omega_{m})/\omega_{m}$ and $\lambda=(1-1/N)/(1-\gamma)$ .

A non-zero solution, $\Delta_{\infty}(\omega_{m})$ , exists for $N<N_{cr}$ , where

N_{cr}=\frac{1-\gamma}{2}\frac{\Gamma^{2}(\gamma/2)}{\Gamma(\gamma)}\left(1+\frac{1}{\cos(\pi\gamma/2)}\right),

(46)

We plot $N_{cr}$ vs $\gamma$ in Fig.8.

For all $\gamma<1$ , $N_{cr}>1$ . Similar to the case $N=1$ , $\Delta_{\infty}(\omega_{m})$ undergoes log-oscillations at $\omega_{m}<{\bar{g}}$ : $\Delta_{\infty}(\omega_{m})\propto|\omega_{m}|^{\gamma/2}\cos(\beta_{N}\log{(|\omega_{m}|/{\bar{g}})^{\gamma}}+\phi)$ , where $\beta_{N}$ is the solution of $\epsilon_{i\beta_{N}}=N$ and $\epsilon_{i\beta}$ is given by (12). A non-zero $\beta_{N}$ exists for $N<N_{cr}$ . For $N\lesssim N_{cr}$ , $\beta_{N}\propto(N_{cr}-N)^{1/2}$ .

IV.1.2 Linearized gap equation, $T\neq 0$ .

At a finite $T$ , the solution of the linearized gap equation exists for a set of critical temperatures, $T_{p,n}$ , like for $N=1$ . An eigenfunction $\Delta_{n}(\omega_{m})$ changes sign $n$ times as a function of discrete Matsubara frequency $\omega_{m}=\pi T(2m+1)$ . All critical lines $T_{p,n}(N)$ for $n>0$ terminate at $T=0$ at $N=N_{cr}$ , while $T_{p,0}$ scales as $1/N^{1/\gamma}$ for large $N$ (Ref. Wang et al. (2016)). At $N=O(1)$ , $T_{p,n}\propto e^{-An}$ for large $n$ .

IV.2 The limit $\gamma\to 1$

IV.2.1 Linearized gap equation, $T=0$ .

In the limit $\gamma\to 1$ , $N_{cr}=1+(\pi/2+\log(4))(1-\gamma)+O(1-\gamma)^{2}$ tends to $1$ , i.e., relevant $N<N_{cr}$ become $N\leq 1$ . Simultaneously, the function $\epsilon_{i\beta}$ becomes flat:

\epsilon_{i\beta_{N}}\approx N_{cr}-(1-\gamma)R(\beta_{N})=1+\left(\frac{\pi}{2}+\log(4)-R(\beta_{N})\right),

where

R(\beta_{N})=\frac{1}{2}\left(\Psi{(1/2+i\beta_{N})}+\Psi{(1/2-i\beta_{N})}\right)-\Psi(1/2)-\frac{\pi}{2}\left(1-\frac{1}{\cosh{\pi\beta_{N}}}\right)

(47)

We plot $R(\beta_{N})$ in Fig. 9. At large $\beta_{N}$ , $R(\beta_{N})\approx\log{|\beta_{N}|}$ . Because $\epsilon_{i\beta_{N}}$ becomes flat, $\beta_{N}$ remains finite for $N=1$ , but exponentially grows for any $N<1$ and becomes infinite at $\gamma=1$ . This implies that the system behavior at $N=1$ and $N<1$ changes discontinuously at $\gamma=1$ . To understand this change, it is instructive to consider the double limit when both $\gamma$ and $N$ tend to one, and $\beta_{N}$ is a continuous function of the ratio $(1-N)/(1-\gamma)$ , or, equivalently, of $(N_{cr}-N)/(1-\gamma)=R(\beta_{N})$ . For $N=N_{cr}$ , $\beta_{N_{cr}}=0$ , for $N=1$ , $\beta_{N=1}=\beta$ tends to $0.792$ , and for $1-N\gg(1-\gamma)$ , $\beta_{N}\approx 0.561/N^{1/(1-\gamma)}\gg 1$ . The case $N<1$ and $\gamma\to 1$ corresponds to the limit $\beta_{N}\rightarrow\infty$ . We emphasize that a continuous evolution is only possible if we keep $N$ as a continuous parameter.

The exact solution for $\Delta_{\infty}(\omega_{m})$ can be obtained for any $\beta_{N}$ . We plot $\Delta_{\infty}(\omega_{m})$ for different $\beta_{N}$ in Fig. 10. To demonstrate the behavior over a large range of frequencies, we use $\log(\omega_{m}/{\bar{g}})$ as a variable. We see that for $\beta_{N}=O(1)$ , $\Delta_{\infty}(\omega_{m})$ oscillates on the logarithmic scale for $\omega_{m}\leq{\bar{g}}$ and decreases as $1/|\omega_{m}|$ at larger frequencies. This agrees with our earlier result for $N\equiv 1$ . However, as $\beta_{N}$ increases, new non-logarithmic oscillations develop at $\omega_{m}\geq 1$ and extend up $\omega_{max}$ . Numerical results strongly indicate that for large enough $\beta_{N}$ , $\omega_{max}\sim{\bar{g}}\log{\beta_{N}}$ , see Fig.11.

This is expected on general grounds because ${\bar{g}}\log{\beta_{N}}\sim{\bar{g}}(1-N)/(1-\gamma)$ , and the latter is the scale at which divergencies in the gap equation are cut when $N\leq 1$ . We also see from Fig.10(b) the overall magnitude of $\Delta_{\infty}(\omega_{m})$ decreases with ${\omega_{m}}$ , while the period of oscillations increases.

To rationalize this observation we again compute the local series $\Delta_{\infty,L}(\omega_{m})$ . We have

\Delta_{L}(\omega_{m})\propto|\omega_{m}|^{1/2}Re\left[e^{i(\beta_{N}\log{|\omega_{m}|/{\bar{g}}}+\phi)}\sum_{m=0}^{\infty}{\tilde{C}}^{N}_{m}\left(\frac{|\omega_{m}|}{\bar{g}}\right)^{m}\right]

(48)

where ${\tilde{C}}^{N}_{m}=\displaystyle\prod_{m^{\prime}=1}^{m}\frac{1}{{\bar{I}}^{N}_{m^{\prime}}}$ , and

{\bar{I}}^{N}_{m^{\prime}}=\frac{((-1)^{m^{\prime}}-1)\pi}{2\cosh{(\pi\beta_{N})}}-\sum_{p=0}^{m^{\prime}-1}\frac{1}{1/2+i\beta_{N}+p}

(49)

The series again converge absolutely, i.e., $\Delta_{L}(\omega_{m})$ can be obtained for any $\omega_{m}$ by summing up enough terms in the series. In Fig. 12 we show both the exact $\Delta_{\infty}(\omega_{m})$ and $\Delta_{\infty,L}(\omega_{m})$ . We see that they nearly coincide over the full range where $\Delta_{\infty}(\omega_{m})$ oscillates. We can also expand the series in (48) in $1/\beta_{N}$ and obtain the analytical expansion in $\omega_{m}/{\bar{g}}$ for the overall factor of $\Delta_{\infty,L}(\omega)$ and the period of oscillations. To leading order in $\beta_{N}$ we find after straightforward but lengthy calculation:

\Delta_{\infty,L}(\omega_{m})\propto|\omega_{m}|^{1/2}f_{1}\left(\frac{|\omega_{m}|}{\bar{g}}\right)\cos{\left(\beta_{N}\left(\log{|\omega_{m}|/{\bar{g}}}+f_{2}\left(\frac{|\omega_{m}|}{\bar{g}}\right)+\phi\right)\right)}

(50)

where

\displaystyle f_{1}(x)=1-\frac{x}{2}+\frac{x^{2}}{8}+....,\qquad f_{2}(x)=-x+\frac{x^{2}}{4}+....

(51)

We see that the envelop of $\Delta_{\infty,L}(\omega_{m})$ varies at $|\omega_{m}|=O({\bar{g}})$ , while the argument of $\cos[...]$ deviates from the low-frequency form $\beta_{N}\log{|\omega_{m}|/{\bar{g}}}+\phi$ already at much smaller $|\omega_{m}|\sim{\bar{g}}/\beta_{N}$ .

In Fig.13 we compare $\Delta_{\infty}(\omega_{m})$ with Eq. (50). We see that the envelope of $\Delta_{\infty}(\omega_{m})$ is well described by Eq. (50), while oscillations become non-logarithmic already at small $|\omega_{m}|\sim{\bar{g}}/\beta_{N}$ and are captured by Eq. (50) up to $\omega_{m}\leq\bar{g}/\beta^{1/3}_{N}$ .

IV.2.2 Non-linear gap equation, $T=0$ .

From a generic point of view, the behavior of $\Delta_{n}(\omega_{m})$ for $\gamma\leq 1$ qualitatively similar to that for smaller $\gamma$ . Namely, $\Delta_{n}(\omega_{m})$ form a discrete, infinite set. A function $\Delta_{n}(\omega_{m})$ behaves as $1/|\omega_{m}|^{\gamma}$ at the highest frequencies, oscillates $n$ times at smaller $\omega_{m}$ , and at even smaller $\omega_{m}$ approaches a finite $\Delta_{n}(0)$ . The condensation energy $E_{c,n}$ is different for different $n$ and is the largest for $n=0$ .

On a more careful look, we find that new features in $\Delta_{n}(\omega_{m})$ gradually develop as $\gamma$ approaches one. To see this, consider the non-linear gap equation at $N\neq 1$ :

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

(52)

For $\gamma\to 1$ and $N\neq 1$ , the dominant contribution to the r.h.s. of (52) comes from $\omega^{\prime}_{m}\approx\omega_{m}$ . Keeping only this contribution, we obtain

\Delta^{2}_{n}(\omega_{m})+\omega^{2}_{m}=\Delta^{2}(0)

(53)

where $\Delta(0)\approx{\bar{g}}^{\gamma}((1-N)/(N(1-\gamma))\sim\omega_{max}$ .

We see that at $\omega_{m}\ll\Delta(0)$ , $\Delta_{n}(\omega_{m})\approx\Delta(0)$ is nearly independent on frequency and is also independent on $n$ . The corrections to Eq. 53 do depend on $n$ , but these corrections are small in $(1-\gamma)/(1-N)\sim{\bar{g}}/\omega_{max}$ . At $\omega_{m}\geq\Delta(0)$ , $\Delta_{n}(\omega_{m})$ oscillates $n$ times and then decreases as $1/|\omega_{m}|$ . Because in this range $\Delta_{n}(\omega_{m})<\omega_{m}$ , the oscillating term is the same as for $\Delta_{\infty,L}(\omega_{m})$ in Eq. (50). A simple analysis then shows that the relative width of the frequency range for $n$ oscillations compare to $\Delta(0)$ is $n/n_{max}$ , where $n_{max}\sim\beta_{N}$ up to a prefactor, which depends on $\log\beta_{N}$ . As long as $n<n_{max}$ , this width is smaller that $\omega_{max}$ , although the upper boundary of oscillations increases with $n$ . We illustrate this in Fig.14

As a result, the range, where $\Delta_{n}(\omega_{m})$ oscillates, accounts only for a subleading contribution to the condensation energy $E_{c,n}$ , the leading one comes from frequencies $\omega_{m}\leq\omega_{max}$ , where $\Delta_{n}\sim\omega_{max}$ is independent on $n$ , up to corrections of order $1/\log{\beta_{N}}\sim(1-\gamma)/(1-N)$ . At $\gamma\to 1$ and $N<1$ , $n_{max}$ tends to infinity and the ratio $E_{c,n}/E_{c,0}$ tends to one for all finite $n$ . At $n\to\infty$ , the result for the condensation energy depends on how the limit $n\to\infty$ and $n_{max}\to\infty$ is taken. If $b=n_{max}/n$ is large, $E_{c,n}\approx E_{c,0}$ . In the opposite limit $b\ll 1$ , oscillations start at a frequency much smaller that $\omega_{max}$ and run up to $\omega_{max}$ . In this case, the corrections to Eq. (53) are no longer small, and the analysis has to be modified. Obviously, at such large $n$ , $\Delta_{n}(0)$ become smaller, and $E_{c,n}$ drops. The outcome of this consideration is that at $\gamma\to 1$ , the spectrum of the condensation energy becomes continuous: $E_{c,n}$ for all finite $n$ becomes equal to $E_{c,0}$ , while at $n\to\infty$ , $E_{c,n}=E_{c}(b)=E_{c,\infty}f(b)$ becomes a continuous variable, ranging between $f(0)=0$ and $f(\infty)=1$ . We illustrate this in Fig.15

There is a certain similarity between our case and how a continuum spectrum develops for lattice vibrations, when the system size becomes infinite and a momentum becomes a continuous variable.

The transformation of the spectrum of $E_{c,n}$ from a discrete one to continuous represents the major qualitative difference between $\gamma=1$ and $\gamma<1$ . For $\gamma<1$ , the set of $E_{c,n}$ is discrete and $E_{c,0}$ is the largest. Because the condensation energy is proportional to the total number of particles, other $E_{c,n}$ are only relevant for spatially inhomogeneous fluctuations at a finite $T$ . When the spectrum of $E_{c,n}$ becomes a continuous one for spatially homogeneous $\Delta(\omega_{m})$ , fluctuations become one-dimensional and will likely restore $U(1)$ phase symmetry. This issue requires further study because although the spectrum becomes continuous, $E_{c,0}$ by itself tends to infinity at $\gamma=1$ . We show in a subsequent publication that a continuous spectrum of condensation energies develops also in more physically transparent case of $N=1$ and $\gamma=2$ . In this last case, $E_{c,0}$ remains finite.

IV.2.3 Linearized gap equation, finite $T$

We recall that for $N=1$ , system behavior evolves smoothly through $\gamma=1$ . Namely, the onset temperature $T_{p,n}$ is of order $\Delta_{n}(0)$ at $T=0$ , and both scale as ${\bar{g}}e^{-An}$ , where $A\sim 1/\beta_{N=1}$ remains $O(1)$ for $\gamma=1$ . An eigenfunction $\Delta_{n}(\omega_{m})$ has the same structure as $\Delta_{\infty}(\omega_{m})$ down to $\omega_{m}\sim T_{p,n}$ , and saturates at smaller frequencies. For $N<1$ , we find two new features, both are consistent with the results at $T=0$ . First, the scale, up to which $\Delta_{n}(\omega_{m})$ oscillates, increases with $n$ (see Fig.16(a)) Second, $T_{p,n}$ at large $n$ still scales as ${\bar{g}}e^{-An}$ , but $A\propto 1/\beta_{N}$ increases as $\gamma\to 1$ , hence $T_{p,n}$ also increases (see Fig.16(b)). We show more numerical results later, in Sec. IV.3.3 and Appendix B.

IV.2.4 Non-linear gap equation, finite $T$

We didn’t attempt to solve the non-linear equation for $\Delta_{n}(\omega_{m})$ at $T<T_{p,n}$ . Like for $N=1$ , we expect, based on the analysis in Sec. IV.2.2, that all $\Delta_{n}(\omega_{m})$ with finite $n>0$ rapidly increase below $T_{p,n}$ and at $T\to 0$ merge with $\Delta_{0}(\omega_{m})$ , which, we recall, develops at a much larger $T_{p,0}$ . We illustrate this in Fig. 21.

IV.3 Case $\gamma\geq 1$

IV.3.1 Linearized gap equation, $T=0$

For $\gamma>1$ , a simple analysis of the linearized gap equation (IV) shows that there is no solution with $\Delta_{n}(\omega_{m})\neq 0$ . Indeed, for $N\neq 1$ , the integral over $\omega^{\prime}_{m}$ diverges at $\omega^{\prime}_{m}=\omega_{m}$ , leaving $\Delta(\omega_{m})=0$ as the only option.

IV.3.2 Nonlinear gap equation, $T=0$

The solution of the non-linear gap equation does not exist for $N>1$ and is singular for $N<1$ . Namely, all $\Delta_{n}(\omega_{m})$ with finite $n$ tend to infinity at any finite $\omega_{m}$ , while the solutions with $n\to\infty$ form a continuous spectrum of the condensation energies. The way to see this is to consider the $T=0$ case as the limit $T\to 0$ . This is what we do in Sec. IV.3.4 below.

IV.3.3 Linearized gap equation, $T\neq 0$

At a finite $T$ the sum over $m^{\prime}\neq m$ in (44) does not diverge. In this situation, it is natural to expect that $\Delta_{n}(\omega_{m})$ is non-zero and finite below a certain $T_{p,n}$ , which, we recall, remains finite for $\gamma=1$ and $N<1$ .

Like for $N=1$ , the calculations are more straightforward, when done for the pairing vertex $\Phi(\omega_{m})$ , expressed via the normal state ${\tilde{\Sigma}}_{\text{norm}}(\omega_{m})$ . The gap function $\Delta(\omega_{m})=\omega_{m}\Phi(\omega_{m})/{\tilde{\Sigma}}_{\text{norm}}(\omega_{m})$ . We have from (IV)

\Phi(\omega_{m})=\frac{{\bar{g}}^{\gamma}}{N}\pi T\sum_{m^{\prime}\neq m}\frac{\Phi(\omega_{m^{\prime}})}{|{\tilde{\Sigma}}_{\text{norm}}(\omega_{m^{\prime}})|}~{}\frac{1}{|\omega_{m}-\omega_{m^{\prime}}|^{\gamma}},

(54)

where ${\tilde{\Sigma}}_{\text{norm}}(\omega_{m^{\prime}})$ is given by Eq. (28) and, we remind, $K=({\bar{g}}/(2\pi T))^{\gamma}$ . For $\gamma>1$ , $A(m\to\infty)=2\zeta(\gamma)$ . Substituting the self-energy into the equation for $\Phi(\omega_{m})=\Phi(m)$ and eliminating the term with $m=0$ , we obtain

	$\displaystyle N\Phi(m>0)$	$\displaystyle=$	$\displaystyle\sum_{n=1,n\neq m}^{\infty}\frac{\Phi(n)}{A(n)+\frac{2n+1}{K}}\frac{1}{\|n-m\|^{\gamma}}+\sum_{n=1}^{\infty}\frac{\Phi(n)}{A(n)+\frac{2n+1}{K}}\frac{1}{(n+m+1)^{\gamma}}+$		(55)
			$\displaystyle\frac{K}{N-K}\sum_{n=1}^{\infty}\frac{\Phi(n)}{A(n)+\frac{2n+1}{K}}\left(\frac{1}{n^{\gamma}}+\frac{1}{(n+1)^{\gamma}}\right)\left(\frac{1}{m^{\gamma}}+\frac{1}{(m+1)^{\gamma}}\right)$		(55)

Consider the limit of small $T$ , when $K\gg 1$ . Like for $\gamma<1$ (Paper II), one solution of (55) exists for $N\approx K\gg 1$ . We express $N$ as $N=K+b_{\gamma}$ , where $b_{\gamma}=O(1)$ and substitute into (55). The divergent $K$ cancels out, and in the remaining equation the kernel factorizes between internal $m^{\prime}$ and external $m$ . Solving the equation, we obtain

	$\displaystyle\Phi(m>0)$	$\displaystyle=$	$\displaystyle C\left(\frac{1}{m^{\gamma}}+\frac{1}{(m+1)^{\gamma}}\right),$
	$\displaystyle\Phi(0)$	$\displaystyle=$	$\displaystyle\frac{1}{b_{\gamma}}\sum_{n=1}^{\infty}\frac{\Phi(n)}{A(n)}\left(\frac{1}{n^{\gamma}}+\frac{1}{(n+1)^{\gamma}}\right)$		(56)

and

b_{\gamma}=\sum_{n=1}^{\infty}\frac{1}{A(n)}\left(\frac{1}{n^{\gamma}}+\frac{1}{(n+1)^{\gamma}}\right)^{2}

(57)

Both $\Phi(m)$ and $b_{\gamma}$ evolve smoothly through $\gamma=1$ . The pairing vertex $\Phi(m)$ and the gap function $\Delta(m)$ do not have nodes and in our classification are $\Phi_{n=0}(m)$ and $\Delta_{n=0}(m)$ . The corresponding $T_{p,0}\approx({\bar{g}}/2\pi)/N^{1/\gamma}$ . We discussed the $n=0$ solution for $N\gg 1$ in length in earlier papers Abanov et al. (2019); *Wu_19_1; Chubukov et al. (2020a). In short: for both $\gamma<1$ and $\gamma>1$ , $\Delta_{n=0}(m)$ displays a re-entrant behavior, i.e., it emerges at a finite $T_{p,0}$ and vanishes at $T=0$ . We verified that for $\gamma>1$ this behavior holds for all $N>1$ .

We now turn to $N<1$ , where, as we will see, system behavior differs qualitatively between $\gamma<1$ and $\gamma>1$ . For $N<1$ and $K\to\infty$ , we obtain from (55):

	$\displaystyle N\Phi(m>0)$	$\displaystyle=$	$\displaystyle\sum_{n=1,n\neq m}^{\infty}\frac{\Phi(n)}{A(n)}\frac{1}{\|n-m\|^{\gamma}}+\sum_{n=1}^{\infty}\frac{\Phi(n)}{A(n)}\frac{1}{(n+m+1)^{\gamma}}-$		(58)
			$\displaystyle\sum_{n=1}^{\infty}\frac{\Phi(n)}{A(n)}\left(\frac{1}{n^{\gamma}}+\frac{1}{(n+1)^{\gamma}}\right)\left(\frac{1}{m^{\gamma}}+\frac{1}{(m+1)^{\gamma}}\right)$		(58)

For $m,n\gg 1$ , the last term becomes irrelevant, and Eq. (58) reduces to

N\Phi(m>0)=\frac{1}{2\zeta(\gamma)}\sum_{n=1,n\neq m}^{\infty}\left(\frac{\Phi(n)}{|n-m|^{\gamma}}+\frac{\Phi(n)}{(n+m)^{\gamma}}\right)

(59)

We search for the solution in the form

\Phi(m>0)\propto\cos\left({\bar{\beta}}_{N}m+{\bar{\phi}}\right)

(60)

where $0<{\bar{\beta}}_{N}<\pi$ . This corresponds to

\Delta(\omega_{m})\propto\omega_{m}\cos\left(\frac{{\bar{\beta}}_{N}}{2\pi T}\omega_{m}+{\bar{\phi}}\right)

(61)

Substituting (60) into (59), we find after simple algebra that Eq. (59) is satisfied if ${\bar{\beta}}_{N}$ satisfies ${\bar{\epsilon}}_{{\bar{\beta}}_{N}}=N$ , where

{\bar{\epsilon}}_{{\bar{\beta}}}=\frac{1}{\zeta(\gamma)}\sum_{n=1}^{\infty}\frac{\cos({\bar{\beta}}n)}{n^{\gamma}}=\frac{1}{\zeta(\gamma)}{\text{Re}}\left[Li_{\gamma}\left(e^{-i{\bar{\beta}}}\right)\right]

(62)

where $Li$ is the polylogarithm function. In Fig.17(a) we plot ${\bar{\beta}}_{N}$ as a function of $N$ for several $\gamma>1$ .

The solution exists for $N$ between maximal $N_{max}\to 1$ and negative minimal $N_{min}=2^{1-\gamma}-1$ . At $N\to N_{max}$ , ${\bar{\beta}}_{N}\to 0$ at $N=N_{min}$ , ${\bar{\beta}}_{N}=\pi$ . For $\gamma\geq 1$ , ${\bar{\beta}}_{N}\sim(1-N)^{1/(\gamma-1)}$ is small for all $N>0$ . At $N=1$ , ${\bar{\beta}}_{N}$ vanishes. In this case, the solution (log-oscillating function) has to be obtained as in Sec. III.2.3.

Like we did for $\gamma<1$ , we interpret ${\bar{\epsilon}}_{{\bar{\beta}}_{N}}=N$ as the dispersion relation and identify ${\bar{\beta}}_{N}$ with the effective momentum and $N$ with the effective energy. Then one can define the density of eigenvalues (DoE) as

\nu(N)\propto\left.\frac{d{\bar{\beta}}}{d{\bar{\epsilon}}_{{\bar{\beta}}}}\right|_{{\bar{\beta}}={\bar{\beta}}_{N}}.

(63)

We plot this function in Fig.17(b) along with the DoE obtained numerically by solving the full Eq. (58) as an eigenvalue/eigenfunction equation. We see that the analytic and numerical DoE are quite similar. Both show divergence at $\gamma-$ dependent $N_{min}$ and vanish at $N_{max}=1$ as $\nu(N)\propto(1-N)^{\frac{2-\gamma}{1-\gamma}}$ . Note that the behavior of $\nu(N)$ near $N=1$ changes at $\gamma=2$ .

We now use the form of $\Phi(m)$ to obtain $T_{p,n}$ . As before, we use the initially free parameter ${\bar{\phi}}$ to match with $\Phi(m)$ at $m=O(1)$ , and match with the power-law form at $\Sigma(\omega_{m})\sim\omega_{m}$ , i.e., at $m\sim K=({\bar{g}}/2\pi T)^{\gamma}$ . In more precise form, we have $m/K\sim(1-N)$ , where $(1-N)$ appears because the constant term in the self-energy $2\pi TK\zeta(\gamma)$ cancels out for $N=1$ . The matching condition is ${\bar{\beta}}_{N}({\bar{g}}/2\pi T)^{\gamma}(1-N)=n\pi+O(1)$ Solving for $T=T_{p,n}$ , we obtain

T_{p,n}\sim\frac{{\bar{g}}}{2\pi}\left(\frac{{\bar{\beta}}_{N}(1-N)}{n\pi}\right)^{1/\gamma}

(64)

We see that the $n-$ dependence of $T_{p,n}$ is now $1/n^{1/\gamma}$ rather than $e^{-An}$ . This implies that for a given $n$ and $N$ , $T_{p,n}$ rapidly increases as $\gamma$ crosses one. For $\gamma\geq 1$ ,

T_{p,n}\sim{\bar{g}}(1-N)^{\frac{1}{\gamma-1}}\left(\frac{1}{n}\right)^{1/\gamma}

(65)

In Fig.18 we plot numerical results for $T_{p,n}$ as a function of $N$ for a given $n$ for several $\gamma\geq 1$ and as a function of $n$ for given $\gamma$ and $N$ . We see that the agreement is quite good.

We now look at the eigenfunctions $\Phi_{n}(m)$ , or, equivalently, $\Delta_{n}(m)$ . We use $\Delta_{n}(m)$ for easier comparison with the results for $\gamma<1$ . The eigenfunctions behave as $\Delta_{n}(m)\sim m\cos({\bar{\beta}}_{N}m+{\bar{\phi}})$ up to $n-$ dependent $m_{max}(n)\sim({\bar{g}}/(2\pi T_{p,n}))^{\gamma}$ . At larger $m$ , oscillations end and each $\Delta_{n}(m)$ decays as $1/m^{\gamma}$ . Comparing with the form of $\Delta_{n}(m)$ for a generic $\gamma<1$ , we see two key differences. First, for $\gamma>1$ , the period of oscillations is set by $m$ rather than by $\log m$ . Second, for $\gamma<1$ , $m_{max}\propto{\bar{g}}/T_{p,n}$ , hence the frequency, where oscillations end, $\omega_{max}=2\pi Tm_{max}$ , is $O({\bar{g}}$ for all $n$ . For $\gamma>1$ , $\omega_{max}\sim{\bar{g}}({\bar{g}}/T_{p,n})^{\gamma-1}$ becomes n-dependent ( $\omega_{max}=\omega_{max}(n)$ ), i.e., the larger is $n$ , the larger is the range of frequencies where $\Delta_{n}(m)$ oscillates as $\cos{({\bar{\beta}}^{eff}_{N}(\omega_{m}/{\bar{g}})+{\bar{\phi}})}$ , where ${\bar{\beta}}^{eff}_{N}={\bar{\beta}}_{N}{\bar{g}}/(2\pi T_{p,n})\sim n^{1/\gamma}$ . At $n\to\infty$ , oscillations extend to $\omega_{m}\to\infty$ . We earlier found the precursor to this behavior for $\gamma\to 1$ . In Fig. 19 we present numerical results for the eigenfunctions $\Delta_{n}(\omega_{m})$ for two different $\gamma\geq 1$ .

We see that the eigenfunctions indeed oscillate with the period set by $\omega_{m}$ rather than $\log{\omega_{m}}$ , and that as $n$ increases, oscillations extend to larger $\omega_{max}$ . These numerical results confirm that there is indeed a qualitative change of system behavior for $N<1$ between $\gamma<1$ and $\gamma>1$ . We also note that the divergence of ${\bar{\beta}}^{eff}_{N}\propto n^{1/\gamma}$ at $n\to\infty$ is consistent with the divergence of $\beta_{N}$ as $\gamma\to 1$ from below.

The crossover from log-oscillations of $\Delta_{n}(\omega_{m})$ for $\gamma<1$ to oscillations with a period set by $\omega_{m}$ for $\gamma>1$ is sharp at $n\gg 1$ , when $T_{c,n}$ is small and relevant Mastubara numbers are large. For smaller $n$ , the crossover gets smoothen up. In numerical calculations, there is an additional smothering due to sampling of a finite number of Matsubara points. In Appendix B we show the numerical results for the crossover behavior and its dependence on the number of Matsubara points, sampled in numerical calculations.

Finally, in Fig.20

we show the dependence of $T_{p,n=1}$ on $N$ near $N=1$ (or, equivalently, the temperature dependence of the second eigenvalue of the gap equation). The onset temperature $T_{p,1}(N)$ decreases as $N$ approaches one from below, but because $T_{p,1}(N=1)$ is finite, it has to remain finite also for $N>1$ . We see that $T_{p,1}(N)$ continuous, as a function of $N$ , into the range $N>1$ , but then reverses trend, such that smaller $T_{p,1}$ correspond to $N$ closer to $N=1$ . This reentrant behavior is the consequence of the fact that at $T=0$ there is no solution of the linearized gap equation for any $N$ , except for $N=1$ .

IV.3.4 Non-linear gap equation, $T\neq 0$

We analyze non-linear equations for the pairing vertex and the self-energy, Eqs. (IV), at small but finite $T$ . It is convenient to introduce dimensionless ${\bar{\Phi}}=(2\pi T)^{\gamma-1}\Phi/{\bar{g}}^{\gamma}$ and ${\bar{\tilde{\Sigma}}}=(2\pi T)^{\gamma-1}{\tilde{\Sigma}}/{\bar{g}}^{\gamma}$ . In these variables, Eqs. (IV) become, for the $n-th$ solution

	$\displaystyle{\bar{\Phi}}_{n}(m)$	$\displaystyle=$	$\displaystyle\frac{1}{2N}\sum_{m^{\prime}\neq m}\frac{{\bar{\Phi}}_{n}(m^{\prime})}{\sqrt{{\bar{\tilde{\Sigma}}}^{2}_{n}(m^{\prime})+{\bar{\Phi}}^{2}_{n}(m^{\prime})}}~{}\frac{1}{\|m-m^{\prime}\|^{\gamma}},$
	$\displaystyle{\bar{\tilde{\Sigma}}}_{n}(m)$	$\displaystyle=$	$\displaystyle\left(\frac{2\pi T}{{\bar{g}}}\right)^{\gamma}(m+1/2)+\frac{1}{2}\sum_{m^{\prime}\neq m}\frac{{\bar{\tilde{\Sigma}}}_{n}(m^{\prime})}{\sqrt{{\bar{\tilde{\Sigma}}}^{2}_{n}(m^{\prime})+{\bar{\Phi}}^{2}_{n}(m^{\prime})}}~{}\frac{1}{\|m-m^{\prime}\|^{\gamma}}$		(66)

Based on our earlier analysis of the case $\gamma\to 1$ from below, we expect that at small $T<T_{p,n}$ , $\Delta_{n}(m)=\pi T(2m+1){\bar{\Phi}}_{n}(m)/{\bar{\tilde{\Sigma}}}_{n}(m)$ is large and weakly dependent on $n$ , up to large $n$ . This holds if ${\bar{\Phi}}_{n}(m)\gg{\bar{\tilde{\Sigma}}}_{n}(m)$ . Using this inequality, we obtain from (IV.3.4)

{\bar{\Phi}}_{n}(m)\approx\frac{\zeta(\gamma)}{N},~{}~{}{\bar{\tilde{\Sigma}}}_{n}(m)\approx\left(\frac{2\pi T}{\bar{g}}\right)^{\gamma}\frac{m+1/2}{1-N}

(67)

and, hence,

\Delta_{n}(m)={\bar{g}}\left(\frac{\bar{g}}{2\pi T}\right)^{\gamma-1}\frac{1-N}{N}+...

(68)

where dots stand for subleading corrections, which depend on $n$ and $m$ . We see that $\Delta_{n}(m)$ diverges as $1/T^{\gamma-1}$ and all $\Delta_{n}(m)$ merge into the same gap function at $T\to 0$ . This holds for $n$ up to some $n_{max}(T)$ , for which $T_{p,{n_{c}}}\sim T$ . Using $T_{p,n}\propto 1/n^{1/\gamma}$ , we obtain $n_{max}\sim({\bar{g}}/T)^{\gamma}$ . As $T\to 0$ , $n_{max}\to\infty$ , hence $\Delta_{n}(m)$ becomes independent on $n$ for all finite $n$ , despite that $T_{p,n}$ are all different. We note in this regard that at $T\leq T_{p,n}$ , $\Delta_{n}(m)$ is of order ${\bar{g}}n^{(\gamma-1)/\gamma}$ , i.e., it increases below $T_{p,n}$ with a slope, which increases with $n$ . We illustrate this in Fig.21

As the consequence, at $T\to 0$ , the condensation energy $E_{c,n}$ becomes equal for all finite $n$ , as we anticipated in Sec. IV.3.2. The gap functions $\Delta_{n}(m)$ with $n>n_{max}$ do not have a $T$ range to develop into Eq. (68), and have smaller condensation energy at $T\to 0$ . The condensation energy for these solutions depends on $b=n_{max}/n=({\bar{g}}/(n^{1/\gamma}T))^{\gamma}$ . At $T\to 0$ , $n_{max}$ tends to infinity, and at $n\to\infty$ , $b$ becomes a continuous variable. The condensation energy $E_{c}(b)=E_{c}(\infty){\tilde{f}}(b)$ , where ${\tilde{f}}(0)=0$ and ${\tilde{f}}(\infty)=1$ . This is consistent with the results in Sec. IV.2.2 on $T=0$ and $\gamma\to 1$ . The behavior of the condensation energy is illustrated in Fig.15 At small $b$ , we used $\Delta_{n}\propto 1/n^{1/\gamma}$ and the expression for the condensation energy from Refs. Chubukov et al. (2020a); Haslinger and Chubukov (2003); Yuzbashyan et al. and obtained $E_{c}\propto(b)^{2(2-\gamma)/\gamma}$ . For $\gamma\geq 1$ , this reduces to $E_{c}\propto b^{2}$ .

We emphasize again that this behavior is qualitatively different from the one in a non-critical BCS/Eliashberg superconductor, where there are at most a few different solutions of the gap equation for ant given $N$ and from quantum-critical superconductivity for $\gamma<1$ , where there exists an infinite set of gap functions for $N<N_{cr}$ , but the spectrum of the condensation energy is discrete. We also emphasize that this behavior does not extend to $N=1$ , for which a discrete set of $\Delta_{n}(m)$ and $E_{c,n}$ survives for $\gamma>1$ . The difference with $N=1$ is obvious from Eq. (68), which shows that the divergent term cancels out for $N=1$ .

V Another extension of the $\gamma$ model

We now propose another extension of the original model, which does not introduce divergencies. For this we return to the original model with $N=1$ and re-express Eqs. (IV) for the pairing vertex and the self-energy by pulling out the divergent terms from the r.h.s., like we did in Sec. III.2.3. We obtain

	$\displaystyle\Phi(\omega_{m})\left(1-{\bar{g}}^{\gamma}\frac{\zeta(\gamma)}{(2\pi T)^{\gamma-1}}\frac{1}{\sqrt{\tilde{\Sigma}^{2}(\omega_{m})+\Phi^{2}(\omega_{m})}}\right)=$
	$\displaystyle{\bar{g}}^{\gamma}\pi T\sum_{m^{\prime}\neq m}\left(\frac{\Phi(\omega_{m^{\prime}})}{\sqrt{{\tilde{\Sigma}}^{2}(\omega_{m^{\prime}})+\Phi^{2}(\omega_{m^{\prime}})}}-\frac{\Phi(\omega_{m})}{\sqrt{{\tilde{\Sigma}}^{2}(\omega_{m})+\Phi^{2}(\omega_{m})}}\right)~{}\frac{1}{\|\omega_{m}-\omega_{m^{\prime}}\|^{\gamma}},$
	$\displaystyle{\tilde{\Sigma}}(\omega_{m})\left(1-{\bar{g}}^{\gamma}\frac{\zeta(\gamma)}{(2\pi T)^{\gamma-1}}\frac{1}{\sqrt{\tilde{\Sigma}^{2}(\omega_{m})+\Phi^{2}(\omega_{m})}}\right)=\omega_{m}+$
	$\displaystyle{\bar{g}}^{\gamma}\pi T\sum_{m^{\prime}\neq m}\left(\frac{{\tilde{\Sigma}}(\omega_{m^{\prime}})}{\sqrt{{\tilde{\Sigma}}^{2}(\omega_{m^{\prime}})+\Phi^{2}(\omega_{m^{\prime}})}}-\frac{{\tilde{\Sigma}}(\omega_{m})}{\sqrt{{\tilde{\Sigma}}^{2}(\omega_{m})+\Phi^{2}(\omega_{m})}}\right)~{}\frac{1}{\|\omega_{m}-\omega_{m^{\prime}}\|^{\gamma}}$		(69)

We then introduce

	$\displaystyle{\bar{\Phi}}(\omega_{m})$	$\displaystyle=$	$\displaystyle\Phi(\omega_{m})\left(1-{\bar{g}}^{\gamma}\frac{\zeta(\gamma)}{(2\pi T)^{\gamma-1}}\frac{1}{\sqrt{\tilde{\Sigma}^{2}(\omega_{m})+\Phi^{2}(\omega_{m})}}\right)$
	$\displaystyle{\bar{\tilde{\Sigma}}}(\omega_{m})$	$\displaystyle=$	$\displaystyle{\tilde{\Sigma}}(\omega_{m})\left(1-{\bar{g}}^{\gamma}\frac{\zeta(\gamma)}{(2\pi T)^{\gamma-1}}\frac{1}{\sqrt{\tilde{\Sigma}^{2}(\omega_{m})+\Phi^{2}(\omega_{m})}}\right)$		(70)

Because $\Phi(\omega_{m})/{\tilde{\Sigma}}(\omega_{m})={\bar{\Phi}}(\omega_{m})/{\bar{\tilde{\Sigma}}}(\omega_{m})$ , Eqs. (V) can be re-expressed solely in terms of ${\bar{\Phi}}(\omega_{m})$ and ${\bar{\tilde{\Sigma}}}(\omega_{m})$ :

	$\displaystyle{\bar{\Phi}}(\omega_{m})={\bar{g}}^{\gamma}\pi T\sum_{m^{\prime}\neq m}\left(\frac{{\bar{\Phi}}(\omega_{m^{\prime}})}{\sqrt{{\bar{\tilde{\Sigma}}}^{2}(\omega_{m^{\prime}})+{\bar{\Phi}}^{2}(\omega_{m^{\prime}})}}-\frac{{\bar{\Phi}}(\omega_{m})}{\sqrt{{\bar{\tilde{\Sigma}}}^{2}(\omega_{m})+{\bar{\Phi}}^{2}(\omega_{m})}}\right)~{}\frac{1}{\|\omega_{m}-\omega_{m^{\prime}}\|^{\gamma}},$
	$\displaystyle{\bar{\tilde{\Sigma}}}(\omega_{m})=\omega_{m}$
	$\displaystyle+{\bar{g}}^{\gamma}\pi T\sum_{m^{\prime}\neq m}\left(\frac{{\bar{\tilde{\Sigma}}}(\omega_{m^{\prime}})}{\sqrt{{\bar{\tilde{\Sigma}}}^{2}(\omega_{m^{\prime}})+{\bar{\Phi}}^{2}(\omega_{m^{\prime}})}}-\frac{{\bar{\tilde{\Sigma}}}(\omega_{m})}{\sqrt{{\bar{\tilde{\Sigma}}}^{2}(\omega_{m})+{\bar{\Phi}}^{2}(\omega_{m})}}\right)~{}\frac{1}{\|\omega_{m}-\omega_{m^{\prime}}\|^{\gamma}}$		(71)

These equations are now free from singularities, even if we replace a summation over Mascara numbers by an integration over Matsubara frequencies.

We now extend the modified Eliashberg equations (V) in the same way as before, by multiplying the interaction in the particle-particle channel by a factor $1/M$ :

{\bar{\Phi}}(\omega_{m})=\frac{{\bar{g}}^{\gamma}}{M}\pi T\sum_{m^{\prime}\neq m}\left(\frac{{\bar{\Phi}}(\omega_{m^{\prime}})}{\sqrt{{\bar{\tilde{\Sigma}}}^{2}(\omega_{m^{\prime}})+{\bar{\Phi}}^{2}(\omega_{m^{\prime}})}}-\frac{{\bar{\Phi}}(\omega_{m})}{\sqrt{{\bar{\tilde{\Sigma}}}^{2}(\omega_{m})+{\bar{\Phi}}^{2}(\omega_{m})}}\right)~{}\frac{1}{|\omega_{m}-\omega_{m^{\prime}}|^{\gamma}}

(72)

The gap function $\Delta(\omega_{m})$ is expressed via ${\bar{\Phi}}(\omega_{m})$ and ${\bar{\tilde{\Sigma}}}(\omega_{m})$ in the same way as via the original $\Phi(\omega_{m})$ and ${\tilde{\Sigma}}(\omega_{m})$ : $\Delta(\omega_{m})=\omega_{m}{\bar{\Phi}}(\omega_{m})/{\bar{\tilde{\Sigma}}}(\omega_{m})$ . The gap equation becomes

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

(73)

V.1 Linearized gap equation, $T=0$

For infinitesimally small ${\bar{\Phi}}(\omega_{m})$ , the self-energy coincides with that in the normal state. Converting $\pi T\sum_{m^{\prime}}$ into $\int d\omega^{\prime}_{m}/2$ and evaluating the frequency integral in (V), we obtain at $T=0$ ,

{\bar{\tilde{\Sigma}}}(\omega_{m})=\omega_{m}B(|\omega_{m}|)

(74)

where

B(|\omega_{m}|)=1-\left(\frac{\bar{g}}{|\omega_{m}|}\right)^{\gamma}\frac{1}{\gamma-1}

(75)

Substituting ${\bar{\tilde{\Sigma}}}(\omega_{m})$ into the equation for ${\bar{\Phi}}(\omega_{m})$ , we obtain

{\bar{\Phi}}_{\infty}(\omega_{m})=\frac{{\bar{g}}^{\gamma}}{2M}\int d\omega^{\prime}_{m}\left(\frac{{\bar{\Phi}}_{\infty}(\omega^{\prime}_{m})}{|\omega^{\prime}_{m}|B(|\omega^{\prime}_{m}|)|}-\frac{{\bar{\Phi}}_{\infty}(\omega_{m})}{|\omega_{m}|B(|\omega_{m}|)|}\right)\frac{1}{|\omega^{\prime}_{m}-\omega_{m}|^{\gamma}}

(76)

we label infinitesimally small ${\bar{\Phi}}(\omega_{m})$ as ${\bar{\Phi}}_{\infty}(\omega_{m})$ , like in earlier analysis, anticipating that the non-linear equation for the pairing vertex will have a discrete set of solutions ${\bar{\Phi}}_{n}(\omega_{m})$ . At small $\omega_{m}$ , the bare $\omega_{m}$ term in ${\bar{\tilde{\Sigma}}}(\omega_{m})$ can be neglected, and (76) reduces to

{\bar{\Phi}}_{\infty}(\omega_{m})=\frac{\gamma-1}{2M}\int d\omega^{\prime}_{m}\frac{{\bar{\Phi}}_{\infty}(\omega_{m})|\omega_{m}|^{\gamma-1}-{\bar{\Phi}}_{\infty}(\omega^{\prime}_{m})|\omega^{\prime}_{m}|^{\gamma-1}}{|\omega^{\prime}_{m}-\omega_{m}|^{\gamma}}

(77)

This equation is similar, but not equivalent, to Eq. (7) for the pairing vertex for $\gamma<1$ . In both cases, the kernel is marginal, and we search for the solution in the form ${\bar{\Phi}}(\omega_{m})\propto|\omega_{m}|^{-\gamma/2+b}$ . Like before, a normalizable solution exists when $b=\pm i\beta_{M}$ is purely imaginary. Substituting power-law form with the complex exponent, we find that (77) is satisfied if $\epsilon_{i\beta_{M}}=M$ , where $\epsilon_{i\beta}$ is exactly the same function as in Eq. (12). In this respect, the extension to $M\neq 1$ for $\gamma>1$ is quite similar to the extension to $N\neq 1$ for $\gamma<1$ . The similarity is particularly transparent for the linearized equation for $D(\omega_{m})=\Delta(\omega_{m})/\omega_{m}$ . From (77) we obtain

D_{\infty}(\omega_{m})\omega_{m}\left(1+{\bar{\lambda}}\left(\frac{\bar{g}}{|\omega_{m}|}\right)^{\gamma}\right)=\frac{{\bar{g}}^{\gamma}}{2M}\int d\omega^{\prime}_{m}\frac{D_{\infty}(\omega^{\prime}_{m})-D_{\infty}(\omega_{m})}{|\omega_{m}-\omega^{\prime}_{m}|^{\gamma}}{\text{sign}}\omega^{\prime}_{m},

(78)

where ${\bar{\lambda}}=(1/M-1)/(\gamma-1)$ . This equation is identical to Eq. (45) once we replace $N$ by $M$ .

Because $d\epsilon_{i\beta}/d\beta$ is positive for $\gamma>1$ , a normalizable solution of the gap equation exists for $M>M_{cr}$ , where $M_{cr}$ satisfies $\epsilon_{0}=M_{cr}$ . We plot $M_{cr}$ as a function of $\gamma$ in Fig.22. We see that $M_{cr}$ gradually decreases as $\gamma$ increases, from $M_{cr}=1$ at $\gamma=1$ to $M_{cr}=0$ at $\gamma=2$ .

V.2 Linearized gap equation at a finite $T$

At a finite $T$ , we obtain from (73) for infinitesimally small $\Delta_{\infty}(\omega_{m})$ = $\Delta_{\infty}(m)$ and $m>0$

S_{m}\left(|2m+1|-2K(\zeta(\gamma)-H(m,\gamma))\right)=\frac{K}{M}\sum_{m^{\prime}\neq m}\frac{S_{m^{\prime}}-S_{m}}{|m^{\prime}-m|^{\gamma}}

(79)

where, $S_{m}=\Delta(m)/|2m+1|$ , $K=({\bar{g}}/(2\pi T))^{\gamma}$ , and $H(m,\gamma)=\sum_{1}^{m}1/p^{\gamma}$ is the the Harmonic number. At small $T$ , (79) simplifies to

S_{m}=\frac{1}{2M(H(m,\gamma)-\zeta(\gamma))}\sum_{m^{\prime}\neq m}\frac{S_{m}-S_{m^{\prime}}}{|m^{\prime}-m|^{\gamma}}

(80)

The solution of this equation exists at a particular $T$ , which determines the onset temperature for the pairing.

We show results of the numerical solution of Eq. (79) in Fig.23. Like for $\gamma<1$ , we find that there exists a discrete set of onset temperatures $T_{p,n}$ , and an eigenfunction $\Delta_{n}(m)$ changes sign $n$ times as a function of Matsubara number $m$ . Different $T_{p,n}(M)$ all approach $M=M_{cr}$ at $T=0$ , although for larger $n$ one has to go to very low $T$ to see this. Such behavior is similar to the one for $\gamma<1$ , the only distinction is that now there is no special behavior of $T_{p,0}$ because ${\bar{\tilde{\Sigma}}}(\omega_{m})$ does not vanish at the first two Matsubara frequencies $\omega_{m}=\pm\pi T$ .

Because the behavior of the gap function is the same in all models with $M>M_{cr}$ , including the original model with $M=1$ , the extension to $M\neq 1$ allows one to extract this behavior by focusing at either $M\geq M_{cr}$ , where the analytical analysis is simplified because relevant frequencies are small, or at $M>1$ , where $T_{p,n}$ are larger and can be detected more easily in numerical studies.

VI Conclusions

In this paper, we continued our analysis of the interplay between the pairing and the non-Fermi liquid behavior in a metal for a set of quantum-critical systems, which at low-energies are described by a model of fermions with an effective dynamical electron-electron interaction $V(\Omega_{m})\propto 1/|\Omega_{m}|^{\gamma}$ between fermions at the Fermi surface (the $\gamma$ -model). We analyzed both the original model and its extension, in which we introduce an extra parameter $N$ to account for non-equal interactions in the particle-hole and particle-particle channel. In the two previous papers we considered the case $0<\gamma<1$ and argued that (i) at $T=0$ , there exists an infinite discrete set of topologically different gap functions, $\Delta_{n}(\omega_{m})$ , all with the same spatial symmetry, and (ii) each $\Delta_{n}$ evolves with temperature and terminates at a particular $T_{p,n}$ . In this paper we analyze how the system behavior changes between $\gamma<1$ and $\gamma>1$ , both at $T=0$ and a finite $T$ . We show that the limit $\gamma\to 1$ is singular due to infra-red divergence of $\int d\omega_{m}V(\Omega_{m})$ , and the system behavior is highly sensitive to how this limit is taken. We showed that in the original model with $N=1$ the divergencies cancel out in the gap equation, and the gap functions $\Delta_{n}(\omega_{m})$ smoothly evolve through $\gamma=1$ both at $T=0$ and a finite $T$ . However, for $N\neq 1$ , the evolution through $\gamma=1$ is not smooth, and qualitatively new behavior emerges for $\gamma\geq 1$ . Namely, there still exists a discrete set of $T_{p,n}$ , below which $\Delta_{n}(\omega_{m})$ appears, but (i) the functional forms of $T_{p,n}$ and $\Delta_{n}(\omega_{m})$ change qualitatively, and (ii) at $T\to 0$ all $\Delta_{n}(\omega_{m})$ with $n<n_{max}\sim({\bar{g}}/T)^{\gamma}$ tend to the same gap function. At $T\to 0$ , $n_{max}$ tends to infinity, and the spectrum of the condensation energy $E_{c,n}$ becomes a continuous one. This opens up the new channel of one-dimensional gap fluctuations. We also discussed another extension of the $\gamma$ -model for $\gamma>1$ , to $M\neq 1$ , for which the extended model is free from singularities, and displays the same behavior as the original model with $M=1$ . This allows one to better understand the physics of the original model by zooming into ranges of $M$ where either analytical or numerical analysis is simplified.

In the next paper in the series, Paper IV, we consider the original $\gamma$ -model in the range $1<\gamma<2$ in more detail. We argue that dynamical vortices appear one-by-one in the upper half-plane of frequency as $\gamma$ increases between one and two, and the new physics emerges at a finite $T$ . In Paper V we show that for $\gamma=2$ , the number of these vortices becomes infinite, and the new physics extends down to $T=0$ .

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 AVC was supported by the NSF DMR-1834856.

Appendix A The exact solution of the linearized gap equation for $\gamma<1$ .

In Paper I we obtained the exact solution of the linearized equation for $T=0$ , $\gamma<1$ and any $N$ . The solution has the form:

\Delta_{\infty}(\omega_{m})\propto|\omega_{m}|^{-\gamma}\tilde{f}(\log|\omega_{m}/\omega_{0}|^{\gamma})

(81)

where $\omega_{0}={\bar{g}}/(1-\gamma)^{1/\gamma}$ and

\tilde{f}(x)=\int_{-\infty}^{\infty}\tilde{b}(\beta)e^{-i\beta x}d\beta,

(82)

where

\tilde{b}(\beta)=\frac{\sinh(\pi\beta_{N})}{\sqrt{\cosh(\pi(\beta-\beta_{N})\cosh(\pi(\beta+\beta_{N}))}}e^{-iI(\beta)}.

(83)

Here, $\beta_{N}$ is the solution of $\epsilon_{i\beta_{N}}=N$ , $\epsilon_{i\beta}$ is given by (12), and

I(\beta)=\frac{1}{2}\int_{-\infty}^{\infty}\log|1-\frac{1}{N}\epsilon_{i\beta^{\prime}}|\left[\tanh(\pi(\beta^{\prime}-\beta))-\tanh(\pi\beta^{\prime})\right]d\beta^{\prime}

(84)

Notice that $I(\beta)$ is real and antisymmetric.

In the limit $\gamma\to 1$ ,

\epsilon_{i\beta}=N_{cr}+(1-\gamma)R(\beta),

where $R(\beta)$ is given by (47). Then

1-\frac{\epsilon_{i\beta}}{N}\approx(1-\gamma)\left(\frac{N-N_{cr}}{1-\gamma}-R(\beta)\right)=(1-\gamma)\left(R(\beta_{N})-R(\beta)\right),

(85)

and the function $I(\beta)$ in (84) becomes

I(\beta)=-\beta\log(1-\gamma)+J(\beta)

(86)

where

J(\beta)\equiv\frac{1}{2}\int_{-\infty}^{\infty}\log|R(\beta_{N})-R(\beta^{\prime})|\left[\tanh(\pi(\beta^{\prime}-\beta))-\tanh(\pi\beta^{\prime})\right]d\beta^{\prime}

(87)

does not depend on $\gamma$ . We see that the function $\tilde{f}(x)$ in (82) is in fact the function of $x-\log(1-\gamma)$ . Substituting into (81) we find that at $\gamma\to 1$ , the argument of $\tilde{f}$ is $\log{|\omega_{m}|/\omega_{0}}-\log{(1-\gamma)}=\log{|\omega_{m}|/{\bar{g}}}$ , i.e., relevant scale for $\Delta_{\infty}(\omega_{m})$ is ${\bar{g}}$ rather than the divergent $\omega_{0}$ :

\Delta_{\infty}(\omega_{m})=|\omega_{m}|^{-1}\tilde{f}(\log|\omega_{m}|/\bar{g}).

Appendix B Numerical results for the crossover from logarithmic to power-law oscillations

In this Appendix, we present the results of a detailed numerical analysis of the crossover from log-oscillations of $\Delta_{n}(\omega_{m})$ for $\gamma<1$ to oscillations with a period set by $\omega_{m}$ for $\gamma>1$ . Like we said in Sec. IV.3.3, the transformation at $\gamma=1$ is sharp at $n\gg 1$ , when $T_{c,n}$ is small and relevant Matsubara numbers are large. For smaller $n$ , the crossover gets smoothen up. In numerical calculations, there is an additional smothering due to sampling of a finite number of Matsubara points.

We show the results in Figs. 24 and 25. We see from Fig. 24 that as $\gamma$ approaches one, log-oscillations of $\Delta_{n}(\omega_{m})$ at a given $n$ progressively get replaced by power-law oscillations. The period of power-law oscillations and the envelope are best fitted by $(\omega_{m}/{\bar{g}})^{t}$ and $(\omega_{m}/{\bar{g}})^{r}$ , respectively. For an infinite number of sampling points, we expect a sharp crossover at $\gamma=1$ and $n\to\infty$ between logarithmic and power-law oscillations.

In Fig.25 we show the results for the exponents $r$ and $t$ , extracted from Fig. 24. For a given $\gamma$ , the values of $r$ and $t$ vary with $n$ and temperature. Analytically, we obtained in Sec. IV.3.3 $r=t=1$ for large $n$ , when $T_{p,n}\ll{\bar{g}}$ . We see that $r$ and $t$ are different from one for a given $n$ , but both tend to $1$ when $n$ becomes large enough and $T\to 0$ . This result confirms our analysis in Eq.(60) for $\gamma\geq 1$ . Overall, the numerical results clearly show the main result of our analysis — the transformation of log-oscillations for $\gamma<1$ into power-law oscillations for $\gamma>1$ .

References

Abanov and Chubukov (2020) A. Abanov and A. V. Chubukov, “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$ ,” (2020), arXiv:2004.13220 [cond-mat.str-el] .
Wu et al. (2020) Y. Wu, A. Abanov, Y. Wang, and A. V. Chubukov, “Interplay between superconductivity and non-Fermi liquid at a quantum critical point in a metal. II. The $\gamma$ -model at a finite $T$ for $0<\gamma<1$ ,” (2020), arXiv:2006.02968 [cond-mat.supr-con] .
Wang et al. (2016) Y. Wang, A. Abanov, B. L. Altshuler, E. A. Yuzbashyan, and A. V. Chubukov, Phys. Rev. Lett. 117, 157001 (2016).
Chubukov et al. (2020a) A. V. Chubukov, A. Abanov, Y. Wang, and Y.-M. Wu, Annals of Physics , 168142 (2020a).
Abanov et al. (2019) A. Abanov, Y.-M. Wu, Y. Wang, and A. V. Chubukov, Phys. Rev. B 99, 180506 (2019).
Wu et al. (2019a) Y.-M. Wu, A. Abanov, Y. Wang, and A. V. Chubukov, Phys. Rev. B 99, 144512 (2019a).
Raghu et al. (2015) S. Raghu, G. Torroba, and H. Wang, Phys. Rev. B 92, 205104 (2015).
Wang et al. (2017) H. Wang, S. Raghu, and G. Torroba, Phys. Rev. B 95, 165137 (2017).
Wang et al. (2018) H. Wang, Y. Wang, and G. Torroba, Phys. Rev. B 97, 054502 (2018).
Abanov et al. (2001) A. Abanov, A. V. Chubukov, and A. M. Finkel’stein, EPL (Europhysics Letters) 54, 488 (2001).
Abanov et al. (2003) A. Abanov, A. V. Chubukov, and J. Schmalian, Advances in Physics 52, 119 (2003).
Moon and Chubukov (2010) E.-G. Moon and A. Chubukov, Journal of Low Temperature Physics 161, 263 (2010).
Metlitski and Sachdev (2010) M. A. Metlitski and S. Sachdev, Phys. Rev. B 82, 075127 (2010).
Mross et al. (2010) D. F. Mross, J. McGreevy, H. Liu, and T. Senthil, Phys. Rev. B 82, 045121 (2010).
Monthoux et al. (2007) P. Monthoux, D. Pines, and G. G. Lonzarich, Nature 450, 1177 (2007).
Efetov et al. (2013) K. B. Efetov, H. Meier, and C. Pepin, Nature Physics 9, 442 (2013).
Metlitski et al. (2015) M. A. Metlitski, D. F. Mross, S. Sachdev, and T. Senthil, Phys. Rev. B 91, 115111 (2015).
Lee et al. (2018) T.-H. Lee, A. Chubukov, H. Miao, and G. Kotliar, Phys. Rev. Lett. 1805, 10280 (2018).
Chubukov et al. (2020b) A. V. Chubukov, A. Abanov, I. Eserlis, and S. A. Kivelson, Annals of Physics (2020b).
Abrikosov et al. (1965) A. A. Abrikosov, L. P. Gorkov, and I. E. Dzyaloshinski, Methods of Quantum Feld Theory in Statistical Physics (Pergamon Press, Oxford, 1965).
Wu et al. (2019b) Y.-M. Wu, A. Abanov, and A. V. Chubukov, Phys. Rev. B 99, 014502 (2019b).
Haslinger and Chubukov (2003) R. Haslinger and A. V. Chubukov, Phys. Rev. B 68, 214508 (2003).
(23) E. A. Yuzbashyan, A. V. Chubukov, A. Abanov, and B. L. Altshuler, In preparation.

Interplay between superconductivity and non-Fermi liquid at a quantum-critical point in a metal. III: The γ\gamma model and its phase diagram across γ=1\gamma=1.

Abstract

I Introduction.

II γ\gamma-model, Eliashberg equations

III Transformation from γ≤1\gamma\leq 1 to γ≥1\gamma\geq 1

III.1 A generic γ<1\gamma<1

III.2 γ≈1\gamma\approx 1

III.2.1 Linearized gap equation, T=0T=0

III.2.2 Nonlinear gap equation, T=0T=0

III.2.3 Linearized gap equation, finite TT

III.2.4 Nonlinear gap equation, finite TT

IV Extension to N≠1N\neq 1

IV.1 A generic γ<1\gamma<1

IV.1.1 Linearized gap equation, T=0T=0.

IV.1.2 Linearized gap equation, T≠0T\neq 0.

IV.2 The limit γ→1\gamma\to 1

IV.2.1 Linearized gap equation, T=0T=0.

IV.2.2 Non-linear gap equation, T=0T=0.

IV.2.3 Linearized gap equation, finite TT

IV.2.4 Non-linear gap equation, finite TT

IV.3 Case γ≥1\gamma\geq 1

IV.3.1 Linearized gap equation, T=0T=0

IV.3.2 Nonlinear gap equation, T=0T=0

IV.3.3 Linearized gap equation, T≠0T\neq 0

IV.3.4 Non-linear gap equation, T≠0T\neq 0

V Another extension of the γ\gamma model

V.1 Linearized gap equation, T=0T=0

V.2 Linearized gap equation at a finite TT

VI Conclusions

Acknowledgements.

Appendix A The exact solution of the linearized gap equation for γ<1\gamma<1.

Appendix B Numerical results for the crossover from logarithmic to power-law oscillations

References

Interplay between superconductivity and non-Fermi liquid at a quantum-critical point in a metal.
III: The $\gamma$ model and its phase diagram across $\gamma=1$ .

II $\gamma$ -model, Eliashberg equations

III Transformation from $\gamma\leq 1$ to $\gamma\geq 1$

III.1 A generic $\gamma<1$

III.2 $\gamma\approx 1$

III.2.1 Linearized gap equation, $T=0$

III.2.2 Nonlinear gap equation, $T=0$

III.2.3 Linearized gap equation, finite $T$

III.2.4 Nonlinear gap equation, finite $T$

IV Extension to $N\neq 1$

IV.1 A generic $\gamma<1$

IV.1.1 Linearized gap equation, $T=0$ .

IV.1.2 Linearized gap equation, $T\neq 0$ .

IV.2 The limit $\gamma\to 1$

IV.2.1 Linearized gap equation, $T=0$ .

IV.2.2 Non-linear gap equation, $T=0$ .

IV.2.3 Linearized gap equation, finite $T$

IV.2.4 Non-linear gap equation, finite $T$

IV.3 Case $\gamma\geq 1$

IV.3.1 Linearized gap equation, $T=0$

IV.3.2 Nonlinear gap equation, $T=0$

IV.3.3 Linearized gap equation, $T\neq 0$

IV.3.4 Non-linear gap equation, $T\neq 0$

V Another extension of the $\gamma$ model

V.1 Linearized gap equation, $T=0$

V.2 Linearized gap equation at a finite $T$

Appendix A The exact solution of the linearized gap equation for $\gamma<1$ .