Interplay between superconductivity and non-Fermi liquid at a quantum-critical point in a metal: V. The γ model and its phase diagram. The case γ=2
Yi-Ming Wu, Shang-Shun Zhang, Artem Abanov, Andrey V. Chubukov
IInterplay between superconductivity and non-Fermi liquid at aquantum-critical point in a metal.V. The γ model and its phase diagram. The case γ = 2 . Yi-Ming Wu, Shang-Shun Zhang, Artem Abanov, and Andrey V. Chubukov School of Physics and Astronomy and William I. Fine Theoretical Physics Institute,University of Minnesota, Minneapolis, MN 55455, USA Department of Physics, Texas A&M University, College Station, USA (Dated: February 9, 2021) a r X i v : . [ c ond - m a t . s up r- c on ] F e b bstract This paper is a continuation and a partial summary of our analysis of the pairing at a quantum-critical point (QCP) in a metal for a set of quantum-critical systems, whose low-energy physics isdescribed by an effective model with dynamical electron-electron interaction V (Ω m ) = (¯ g/ | Ω m | ) γ (the γ -model). Examples include pairing at the onset of various spin and charge density-wave andnematic orders and pairing in SYK-type models. In previous papers, we analyzed the physics for γ <
2. We have shown that the onset temperature for the pairing T p is finite, of order ¯ g , yet thegap equation at T = 0 has an infinite set of solutions within the same spatial symmetry. As theconsequence, the condensation energy E c has an infinite number of minima. The spectrum of E c is discrete, but becomes more dense as γ increases. Here we consider the case γ = 2. The γ = 2model attracted special interest in the past as it describes the pairing by an Einstein phonon inthe limit when the dressed phonon mass ω D vanishes. We show that for γ = 2, the spectrum of E c becomes continuous. We argue that the associated gapless ”longitudinal” fluctuations destroysuperconducting phase coherence at a finite T , such that at 0 < T < T p the system displayspseudogap behavior. We show that for each gap function from the continuum spectrum, there isan infinite array of dynamical vortices in the upper half-plane of frequency. For the electron-phononcase at a finite ω D , our results show that T p = 0 . g , obtained in earlier studies, marks the onsetof the pseudogap behavior, while the actual superconducting T c scales with ω D and vanishes at ω D → I. INTRODUCTION.
This work presents a continuation and a partial summary of our analysis of the inter-play between non-Fermi liquid (NFL) physics and superconductivity for a set of quantum-critical (QC) itinerant systems, whose low-energy dynamics can be described by an effectivemodel of spin-full electrons with dynamical interaction on the Matsubara axis V (Ω m ) =(¯ g/ | Ω m | ) γ (the γ -model). Examples include pairing near spin-density-wave, charge-densitywave, and Ising-nematic instabilities in isotropic and anisotropic 3D and 2D systems, pair-ing of fermions at a half-filled Landau level, pairing of dispersion-less fermions randomlycoupled to phonons, and so on. We listed the examples in the first paper in the series, Ref.[1], and discussed earlier works. In that and subsequent papers [1–4], hereafter called Papers2-IV, we analyzed the γ -model with exponents 0 < γ <
2. We rationalized Eliashberg-typeapproach, solved generalized Eliashberg equations, and found that the solution with a non-zero pairing gap develops below a finite T p , which for a generic γ = O (1) is of order ¯ g . Thecorresponding gap function ∆( k, ω m ) can be roughly approximated as ∆( ω m ) f ( k ), wherea normalized f ( k ) has a particular spatial symmetry ( d − wave, s + − , etc), and ∆( ω m ) is asign-preserving function of Matsubara frequency, whose amplitude increases with decreasing T . At T = 0, ∆(0) is of order T p (the ratio ∆(0) /T p is a γ − dependent number [5, 6]). In thisrespect, the pairing at a QCP is similar to pairing away from a QCP, when V (Ω m ) saturatesat a finite value at Ω m = 0. However, on a more careful look, we found qualitative differencebetween the two cases. Namely, away from a QCP, Eliashberg gap equation at T = 0 has atmost a finite number of solutions with a given spatial symmetry. At a QCP, it has an infinitenumber of solutions for ∆( ω ), i.e., the condensation energy, E c , has an infinite number oflocal minima. The solutions ∆ n ( ω m ), labeled by an integer n , are topologically distinct inthe sense that ∆ n ( ω m ) changes sign n times along the positive Matsubara axis, and eachsuch point is a vortex in the upper half-plane of frequency.The ultimate goal of our studies of the γ − model is to understand how the existence ofan infinite set of solutions affects the interplay between pairing (i.e., the appearance of anon-zero ∆( ω )) and a true superconductivity. In a conventional Eliashberg theory (ET) ofSC out of a non-critical Fermi liquid, phase fluctuations are small in the same parameterthat allows one to neglect vertex corrections. In this situation, the onset temperature forthe pairing, T p , and the actual superconducting T c nearly coincide. To address a possiblereduction of superconducting T c by phase fluctuations, one then has to include the effectsnot incorporated into ET, e.g., a localization of electrons when the interaction exceeds thefermionic bandwidth. We analyze whether the emergence of an infinite number of solutionsfor the gap at a QCP gives rise to a substantial reduction of T c /T p ratio already when theinteraction is smaller than the fermionic bandwidth. For γ <
2, which we analyzed in PapersI-IV, the spectrum of the condensation energies, E c,n , is discrete, and E c, is the largest. Inthis situation, the physics at small T is determined by a single solution ∆ ( ω m ), and phasefluctuations are weak. We showed, however, that the set becomes more dense at γ increases:the ratios ( E c,n +1 − E c,n ) /E c,n get progressively smaller.In this paper, we analyze the case γ = 2. We argue that for this γ , the set of the gapfunctions and the spectrum of the condensation energies become continuous. Specifically, at3 = 2 − δ , and δ = 0+, ∆ n ( ω m ) with n < /δ become equal to ∆ ( ω m ) for all ω m >
0, and E c,n become equal to E c, , while ∆ n ( ω m ) and E c,n with infinite n > /δ form a continuous,one-parameter gapless spectra, ∆ ξ ( ω m ) E c,ξ , where ξ , that runs between zero and infinity,is a continuous variable, which depends on how the double limit n → ∞ and δ → ξ such that the minimum of E c,ξ is at ξ = 0, and E c, ∞ = 0). We argue thatfluctuation corrections to superconducting order parameter from E c,ξ destroy long-rangesuperconducting order at any finite T . We emphasize that this holds for itinerant fermions,in the limit when the interaction is smaller than the bandwidth, and the ET is rigorouslyjustified.We present a corroborative evidence that the γ = 2 model is critical. It comes fromthe analysis of the gap equation on the real frequency axis and in the upper half-plane offrequency, z = ω (cid:48) + iω (cid:48)(cid:48) , ω (cid:48)(cid:48) >
0. The gap function ∆( z ) generally cannot be obtained from∆( ω m ) by just replacing iω m by z as such gap function is not guaranteed to be analytic. Toobtain an analytic function, one has to perform a more sophisticated analysis [7–10]. As aconsequence, ∆( ω ) on the real axis can be quite different from ∆( ω m ). For the γ -model,some difference is expected on general grounds, particularly for γ >
1, because while theinteraction V (Ω m ) on the Matsubara axis is positive (attractive) for all γ , the one on thereal axis is complex, V (Ω) = e iπγ/ / | Ω | γ , and its real part V (cid:48) (Ω) = (¯ g/ | Ω | ) γ cos πγ/ γ > n ( ω m ) and ∆ n ( ω ) for 1 < γ <
2. We foundthat at ω, ω m < ¯ g and at ω, ω m > ¯ g ( | log (2 − γ ) | / (2 − γ )) / , the two gap functions trans-fer into each other under a rotation iω m → ω . However, at intermediate ¯ g < ω, ω m < ( | log (2 − γ ) | / (2 − γ )) / , ∆ n ( ω m ) = a n / | ω m | γ is a sign-preserving function of frequency,while ∆ n ( ω ) = | ∆ n ( ω ) | e iη n ( ω ) oscillates, and its phase η n ( ω ) winds up by 2 πk γ , where k γ isan integer, which depends on γ , but not on n . We extended the analysis to complex z in theupper frequency half-plane and showed that there exists an array of k dynamical vortices,centered at some complex z i .Here we show that for γ = 2, the gap functions on the real axis form a continuous set,and each ∆ ξ ( ω ) oscillates up to an infinite frequency, and its phase winds up by an infinitenumber of 2 π . Accordingly, the number of vortices at z i becomes infinite, and the array of z i extends up to an infinite frequency, where, we argue, each ∆ ξ ( ω ) develops an essentialsingularity. We show that for each gap function from the set, the density of states (DoS)4 ξ ( ω + i
0) has an infinite number of maxima and minima, and does not recover the normalstate form up to ω = ∞ . For the solution with ξ = 0, which was studied before [7–9, 11, 12], N ( ω + iδ ) reduces to a set of δ − functions at some ω i .We combine the results for γ = 2 and earlier results for γ < γ − model for γ ≤
2, Fig.16. For all γ , the ground state is asuperconductor with a finite superfluid stiffness ρ s , and the onset temperature for the pairing, T p , is finite. However, superconducting T c decreases with γ and vanishes for γ = 2. Inbetween T p and T c , the system displays a pseudogap (preformed pairs) behavior. Onefeature of this phase is ”gap filling” behavior, as T increases towards T p . In the next paperwe consider the case γ >
2. We show that the behavior at a finite T remains largely thesame as for γ = 2, however new physics emerges ar T = 0 and gives rise to a reduction andeventual vanishing of ρ s even in the ground state.The model with the pairing interaction V (Ω m ) = (¯ g/ | Ω m | ) attracted a substantial at-tention on its own as it describes the pairing, mediated by an Einstein boson, in the limitwhere the effective (dressed) Debye frequency ω D vanishes. [13]. Electron-phonon model at ω D → λ = (¯ g/ω D ) diverges at ω D →
0. However, the interaction ¯ g is still assumed to be much smaller thanthe Fermi energy, E F . Indeed, ET includes contributions to all orders in λ within the ladderapproximation, but neglects vertex corrections to ladder series. The latter hold in powersof Migdal-Eliashberg parameter λ E = ¯ g N /ω D = λ ( N ω D ), where N ∼ /E F in the DoSper unit volume. For small enough ¯ g/E F , λ E remains small even when λ is large. From thisperspective, the strong coupling limit of the ET is the double limit in which ω D and ¯ g/E F tend to zero simultaneously, such that λ E remains small. In physical terms, the smallnessof λ E (cid:28) λ comes about because in a process that gives rise to a vertex correction, fermionsare forced to vibrate at a phonon frequency, far away from their own resonance, while in theprocesses, which form series in λ , fermions are vibrating near their resonance frequencies.The smallness of λ E also allows one to neglect the renormalization of the bosonic propagatorby fermions, both in the normal and in the superconducting state.Previous studies have found that a non-zero gap function emerges at T p ≈ . ω D e − /λ
5t weak coupling (Refs. [18–24]) and at T p = 0 . g at strong coupling [8, 9, 25, 26]. [27].To understand the interplay between the onset of pairing and T c , one has to also computesuperfluid stiffness, ρ s . At weak coupling, ρ s ∼ E F (cid:29) T p [28–33]. In this situation, T c and T p almost coincide. At strong coupling, the situation is more complex. At T = 0, the ρ s ∼ T p /λ E (see below). Within the validity of ET, this stiffness exceeds T p . If we wereto neglect the continuum spectrum of the condensation energy, we would obtain that T p and T c again also coincide, as thermal corrections to SC order parameter are of order T /ρ s and remain small for all T < T p . Including the additional corrections from the continuumof E c,ξ , we find that thermal corrections actually hold in powers of T / ( ω D λ E ) and becomeof order one at T ∼ ¯ g N ∼ ω D λ E , which we identify with the actual T c . At small ω D / ¯ g ,this T c is much smaller than T p , even if we set λ E = O (1). In between T = T p and T c thesystem displays a preformed pairs behavior. When ω D increases and becomes of order ¯ g ,the pseudogap region shrinks and the system gradually recovers BCS-like behavior (Fig.15).The structure of the paper is the following. In Sec.II we present the Eliashberg gapequations that we use in this paper. In Sec. III we discuss the solution of the gap equationalong the Matsubara axis at T = 0 and γ →
2. We first obtain, in Secs. III A and III B,the exact solution of the linearized gap equation, ∆ ∞ ( ω m ), which changes sign an infinitenumber of times between ω m = 0 and ω m ∼ ¯ g , and sign-preserving solution ∆ ( ω m ), whichtends to a finite value at ω m →
0. At larger ω m > ¯ g , both ∆ ∞ ( ω m ) and ∆ ( ω m ) scaleas 1 / | ω m | . In Sec. III C, we obtain the solutions of the non-linear gap equation in theorder-by-order expansion in the gap magnitude and show that they form a one-parametercontinuum set ∆ ξ ( ω ), for which ∆ ∞ ( ω m ) and ∆ ( ω m ) are the two limiting cases. In Sec.IV, we analyze the properties of the gap function ∆( ω ) along the real frequency axis. Wefirst obtain, in Sec. IV A, the exact solution of the linearized gap equation on the realaxis, ∆ ∞ ( ω ), and show that it oscillates not only at ω < ¯ g , but also at ω > ¯ g , with adifferent period. In Sec. IV B we consider the real-frequency form of ∆ ( ω ), which does notchange sign on the Matsubara axis. We use as an input the results from earlier works [7–9], which demonstrated that ∆( ω ) = | ∆( ω ) | e iη ( ω ) oscillates at ω > ¯ g and argue that thephase η ( ω ) winds up by an infinite number of 2 π between ω = O (¯ g ) and ω = ∞ . In Sec.IV C we present a one-parameter continuum set of ∆ ξ ( ω ), which in the two limits reducesto ∆ ∞ ( ω ) and ∆ ( ω ), respectively. In Sec. V we extend ∆ ξ ( ω ) into the upper frequencyhalf-plane ( ω → z ) and show that for each ξ , there is an infinite array of vortices in the6pper frequency half-plane and an essential singularity at | z | = ∞ . In Sec. VI we considerthe gap equation at a finite ω D . We argue that the number of vortices becomes finite andthe high-frequency behavior of the gap function becomes regular; however, this holds onlyabove a frequency, which scales inversely with ω D . In Sec. VII we consider fluctuationcorrections to superconducting order parameter ∆ (cid:104) e iη (cid:105) We argue that the ground state is asuperconductor, however corrections to (cid:104) e iη (cid:105) become O (1) already at T ≤ ω D . We identifythis scale with the actual superconducting T c and discus pseudogap behavior in between T p ∼ ¯ g and T c . In Sec. VIII we combine the results for γ = 2 and for γ < γ model for γ ≤
2. We present our conclusionsin Sec. IX. Some technical aspects are discussed in the Appendices.
II. ELIASHBERG EQUATIONS
As we said, Eliashberg equations for γ = 2 are the same as for the interaction withan Einstein phonon in the limit of vanishing ω D . The electron-phonon problem has beenanalyzed by many authors [7–9, 11, 12, 14–17, 26, 34–36] and we will be brief. On theMatsubara axis, the equation for the gap function ∆ n ( ω m ) is∆( ω m ) = ¯ g πT (cid:88) m (cid:48) ∆( ω m (cid:48) ) − ∆( ω m ) ω m (cid:48) ω m (cid:112) ( ω m (cid:48) ) + ∆ ( ω m (cid:48) ) 1 | ω m − ω m (cid:48) | + ω D . (1)The self-action term with m (cid:48) = m in the r.h.s. can be safely eliminated because the numer-ator vanishes at m = m (cid:48) . Taking then the limit ω D →
0, we obtain∆( ω m ) = ¯ g πT (cid:88) m (cid:48) (cid:54) = m ∆( ω m (cid:48) ) − ∆( ω m ) ω m (cid:48) ω m (cid:112) ( ω m (cid:48) ) + ∆ ( ω m (cid:48) ) 1 | ω m − ω m (cid:48) | . (2)At T = 0, πT (cid:80) m (cid:48) (cid:54) = m → (1 / (cid:82) dω m .The gap equation on the real axis is more conveniently expressed in terms of D ( ω ) = ∆( ω ) /ω. (3)The equation has the form [9] D ( ω ) ωB ( ω ) = A ( ω ) + C ( ω ) (4)7here A ( ω ) = − ¯ g (cid:90) ∞ dω (cid:48) (cid:60) D ( ω (cid:48) ) (cid:112) − D ( ω (cid:48) ) A T A T = tanh ω (cid:48) T + tanh ω T ( ω (cid:48) + ω ) + tanh ω (cid:48) T − tanh ω T ( ω (cid:48) − ω ) − T cosh ω T ω (cid:48) ( ω (cid:48) ) − ω B ( ω ) = 1 + ¯ g ω (cid:90) ∞ dω (cid:48) (cid:34) (cid:60) (cid:112) − D ( ω (cid:48) ) (cid:35) B T B T = tanh ω (cid:48) T + tanh ω T ( ω (cid:48) + ω ) − tanh ω (cid:48) T − tanh ω T ( ω (cid:48) − ω ) + 1 T cosh ω T ω ( ω (cid:48) ) − ω C ( ω ) = − i ¯ g π (cid:112) − D ( ω ) (cid:34) dD ( ω ) dω tanh ω T − T (cid:32) dD ( ω ) dω + (cid:18) dD ( ω ) dω (cid:19) D ( ω )1 − D ( ω ) (cid:33)(cid:35) (5)where the integrals are principal values. At T = 0, the expressions simplify to A ( ω ) = − ¯ g (cid:90) ∞ dω (cid:48) ( | ω | + ω (cid:48) ) (cid:60) D ( ω (cid:48) ) (cid:112) − D ( ω (cid:48) ) B ( ω ) = 1 + ¯ g | ω | (cid:90) ∞ dω (cid:48) ( | ω | + ω (cid:48) ) (cid:60) (cid:112) − D ( ω (cid:48) ) C ( ω ) = − i π ¯ g (cid:112) − D ( ω ) dD ( ω ) dω sign ω, (6)and the gap equation becomes − i π ¯ g dD ( ω ) dω (cid:112) − D ( ω ) sign ω = (7) D ( ω ) ω (cid:32) g | ω | (cid:90) ∞ dω (cid:48) ( | ω | + ω (cid:48) ) (cid:60) (cid:112) − D ( ω (cid:48) ) (cid:33) + (cid:90) ∞ dω (cid:48) ( | ω | + ω (cid:48) ) (cid:60) D ( ω (cid:48) ) (cid:112) − D ( ω (cid:48) )The functions A ( ω ) and B ( ω ) can be equivalently expressed in terms of the solution ofthe gap equation on the Matsubara axis [8, 9, 36]: A ( ω ) = 2 πT ∞ (cid:88) m =0 ∆( ω m ) (cid:112) Ω m + ∆ ( ω m ) ω m − ω ( ω m + ω ) B ( ω ) = 1 + 4 πT ∞ (cid:88) m =0 (cid:112) Ω m + ∆ ( ω m ) ω m ( ω m + ω ) . (8)This simplifies numerical calculations: the recipe is to first solve for the gap at the Matsubatapoints ω m = πT (2 m + 1) and then use Eqs. (8) as an input for the calculation of D ( ω ) onthe real axis. 8 II. GAP EQUATION ALONG THE MATSUBARA AXIS AT T = 0 In Papers I-IV we analyzed the gap equation for γ < T = 0 it hasan infinite, discrete set of solutions at T = 0, ∆ n ( ω m ). A gap function ∆ n ( ω m ) changessign n times between ω m = 0 and ω m = O (¯ g ) and decays as 1 / | ω m | γ at larger frequencies.The two end points of the set are the sign-preserving solution ∆ ( ω m ) and the solution ofthe linearized gap equation ∆ ∞ ( ω m ), which changes sign an infinite number of times. Theexistence of this infinite set is a distinct feature of the pairing at a QCP. Away from a QCP,the number of solutions becomes finite ( n = 0 , ..n max ), and far away from a QCP only the n = 0 solution remains, like in a conventional BCS/Eliashberg theory.Here we extend this analysis to γ = 2. We show that for this γ , the set of gap functionsbecomes continuous ∆ ξ ( ω m ), where 0 ≤ ξ ≤ ∞ is a continuous variable. We first analyzethe two end points, ∆ ∞ ( ω m ) and ∆ ( ω m ), and then obtain the gap function for arbitrary ξ . A. Linearized gap equation
The linearized gap equation at T = 0 is obtained from (2) by assuming that the gapfunction is infinitesimally small. In terms of D ( ω m ) = ∆( ω m ) /ω m we have D ∞ ( ω m ) = ¯ g ω m (cid:90) dω (cid:48) m D ∞ ( ω (cid:48) m ) − D ∞ ( ω m )( ω m − ω m (cid:48) ) sign ω (cid:48) m . (9)One can verify that the leading term in D ∞ ( ω m ) at small ω m (cid:28) ¯ g is obtained by neglectingthe l.h.s. of (9), i.e., by solving (cid:90) ∞ dω (cid:48) m (cid:20) D ∞ ( ω (cid:48) m ) − D ∞ ( ω m )( ω m − ω m (cid:48) ) + D ∞ ( ω (cid:48) m ) + D ∞ ( ω m )( ω m + ω m (cid:48) ) (cid:21) = 0 (10)This approximation is equivalent to neglecting the bare ω in the fermionic propagator incomparison with the NFL fermionic self-energy without the self-action term, Σ( ω ) = − ¯ g /ω .The solution of (10) is D ∞ ( ω m ) = 2 (cid:15) Re (cid:20) e i ( β log ( | ωm | ¯ g ) + φ ) (cid:21) sign ω = 2 (cid:15) cos (cid:32) β log (cid:18) | ω m | ¯ g (cid:19) + φ (cid:33) sign ω (11)where (cid:15) is an infinitesimally small real overall factor, φ is a phase factor, which is arbitraryat this stage, and β = 0 . πβ tanh( πβ ) = 1 and is the solution of (cid:90) ∞−∞ dx | x | iβ − sign x ( x − = 0 , (12)9 - - -
10 20 30 - - - - - - ω m / g D ∞ ( ω m ) ω m D ∞ ( ω m ) FIG. 1. D ∞ ( ω m ) as a function of ω m / ¯ g . The scale is logarithmic for ω m < ¯ g and linear at ω m > ¯ g . The function D ∞ ( ω m ) is scale-invariant (an arbitrary phase factor φ can be absorbed intothe prefactor for ω under the logarithm). This is the consequence of the fact that ¯ g falls offfrom the gap equation (9), once we neglect D ∞ ( ω m ) in the l.h.s.We now analyze the full gap equation. By power counting, the r.h.s of (9) is of order D ∞ ( ω m )(¯ g/ | ω m | ). This justifies neglecting D ∞ ( ω m ) in the l.h.s. for | ω m | < O (¯ g ), but forlarger frequencies it must be kept.We obtained the exact solution of Eq. (9). The derivation parallels the one for smaller γ in Papers I and IV. We skip the details and present the final result: D ∞ ( ω m ) = (cid:15) ¯ gω m (cid:90) ∞−∞ dkb k e − ik log ( ω m / ¯ g ) , (13)where b k = e − iI k [cosh( π ( k − β )) cosh( π ( k + β ))] / , (14)and I k = 12 (cid:90) ∞−∞ dk (cid:48) log | (cid:15) k (cid:48) − | tanh π ( k (cid:48) − k + iδ ) , (15) (cid:15) k (cid:48) = πk (cid:48) tanh( πk (cid:48) ) , (16)Here β (cid:39) . ω m (cid:28) ¯ g , the exact D ∞ ( ω m ) has the form of Eq. (11) with some particular φ . At ω m (cid:29) ¯ g , D ∞ ( ω m ) does not oscillate and decreases as 1 / ( ω m ) (∆ ∞ ( ω m ) decreases as 1 / ( ω m ) ). Weplot the exact D ∞ ( ω m ) in Fig. 1. The crossover between the two forms occurs at ω m ∼ ¯ g ,as expected. 10he corrections to Eq. (11) at small ω m hold in powers of | ω m | / ¯ g ; the leading correctionscales as ( | ω m | / ¯ g ) . . The corrections to 1 / ( ω m ) at large ω m hold in powers of ¯ g/ | ω m | ; theleading correction scales as (¯ g/ | ω m | ) log( | ω m | / ¯ g ). We present the details of the analysis inAppendix B. There, we also show that at ω m (cid:29) ¯ g there also exists an exponentially small,oscillating component D ∞ ; u ( ω m ) in the form D ∞ ; u ( ω m ) ∝ √ (cid:15)e − ( | ω m | / ¯ g ) cos (cid:34) ( π − π (cid:18) | ω m | ¯ g (cid:19) + π (cid:35) (17)where (cid:15) is an infinitesimally small factor. This term is the contribution to D ∞ from large k and k (cid:48) in Eqs. (14) and (15), It is completely irrelevant on the Matsubara axis, but we willsee that it gives the dominant contribution to D ∞ ( ω ) on the real axis. B. Non-linear gap equation. Sign-preserving solution.
We now analyse the full non-linear gap equation, Eq. (2). We first search for a “conven-tional” sign-preserving solution ∆ ( ω m )The analytical analysis uses the same computational steps as in Paper IV and we will bebrief. We use the identity (cid:90) ∞−∞ dω (cid:48) m − ω (cid:48) m ω m | ω m − ω (cid:48) m | γ = 0 , (18)valid for γ >
1, and re-express Eq. (2) as∆ ( ω m ) (cid:34) − ¯ g (cid:90) ∞−∞ dω (cid:48) m − ω (cid:48) m ω m | ω m − ω (cid:48) m | (cid:32) (cid:112) ∆ ( ω (cid:48) m ) + ( ω (cid:48) m ) − ( ω m ) (cid:33)(cid:35) = ¯ g (cid:90) ∞−∞ dω (cid:48) m ∆ ( ω (cid:48) m ) − ∆ ( ω m ) | ω m − ω (cid:48) m | (cid:112) ∆ ( ω (cid:48) m ) + ( ω (cid:48) m ) (19)Both integrals in (19) are infra-red convergent and are determined by ω (cid:48) m ≤ ∆( ω (cid:48) m ). In thelimit ω m →
0, Eq. (19) reduces to∆ (0) (cid:34) − ¯ g (cid:90) ∞ dω (cid:48) m (cid:112) ∆ ( ω (cid:48) m ) + ( ω (cid:48) m ) − ∆ (0)∆ (0) (cid:112) ∆ ( ω (cid:48) m ) + ( ω (cid:48) m ) | ω (cid:48) m | (cid:35) = ¯ g (cid:90) ∞ dω (cid:48) m (∆ ( ω (cid:48) m ) − ∆ (0)) (cid:112) ∆ ( ω (cid:48) m ) + ( ω (cid:48) m ) | ω (cid:48) m | (20)We assume and then verify that ∆ ( ω (cid:48) m ) ≈ ∆ (0) for ω (cid:48) m ≤ ∆ ( ω (cid:48) m ), relevant for bothintegrals in (20). Substituting into (20), we find1 = ¯ g (cid:90) ∞ dω (cid:48) m (cid:112) ∆ (0) + ( ω (cid:48) m ) − ∆ (0)∆ (0) (cid:112) ∆ (0) + ( ω (cid:48) m ) | ω (cid:48) m | (21)11 -1 /g -4 -2 () / g ω / ¯ g Δ ( ω m ) / ¯ g FIG. 2. Sign-preserving solution ∆ ( ω m ) of the nonlinear gap equation along the Matsubara axis.At ω m < ¯ g , ∆ ( ω m ) remains close to ¯ g ; at larger frequencies it decays as 1 /ω m . The integral can be evaluated analytically and yields ∆ (0) = ¯ g . This is quite close to thenumerical value 1 . g (Ref. [8]). Substituting further ∆ ( ω (cid:48) m ) = ¯ g into the r.h.s. of (19),we find that ∆ ( ω m ) varies quadratically with ω m at small ω m and for ω m ≤ ¯ g remains closeto ∆ (0) (Fig.2). In the opposite limit of large ω m , the prefactor for ∆ ( ω ) in the l.h.s. of(19) is approximately 1, and in the r.h.s. of this equation 1 / | ω m | can be pulled out fromthe integral. This yields ∆ ( ω m ) ≈ Q (cid:18) ¯ g | ω m | (cid:19) (22)where Q = (cid:90) ∞ dω (cid:48) m ∆ ( ω (cid:48) m ) (cid:112) ∆ ( ω (cid:48) m ) + ( ω (cid:48) m ) (23)The integral is determined by ω (cid:48) m ∼ ¯ g and is of order ¯ g . Then ∆ ( ω m ) ∼ ¯ g / | ω m | at highfrequencies. The full gap function is sign-preserving (Fig.2) C. Continuous set of solutions. Expansion in the gap magnitude
The solutions ∆ ∞ ( ω m ) and ∆ ( ω m ) (or, equivalently, D ∞ ( ω m ) and D ( ω m )) also exist for γ <
2. For such γ , these two solutions are the end points of a discrete set of topologicallydistinct solutions ∆ n ( ω ). For a continuous set there is no one-to-one correspondence betweena particular member of the set and integer n , and we will show how this correspondence gets12ost at γ = 2 − D ∞ ( ω m ) and D ( ω m ), we see that they have the same form 1 / ( ω m ) for ω m > ¯ g , but are very different for ω m < ¯ g . We therefore focus on the range ω m < ¯ g anduse to our advantage the fact that we know the analytic form of D ∞ ( ω m ) in this range, Eq.(11). We use this D ∞ ( ω m ) as an input and expand it in powers of D ( ω (cid:48) m ) in the r.h.s. ofthe gap equation (2). Specifically, we will be searching for the solution of (2) in the form D ( ω m ) = ∞ (cid:88) j =0 D (2 j +1) ( ω m ) (24)where D (1) ( ω m ) = D ∞ ( ω m ) = 2 (cid:15) cos f ( ω m ) sign ω m (25)with f ( ω m ) = β log (cid:18) | ω m | ¯ g (cid:19) + φ (26)We will see that D (2 j +1) ∼ (cid:15) j +1 .Substituting D ( ω m ) from (24) into (2) and expanding in D ( ω (cid:48) m ) in the r.h.s. of (2), weobtain the set of equations, which express D (2 j +1) for a given j in terms of D (2 j +1) withsmaller j . For j = 1, we have D (3) ( ω m ) ω m − ¯ g (cid:90) dω (cid:48) m (cid:0) D (3) ( ω (cid:48) m ) − D (3) ( ω m ) (cid:1) sign ω (cid:48) m ( ω m − ω (cid:48) m ) = K ( ω m ) (27)where the source term is K ( ω m ) = − ¯ g (cid:90) dω (cid:48) m (cid:0) D (1) ( ω (cid:48) m ) − D (1) ( ω m ) (cid:1) [ D (1) ( ω (cid:48) m )] sign ω (cid:48) m ( ω m − ω (cid:48) m ) (28)The source term is of order (cid:15) , hence D (3) ∝ (cid:15) ( D (5) ∝ (cid:15) and so on). Substituting D (1) ( ω m )from Eq. (25) and evaluating the integrals, we find the source term for D (3) as the sum ofthe two terms, K = K a + K b , where K a ( ω m ) = − (cid:15) ¯ gω m cos (3 f ( ω )) (2 πβ coth(2 πβ ) − πβ tanh(3 πβ )) sign ω m (29)and K b ( ω m ) = − (cid:15) ¯ gω m cos ( f ( ω )) 1 + sinh ( πβ )sinh ( πβ ) sign ω m (30)Solving for D (3) we find that the first term gives rise to (cid:15) cos (3 f ( ω )), while the second termaccounts for the renormalization of the prefactor for log( ω m ) in f ( ω m ) in (26) To order (cid:15) ,13he dressed f ( ω m ), which we label f (cid:15) ( ω m ), becomes f (cid:15) ( ω m ) = β (cid:15) log (cid:18) | ω m | ¯ g (cid:19) + φ (cid:15) (31)where β (cid:15) = β (cid:0) − (cid:15) / (cid:1) ≈ β (1 − (cid:15) ) / (32)The full D ( ω m ) to order (cid:15) is D ( ω m ) = 2 (cid:0) (cid:15) cos f (cid:15) ( ω m ) + Q (cid:15) cos 3 f ( ω m ) (cid:1) sign ω m (33)where Q = 2 πβ coth (2 πβ ) − πβ tanh (3 πβ )2(1 − πβ tanh (3 πβ )) = 5 − ( πβ ) ≈ .
222 (34)Expanding to next order, we find (i) (cid:15) cos 5 f ( ω m ) term with the prefactor Q = 0 . O ( (cid:15) ) corrections to β (cid:15) in (32) ( β (cid:15) = 1 − . (cid:15) + 0 . (cid:15) )), and (iii) O ( (cid:15) ) corrections to Q ( Q → Q ,(cid:15) ) and to the argument of cos 3 f ( ω m ) in (33). We verified that the last correctionchanges cos 3 f ( ω m ) to cos 3 f (cid:15) ( ω m ) with the same f (cid:15) as in (31). This is the strong indicationthat the series contain the same fully renormalized f (cid:15) ( ω m ) in each term. Combining theresults, we obtain, for ω m (cid:28) ¯ g , D ( ω m ) = D (cid:15) ( ω m ), where D (cid:15) ( ω m ) = 2 (cid:15) (cid:0) cos f (cid:15) ( ω m ) + Q ,(cid:15) (cid:15) cos 3 f (cid:15) ( ω m ) + Q ,(cid:15) (cid:15) cos 5 f (cid:15) ( ω m ) + ... (cid:1) sign ω m (35)We emphasize that a continuous set of solutions exists only for γ = 2. Applying the sameperturbative analysis for γ <
2, we find that the expansion holds in (cid:15) (¯ g/ | ω m | ) − γ andbreaks at a finite ω min ∼ ¯ g(cid:15) / (2 − γ ) (see Appendix A for more detail). At smaller ω m , ∆( ω m )saturates, and D ( ω m ) ∝ /ω m . The forms of D ( ω m ) at ω m < ω min and ω m > ω min matchonly for a discrete set of (cid:15) = (cid:15) n , which implies that for γ < f (cid:15) ( ω m ) contains log ω m , each D (cid:15) ( ω m ) from (35) changes sign an infinite numberof times down to ω m = 0, i.e., in our original classification the gap functions from the setare different realizations of n = ∞ . At ω m = 0, each D (cid:15) ( ω m ) has an essential singularity asneither lim ω m → D (cid:15) ( ω m ) nor lim ω m → /D (cid:15) ( ω m ) exist.For a generic (cid:15) , Eq. (35) is valid for ω m < ¯ g . At larger ω m , D (cid:15) ( ω m ) = D (cid:15) / | ω m | . Weexpect that for every (cid:15) , the crossover to proper high-frequency behavior can be achieved byfixing the phase factor φ (cid:15) in (31) (see Paper I for a similar analysis for the linearized gapequation for γ < β (cid:15) decreases with increasing (cid:15) . It is natural to expectthat it vanishes at some (cid:15) cr = O (1). The expansion in (35) holds only as long as β (cid:15) isreal, as there is no solution of the linearized gap equation for imaginary β (cid:15) (see Paper I fordetailed discussion on this). For (cid:15) ≤ (cid:15) cr , β (cid:15) is small, and the range, where D ( ω m ) oscillates,is confined to small ω m ≤ ¯ ge − π/β (cid:15) . By properly taking the double limit (cid:15) → (cid:15) cr and ω m → ω m = 0. At (cid:15) = (cid:15) cr all these gap functions coincide with ∆ ( ω m )at any ω m >
0. This agrees with the observation in Paper IV that as γ increases towards 2,the region, where ∆ n ( ω m ) with finite n change sign, gets confined to progressively smaller ω m , while at larger ω m , all ∆ n ( ω m ) with n = 0 , , ... nearly coincide. We illustrate this inFig.3. For consistency with the notations in previous sections, it is convenient to introduce ξ = ( (cid:15) cr − (cid:15) ) /(cid:15) and label the continuum set of the gap functions by ∆ ξ ( ω m ). Then the endpoint solutions (cid:15) → (cid:15) = (cid:15) cr are ∆ ∞ ( ω m ) and ∆ ( ω m ).It is beyond the ability of the order-by-order expansion to determine the form of ∆ ξ ( ω m )near ξ = 0. On general grounds, we expect that corrections to f (cid:15) ( ω m ) → f ξ ( ω m ) in (31)become relevant starting already from small frequencies, and that at ξ = 0, the gap functioncoincides with ∆ ( ω m ), which we found in the previous section. A way to reproduce thisbehavior is to assume that at ξ →
0, the series for D ( ω m ) in (35) become geometrical( Q n +1 ,(cid:15) cr (cid:15) n +1 max ≈ ( − n ). In this case, D ξ ( ω m ) ∼ cos f ξ ( ω m ) α + cos f ξ ( ω m ) sign ω m , (36)where α ∼ ξ and f ξ ( ω m ) ∼ √ α log ¯ g/ | ω m | + f ∗ ( ω m ), where f ∗ ( ω m ) is a regular function of ω m , which at low frequencies reduces to π/ O ( ω m / ¯ g ). For any ξ >
0, this D ξ ( ω m ) changessign an infinite number of times, but at ξ = 0, D ξ =0 ( ω m ) ∼ ¯ g/ω m , as we expect. We alsonote that between the nodes, D ξ ( ω m ) from (36) is large, of order 1 /ξ . Extending this D ( ω m )to complex frequencies, z = ω (cid:48) + iω (cid:48)(cid:48) , we find that there exist anti-vortices at small z in thelower frequency half-plane. At ξ = 0, vortices and anti-vortices annihilate at z = 0, leavinga regular gap function ∆ ( ω m ).In Appendix C we consider the extended γ − model with non-equal interactions in theparticle-particle and particle-hole channels and use M (cid:54) = 1 as a measure of the differenceof the two interactions. For the extended model, there is a critical M max , below which theground state is a non-Fermi liquid with ∆ = 0. For γ = 2, M max = 0. We obtain the set15 ( ω m ) Δ ω m Δ ω m Δ ω m C = C cr Δ ω m C < C cr Δ ( ω m )Δ ( ω m ) Δ ∞ ( ω m ) γ < 2 γ = 2 ϵ = ϵ max ϵ < ϵ max ξ = 0 ξ > 0 FIG. 3. The gap function ∆ n ( ω m ) for γ < γ = 2. For γ <
2, ∆ n ( ω m ) changes sign n times. As γ gets close to 2, the frequency region where these n sign changes happen, shrinks toprogressively smaller ω m = 0, and at γ = 2 −
0, ∆ n ( ω m ) with finite n collapse into ∆ ( ω m ) at all ω m >
0. The continuum set of ∆ ξ ( ω m ) at γ = 2 − n ( ω m ) with n → ∞ , and thecontinuous parameter ξ is determined by how the double limit n → ∞ and γ → ξ ( ω m ) with ξ > ω m = 0 and ω m ∼ ¯ g . The solution of the linearized gap equation is the ξ → ∞ limit of this set. of ∆ ξ ( ω m ) at small ω m at M = 0+ and show that all gap functions from the continuous setappear simultaneously with the overall magnitude M / .We next analyze the condensation energy E c . We define E c as the difference between theactual ground state energy E ∆ at a finite ∆( ω m ) and the would be ground state energy of thenormal state, E ∆=0 . The expression for E c for γ = 2 has been obtained before: [35, 37–39]and we just copy it here: E c = − N (cid:90) ∞ dω m ω m ( (cid:112) D ( ω m ) − (cid:112) D ( ω m ) − N ¯ g (cid:90) ∞ dω m dω (cid:48) m (cid:16)(cid:112) D ( ω m ) − (cid:112) D ( ω (cid:48) m ) (cid:17) (cid:112) D ( ω m ) (cid:112) D ( ω (cid:48) m ) ω m ω (cid:48) m ( ω m − ( ω (cid:48) m ) ) (37)16 E c , n . . . n = 0 n = 1 n = 2 γ < 2 γ = 2 ξE c , ξ , all solutions with fi nite ξ = 0 n correspondingly ξ > 0 n → ∞ (a) (b) , n = ∞ FIG. 4. (a) The condensation energy E c the solutions of the Eliashberg gap equation for γ < E c = E c,n is a discrete function of a number of a solution, n . The largest condensation energyis for n = 0. (b) The condensation energy E c,ξ for γ = 2. E c,ξ is a continuous function of theparameter ξ . The condensation energy at ξ = 0 is the accumulation point of all E c,n from γ < n = 0 , .. . Every other point on the curve E c,ξ comes from the limit n → ∞ , anddifferent ξ correspond to different ways how the double limit n → ∞ and γ → ξ → ∞ , E c is the condensation energy for infinitesimally small gap function ∆ ∞ ( ω m ). This formula has been derived with the use of (2) and is therefore valid only for the solutionsof the gap equation. Both terms in (37) are negative, i.e. any solution of the gap equationlowers the ground state energy compared to the normal state.Substituting (35) into (37), we find that E c = E c,ξ is a continuous function of ξ . At ξ (cid:29) E c,ξ = − aN ¯ g ξ (38)where a = O (1). It is natural to expect that | E c,ξ | increases with decreasing ξ and reachesa maximum at ξ = 0, see Fig. 4. [40]In the next two Sections we present corroborative evidence for the special, critical behav-ior of the γ -model with γ = 2 from the analysis of the gap function on the real frequencyaxis and in the upper half-plane of frequency.17 V. GAP EQUATION ALONG THE REAL FREQUENCY AXIS
As we said in the Introduction, the analysis of the gap equation for the γ modelalong real frequency axis should generally be more revealing than the analysis along theMatsubara axis, because the pairing interaction on the real axis V (Ω) = (cos( πγ/
2) + i sign(Ω) sin( πγ/ g/ | Ω | ) γ is complex. The real part of the interaction becomes repulsivefor γ >
1, and the imaginary part vanishes at γ = 2 for any non-zero Ω. This makes the γ = 2 case special.We present the results for ∆( ω ) on the real axis in the same order as in previous section:we first obtain the solution of the linearized gap equation, which we label ∆ ∞ ( ω ), thenanalyze the solution ∆ ( ω ), and then show that there is a one-parameter continuous set ofsolutions ∆ (cid:15) ( ω ) in between ∆ ∞ ( ω ) and ∆ ( ω ). A. Linearized gap equation in real frequencies
The linearized gap equation in real frequencies is obtained by taking the limit ∆( ω ) → D ∞ ( ω ) = ∆ ∞ ( ω ) /ω and re-write the gap equation as D ∞ ( ω ) = − ¯ g ω (cid:20) i π dD ∞ ( ω ) dω sign ω + D ∞ ( ω ) ω + (cid:90) ∞ dω (cid:48) ( | ω | + ω (cid:48) ) (cid:60) D ∞ ( ω (cid:48) ) (cid:21) , (39)where (cid:60) stands for the real part. The D ∞ ( ω ) term in the l.h.s. of (39) is the analog of D ∞ ( ω m ) in the l.h.s. of the gap equation (9) on the Matsubara axis, and, like there, itoriginates from the bare ω term in the fermionic Green’s function. Neglecting this term, wefind that the solution of (39) is D ∞ ( ω ) = − i(cid:15) cos (cid:34) β (cid:32) log (cid:18) ω ¯ g (cid:19) − iπ sign( ω ) (cid:33) + φ (cid:35) , (40)where β = 0 . (cid:15) is infinitesimally small. We note that this D ∞ ( ω ) can be obtained from D ∞ ( ω m ), Eq. (11), by rotating from iω m to ω + i
0. In explicitform, D (cid:48)∞ ( ω ) = 2 (cid:15) sin (cid:32) β log (cid:18) ω ¯ g (cid:19) + φ (cid:33) sinh( πβ )sign ωD (cid:48)(cid:48) ∞ ( ω ) = − (cid:15) cos (cid:32) β log (cid:18) ω ¯ g (cid:19) + φ (cid:33) cosh( πβ ) (41)18bserve that D (cid:48)∞ ( − ω ) = − D (cid:48)∞ ( ω ) and D (cid:48)(cid:48) ∞ ( − ω ) = D (cid:48)(cid:48) ∞ ( ω ), as it should be. The relation (cid:90) ∞ dx x iβ ( x + 1) = 2 πβ sinh(2 πβ ) = 1sinh ( πβ ) , (42)is useful for the verification that D ∞ ( ω ) satisfies Eq. (39) without D ∞ ( ω ) in the l.h.s. Usinganother relation (cid:90) ∞ dx x iβ x − iπ coth( πβ ) , (43)one can verify that D (cid:48)∞ and D (cid:48)(cid:48) ∞ satisfy KK relations:2 π (cid:90) dx D (cid:48)∞ ( x ) xx − ω = − D (cid:48)(cid:48) ∞ ( ω ) , ωπ (cid:90) dx D (cid:48)(cid:48) ∞ ( x ) x − ω = D (cid:48) ∞ ( ω ) , (44)where the integrals are principle values.We next consider | ω | > ¯ g . To obtain D ∞ ( ω ) in this region, we take as an input the exactsolution on the Matsubara axis and analytically continue it to the real axis. By construction,this can be done by replacing ω m by − iz - the function ∆ ∞ ( z ) is guaranteed to be analyticin the upper half-plane of frequency. However, because we don’t have the exact analyticalexpression for D ∞ ( ω m ) for an arbitrary ω m , have to replace ω m by − i ( ω + i
0) in Eq. (13) andobtain ∆ ∞ ( ω ) by integrating over k . For small ω < ¯ g , we find, after this integration, seriesof corrections to (40) in ω/ ¯ g . For large ω > ¯ g , the largest contribution to ∆ ∞ ( ω ) comesfrom the continuation of the universal oscillating term ∆ ∞ ; u , Eq. (17). Upon rotation tothe real axis, this term splits into two. One remains exponentially small, but in the otherthe exponential factor cancels out. As a result, on the real axis we have (see Appendix Bfor details) D ∞ ; u ( ω ) ∼ √ (cid:15)e iπ (cid:104) ( ω ¯ g ) +log ( ω ¯ g ) (cid:105) . (45)Other contributions contain powers of ¯ g/ | ω | and are smaller. Neglecting them, we obtain D ∞ ( ω ) = D ∞ ; u ( ω ) at ω (cid:29) ¯ g .Comparing this form with (40), we see that both D (cid:48)∞ ( ω ) and D (cid:48)(cid:48) ∞ ( ω ) continue oscillatingat ω > ¯ g , but with the period set predominantly by ( ω/ ¯ g ) rather than by log ( ω/ ¯ g ) . InFig. 5, we plot real and imaginary parts of D ∞ ( ω ) and the phase of the gap, η ∞ ( ω ), definedvia D ∞ ( ω ) = | D ∞ ( ω ) | e iη ∞ ( ω ) , or, equivalently, via η ∞ ( ω ) = Im log D ∞ ( ω ) = Im log ∆ ∞ ( ω ).We see that the phase winds up an infinite number of times between ω = 0 and O (¯ g ) andbetween O (¯ g ) and ∞ . Oscillations at ω < ¯ g are directly related to oscillations of ∆ ∞ ( ω m )on the Matsubara axis, and there is one-to-one correspondence between each phase winding19 b)(a) FIG. 5. Real and imaginary part (a) and the phase η ∞ ( ω ) (b) of D ∞ ( ω ). The periodicity ofoscillation is set by [( ω/ ¯ g ) + log( ω/ ¯ g ) ] /π . by 2 π on a real axis and a vortex on the Matsubara axis. Oscillations and phase windingat ω > ¯ g are present on the real axis, but not on the Matsubara axis. It is natural to relatethis discrepancy to the fact that the pairing interaction is attractive on the Matsubara axis,but on the real axis, V (cid:48) (Ω) is repulsive, and a non-zero D ∞ ( ω ) comes from V (cid:48)(cid:48) (Ω) ∝ δ (Ω )(see more on this below). B. The function D ( ω ) We now consider the opposite limit of the real-axis partner of sign-preserving D ( ω m ).At ω m < ¯ g , D ( ω m ) ≈ ∆ (0) /ω m , and D ( ω ) on the real axis must also be close to D (0) /ω ω < ¯ g . At larger ω m > ¯ g , we will see that D ( ω m ) and D ( ω ) are very different: D ( ω m )decays as 1 /ω m , while D ( ω ) does not decay and oscillates in sign.The solution of the gap equation along the real axis for ω > ¯ g has been found byCombescot [9], who build his analysis on earlier results by Karakozov, Maksimov, andMikhailovsky [7] and by Marsiglio and Carbotte [8]. We follow Ref. [9] below.It is convenient to introduce φ ( ω ) via D ( ω ) = 1 / sin φ ( ω ) and re-express the gapequation (7) at T = 0 as the equation on φ ( ω ). The equation is dφ ( ω ) dω = 2 π ¯ g [ ωB ( ω ) − A ( ω ) sin φ ] (46)where A ( ω ) and B ( ω ) are given by Eq. (6). The initial condition for φ is φ (¯ g ) ≈ ¯ g/ ∆ (0) = O (1), consistent with φ ( ω ) = ω/ ∆ (0) at ω < ¯ g .At ω ≥ ¯ g , B ( ω ) ≈ g /ω and A ( ω ) ≈ − α ¯ g /ω , where α ≈ .
27 (Ref. [9]) . The A ( ω )term can then be neglected if φ ( ω ) is real, as we will assume and then verify. Without thisterm, Eq. (46) can be solved easily, and the result is φ ( ω ) ≈ π (cid:32) log (cid:18) ω ¯ g (cid:19) + (cid:18) ω ¯ g (cid:19) + C (cid:33) (47)where C = ¯ g/ ∆ (0) − /π . We see that φ ( ω ) is real, as we anticipated. We note that this φ ( ω ) coincides with the argument of the exponent for D ∞ ( ω ) in (45)The function D ( ω ) = 1sin φ ( ω + i
0) (48)is a sign-changing function of ω , whose real part almost diverges at a set of frequencieswhere φ ( ω ) = pπ , where p = 1 , .. is an integer. The imaginary component D (cid:48)(cid:48) ( ω ) is a setof δ − functions at these frequencies. We plot the gap function ∆ ( ω ) = ωD ( ω ) in Fig.6.To analyze the phase winding, we again introduce the phase factor via D ( ω ) = | D ( ω ) | e iη ( ω ) and consider how η ( ω ) varies at ω ≥ ¯ g . The imaginary component D (cid:48)(cid:48) ( ω ) in(48) is infinitesimally small, except in the vicinity of ω p , where φ ( ω p ) = pπ . We use Eq.(47) for φ ( ω ) and express D ( ω ) near each such point as D ( ω ) ≈ π ¯ g − p ω p ω p + ¯ g ω − ω p + iδ . (49)Then e iη ( ω ) = ( − p ω − ω p − iδ (cid:112) ( ω − ω p ) + δ (50)21 ω / g - - Δ ( ω ) Δ ' ( ω ) Δ '' ( ω ) FIG. 6. ∆ ( ω ) = ω/ sin φ ( ω ) for φ ( ω ) given by (47). The real part of the gap ∆ (cid:48) ( ω ) diverges atthe set of points where φ ( ω ) = pπ , p = 1 , ... . The imaginary part ∆ (cid:48)(cid:48) ( ω ) is a set of δ − functionsat these points. The behavior of ∆ (cid:48) ( ω ) has been obtained in Refs. [7–9, 11]. ω / g - π - π π πη ( ω ) FIG. 7. Variation of the phase of the gap η ( ω ) (∆ ( ω ) = | ∆ ( ω ) | e iη ( ω ) ). We restrict η ( ω ) to( − π, π ). Phase slips of η ( ω ) continue up to infinite frequency. We plot η ( ω ) in Fig.7. We see that the phase rapidly changes by π around each ω p . Ifwe restrict η ( ω ) to ( − π, π ), we find that the phase jumps by 2 π in between ω p and ω p +1 .The number of ω p is infinite, hence the number 2 π jumps is also infinite. We reiterate thatbehavior has no analog the Matsubara axis, where D ( iω m ) is real and positive for ω m > η = 0. [41] 22 . The one-parameter set of gap functions We follow the same strategy as in the analysis on the Matsubara axis and expand thenon-linear gap equation (7) in powers of D . We search for the solution in the form D ( ω ) = ∞ (cid:88) j =0 D (2 j +1) ( ω ) (51)where D (1) ( ω ) = D ∞ ( ω ) and higher-order terms are obtained by solving Eq. (7) iteratively.For ω < ¯ g , we use Eq. (40) for D ∞ ( ω ). The computational steps are the same as in Sec.III C, and we obtain D (cid:15) ( ω ) = − i(cid:15) (cid:16) cos ˜ f (cid:15) ( ω ) + (cid:15) Q ,(cid:15) cos 3 ˜ f (cid:15) ( ω ) + (cid:15) Q ,(cid:15) cos 5 ˜ f (cid:15) ( ω ) + ... (cid:17) (52)where Q i,(cid:15) are the same as in (35) and˜ f (cid:15) ( ω ) = β (cid:15) (cid:32) log (cid:18) ω ¯ g (cid:19) − iπ sign ω (cid:33) + φ (cid:15) (53)This D (cid:15) ( ω ) could also be obtained directly from (35) by replacing log ω m by log ω − iπ ineach term in (35).We recall that the continuous set exists for (cid:15) ≤ (cid:15) cr . For any (cid:15) < (cid:15) cr , D ( ω ) oscillatesan infinite number of times down to ω = 0. As (cid:15) approaches (cid:15) cr , log-oscillations shift toprogressively smaller frequencies. At (cid:15) = (cid:15) cr , β (cid:15) vanishes and log-oscillations disappear. Thebehavior of D ( ω ) at ω → (cid:15) → (cid:15) cr depends on how the double limit ω → (cid:15) → (cid:15) cr is taken.Like we did in Sec. III C, we introduce ξ = ( (cid:15) cr − (cid:15) ) /(cid:15) and re-express ∆ (cid:15) ( ω ) as ∆ ξ ( ω ).The two limits (cid:15) = 0 and (cid:15) = (cid:15) cr now correspond to ξ = ∞ and ξ = 0, respectively. Thisbrings the notations in line with the ones we used in Secs. IV A and IV B.On the Matsubara frequency, all ∆ ξ ( ω m ) behave in the same way at ω m > ¯ g : ∆ ξ ∝ /ω m .On the real axis, the dependence on ξ is more complex. To see this, we use the solution ofthe linearized gap equation D (1) ∝ ie iφ ( ω ) with φ ( ω ), given by (47), and evaluate D (2 n +1) in order-by-order expansion of the non-liner gap equation in D . Collecting the series, weobtain the closed form expression D ξ ( ω ) = − ie iφ ( ω ) ξ − e iφ ( ω ) / (1 + ξ ) = 1sin[ φ ( ω ) + i log(1 + ξ )] (54)23his expression can be equivalently obtained by solving Eq. (46) for φ ( ω ) with the initialcondition φ (¯ g ) = ¯ g/ ∆ (0) + i log (1 + ξ ).The parameter ξ runs between 0 and ∞ . For ξ = 0, Eq. (54) yields D ( ω ) = 1 / sin φ ( ω ),which agrees with (48) (one should add i ω in this case). At ξ → ∞ we recover, byconstruction, the solution of the linearized gap equation. For any ξ , including ξ = 0, D ξ ( ω )oscillates up to an infinite frequency, and its phase η ξ ( ω ) winds up by an infinite number of2 π between ω ∼ ¯ g and ω = ∞ .We see therefore that in both limits ω (cid:28) ¯ g and ω (cid:29) ¯ g , the solutions of the non-lineargap equation form a continuous one-parameter set, Eqs. (52) and (54). We conjecture thatfor any value of ξ , one can use a free phase factor φ ξ in (53) to merge small- ω and large- ω expressions into a single D ξ ( ω ). D. density of states
The fermionic DoS is defined as N ( ω ) = ( − N /π )Im G l ( ω ), where N is the DoS in thenormal state and G l ( ω ) = − iπ (cid:115) − D ( ω ) (55)is a retarded Green’s function, integrated over the dispersion.In a BCS superconductor, N ( ω ) ∝ Re ω/ √ ω − ∆ vanishes at ω < ∆, and has anintegrable singularity at ω = ∆ + 0, and is non-zero for all ω > ∆ because quasiparticlestates in a BCS superconductor form a continuum ω = ± (cid:112) ∆ + ( (cid:15) k − µ ) . In our case, theform of N ( ω ) = N ξ ( ω ) strongly depends on ξ . At small ω < ¯ g , N ξ ( ω ) remains finite downto ω = 0 for all ξ >
0. In this respect, all such solutions describe gapless superconductivity.The gap function ∆ ( ω ) tends to a finite ∆ (0) at small ω , and the corresponding N ( ω )vanishes, like in BCS superconductor. We show this in Fig.8(a).At ω > ¯ g , ∆ ξ ( ω ) is given by (54), and N ξ ( ω ) = N (cid:61) tan[ φ ( ω ) + i log (1 + ξ )]. For ξ > N ξ ( ω ) oscillates around N up to ω = ∞ . The amplitude of the oscillations increases withdecreasing ξ . For ξ = 0, N ( ω ) = N δ/ (cos φ ( ω ) + δ ), where δ = 0+. This DoS consistsof a set of δ − functions at the points ω k , for which φ ( ω k ) = π/ kπ ( k is an integer)(Refs [9, 11, 12]). We show this in Fig.8(b) In Fig.8(c) we show N ξ ( ω ) in the whole rangeof frequencies. 24 - - ω / g N ( ω ) / N ϵ = ϵ = ϵ = ϵ = ϵ max - - ω / g N ( ω ) / N ϵ = ϵ = ϵ = ϵ = ϵ max ω / g N ( ω ) / N ϵ * = ϵ * = ϵ * = ϵ * = ω / g N ( ω ) / N ϵ * = ϵ * = ϵ * = ϵ * = ω / g N ( ω ) / N ϵ * = ϵ * = ϵ * = ϵ * = - - ω / g N ( ω ) / N ϵ = ϵ = ϵ max ϵ * = ϵ * = - - ω / g N ( ω ) / N ϵ = ϵ = ϵ max ϵ , ϵ * → ϵ * → ϵ→ϵ max (a) (b) (c) ξ = 4 ξ = 5.7 ξ = 11.5 ξ = 0.1 ξ = 9 ξ = 2.3 ξ = 0.43 ξ = 0.1 ξ → ∞ ξ → 0 FIG. 8. DoS N ξ ( ω ) for (a) ω < ¯ g and (b) ω > ¯ g and for different ξ . For all ξ > N ξ ( ω ) remainsfinite down to ω = 0 (a gapless superconductivity). For ξ = 0, the DoS N ( ω ) vanishes at small ω and has a set of δ − functional peaks at ω > ¯ g . In panel (c) we present the schematic plot of theDoS at all frequencies. The function N ( ω ) is the true DoS at T = 0, as the ξ = 0 solution has the lowestcondensation energy. It is different from the DoS in a BCS-type superconductor, which isnon-zero at all ω > ∆ and approaches N at ω → ∞ . We emphasize that a qualitativedistinction holds only for γ = 2. For smaller γ , the DoS for the n = 0 solution evolves as afunction of frequency, but still remains non-zero at all ω > ∆ and approaches N at infinite ω (see Paper IV).In a BCS superconductor, a continuous N ( ω ) at ω > ∆ is the consequence of the factthat fermionic energy E k is a continuous function of the normal state dispersion (cid:15) k , E k = (cid:112) (cid:15) k + ∆ . The form of N ( ω ) as a set of δ -functional peaks raises the issue whether fermionicenergies get quantized at γ = 2. To address this issue, we compute the total weight of eachlevel: N k = (1 / π ) (cid:82) N ( ω ) /N , where the integration is confined to the vicinity of ω k . Using φ ( ω ) ≈ ω /π , we obtain N k = 1 / (cid:112) k ). We see that N k < k . Because ofthis, ω k cannot be viewed as true quantized fermionic energy levels, as a fermion is necessarydistributed between ω k with different k . 25 . GAP FUNCTION IN THE UPPER FREQUENCY HALF-PLANE. Comparing D ξ ( ω ) and D ξ ( ω m ), we see that they are similar at small frequencies, but verydifferent at ω, ω m > ¯ g . Indeed, on the real axis, the phase η ξ ( ω ) winds up by an infinitenumber of 2 π between ω = O (¯ g ) and ω = ∞ , while near the Matsubara axis, η ξ ( ω m ) = 0in this frequency range. The discrepancy implies that phase winding must end somewherebetween the real and the Matsubara axis. We now argue that there is a set of vortices inthe upper frequency half-plane, at | z | ≥ ¯ g and the phase winding drops by 2 π each time theaxes of z passes through a vortex upon rotation away from the real axis.We use the Cauchy relation ∆( z ) = 2 π (cid:90) ∞ dx x ∆ (cid:48)(cid:48) ( x ) x − z (56)to extend the gap function ∆( x ) = xD ( x ) from the real axis to complex z = ω (cid:48) + iω (cid:48)(cid:48) with ω (cid:48)(cid:48) >
0. We use Eq. (54) for the gap function as we expect vortices to be present at | z | > ¯ g .To simplify the presentation, we first consider the case ξ = 0, then the case ξ → ∞ , andthen extend the analysis to arbitrary ξ . A. case ξ = 0 Using the expansion near φ ( ω p ) = pπ , Eq. (49), we approximate ∆ (cid:48)(cid:48) ( ω ) as∆ (cid:48)(cid:48) ( ω ) ≈ π ¯ g ∞ (cid:88) p =1 ( − p +1 ω p ω p + ¯ g δ ( ω − ω p ) . (57)Substituting into (56), we obtain∆ ( z ) = π ¯ g ∞ (cid:88) p =1 ( − p +1 ω p ω p + ¯ g ω p − z . (58)Here ω p is a solution of φ ( ω p ) = πp , where φ ( ω ) is given by (47). We verified numericallythat KK relations on the real axis are satisfied, i.e., if we use (57), we reproduce with highaccuracy ∆ (cid:48) ( ω ). On the Matsubara axis, z = iω m , Eq. (58) yields, at ω m (cid:29) ¯ g ,∆ ( ω m ) = a ¯ g ω m (59)where a = π (cid:80) ∞ p =1 (cid:2) ( − p +1 ω p / (( ω p + ¯ g )¯ g ) (cid:3) . Approximating ω p by ¯ gπ √ p , we find a = 2 . .
27, obtained by solving the gap equation on the26 ′ / ¯ g ω ′ ′ / ¯ g ω ′ ′ / ¯ g log | Δ( z ) | η ( z ) - - - - - - - - - - - - - π - π / π / π - π - π / π / π - π - π / π / π - π - π / π / π - π - π / π / π - π - π / π / π log | Δ ( z ) | η ( z ) FIG. 9. ∆ ( z ) in the upper half plane. Upper panel: log | ∆ ( z ) | . Blue spots mark the locationsof dynamical vortices, where | ∆ ( z ) | = 0. Lower panel: the phase of the gap η ( z ) defined via∆ ( z ) = | ∆ ( z ) | e iη ( z ) . The phase slips by 2 π upon crossing a white line in the direction from near-white to dark-blue. The white lines are the locations of points where ∆ (cid:48)(cid:48) ( z ) = 0 and ∆ (cid:48) ( z ) < Matsubara axis (Ref. [9] and Sec. III B). The difference likely comes from subleading termsin φ ( ω ).We plot ∆ ( z ) for a generic z in the upper half-plane in Fig.9 We clearly see that thereis a set of points, where ∆ (cid:48) ( z ) = ∆ (cid:48)(cid:48) ( z ) = 0. These points are the centra of dynamicalvortices with anti-clockwise circulation 2 π . The vortices are located along a particular linein the complex plane. The set extends to an infinite frequency, i.e., the number of vorticesis infinite. This is consistent with an infinite phase winding along the real axis. We verifiedthat if we use a more accurate form of ω p , the positions of the vortices shift a bit, but theirnumber remains infinite.To see how the winding number changes once we rotate from the real to the Matsubaraaxis, we introduce z = | z | e iψ ( ψ = 0 along the positive real semi-axis and π/ ( z ) = | ∆ ( z ) | e iη ( z ) between | z | ∼ ¯ g and | z | → ∞ along the directions in the upper frequency half-plane, specified by ψ . Weshow the results in Fig.10. We see that for any ψ >
0, the phase η ( z ) winds for | z | below acertain value, and then saturates. At larger | z | , both ∆ (cid:48) ( z ) and ∆ (cid:48)(cid:48) ( z ) scale as 1 / | z | with27
10 15 20 | z |/ g - ππη (| z | , ψ ) ψ = | z |/ g - ππη (| z | , ψ ) ψ = | z |/ g - ππη (| z | , ψ ) ψ = | z |/ g - ππη (| z | , ψ ) ψ = | z |/ g - ππη (| z | , ψ ) ψ = | z |/ g - ππη (| z | , ψ ) ψ = | z |/ g - ππη (| z | , ψ ) ψ = | z |/ g - ππη (| z | , ψ ) ψ = | z |/ g - ππη (| z | , ψ ) ψ = | z |/ g - ππη (| z | , ψ ) ψ = | z |/ g - ππη (| z | , ψ ) ψ = | z |/ g - ππη (| z | , ψ ) ψ = | z |/ g - ππη (| z | , ψ ) ψ = | z |/ g - ππη (| z | , ψ ) ψ = FIG. 10. Phase variation η ( | z | , ψ ) along different paths specified by ψ , defined via z = | z | e iψ .Along real axis, ψ = 0; along Matsubara axis, ψ = π/
2. Along the real axis the phase η ( ω ) windsup an infinite number of times, i.e., the winding number (the number of 2 π phase slips) is infinite.For a finite ψ , phase winding ends at some finite | z | , and the winding number becomes finite. no oscillations. Counting the total phase winding δη between | z | = O (¯ g ) and | z | = ∞ , wesee that δη = 2 πs , where s is an integer. It decreases by one every time the direction setby ψ passes through a vortex. The winding vanishes at some ψ ≤ π/ - - FIG. 11. The phase η ∞ ( z ) of ∆ ∞ ( z ) in the first quarter of the complex plane of frequency ( ω (cid:48) > ω (cid:48)(cid:48) > B. case ξ = ∞ We next consider the opposite limit ξ = ∞ . The form of ∆ ∞ ( z ) can be obtained startingfrom (13) and replacing ω m by | z | e i (2 ψ − π ) . This gives∆ ∞ ( z ) ∝ (cid:90) ∞ dk (cid:16) b k e − ik log | z | / ¯ g +(2 ψ − π ) k + b − k e ik log | z | / ¯ g − (2 ψ − π ) k (cid:17) , (60)where b k is defined in (14). We obtain ∆ ∞ ( z ) by numerical integration. We plot its phase η ∞ ( z ) in Fig. 11. We again see that there is an infinite array of vortices. The array extendsto an infinite frequency, where it approaches the real axis. The vortex arrangement in Fig. 11is remarkable similar to that in Fig. 9 for ξ = 0. Moreover, if we approximate φ ( ω ) by theleading term ( ω/ ¯ g ) /π , we find that the positions of the vortices are at the same z i in bothcases. We can see this by comparing Fig.12(a) where ξ → ∞ and Fig.12(c) where ξ → ξ ( z ) are very similar in these two cases, despite that the overall factorsare different. The vortex positions for these two cases are almost identical, as can be seenfrom Fig.12(d). 29 rinted by Wolfram Mathematica Student Edition Printed by Wolfram Mathematica Student Edition
Printed by Wolfram Mathematica Student Edition ω ' ω '' Printed by Wolfram Mathematica Student Edition - - - - - - - - - - - log | Δ ξ | , ξ = 2 × 10 log | Δ ξ | , ξ = 1log | Δ ξ | , ξ = 0.01 ξ = 0.01 ξ = 1 ξ = 2 × 10 (a)(b)(c)(d) ω ′ ω ′ ′ ω ′ ′ ω ′ ′ ω ′ ′ FIG. 12. (a)-(c)Gap functions ∆ ξ defined in (54) for different ξ in the frequency upper half plane.Here we take φ ( ω ) ≈ ( ω/ ¯ g ) /π (d) Comparison of the vortex positions for different ξ obtained byapproximating φ ( ω ) by ( ω/ ¯ g ) /π . The results suggest that the positions of the vortices almost donot depend on the value of ξ . C. arbitrary ξ The same infinite array of vortices exists for all 0 < ξ < ∞ . As an example, in Fig.12(b)we show the gap function for ξ = 1. We clearly see that there is an infinite array of vortices,similar to the ones for ξ = ∞ and ξ = 0. If we again approximate φ ( ω ) by the leading term( ω/ ¯ g ) /π , we find that the positions of the vortices for ξ = 1 almost coincide with those formuch larger and much smaller ξ . We demonstrate this in Fig.12(d). D. Essential singularity
There is another consequence of the existence of an infinite array of vortices – each the gapfunction ∆ ξ ( z ) has an essential singularity at | z | = ∞ . Indeed, one can reach | z | = ∞ fromthe set of vortex points, where ∆ ξ ( z ) = 0, or from the real axis, where ∆ ξ ( ω ) oscillates, andthe amplitude of the oscillations does not vanish at ω → ∞ , hence neither lim | z |→∞ ∆ ξ ( | z | )nor lim | z |→∞ / ∆ ξ ( | z | ) exist. We emphasize that an essential singularity is only present for30 = 2. For smaller γ , phase winding and associated vortices exist only at | z | smaller than acertain, γ − dependent value. At larger | z | , ∆( z ) scales as 1 / | z | γ and vanishes at | z | = ∞ nomatter how this limit is reached.Further, for γ = 2, the very existence of a non-zero ∆ ξ ( z ) for a generic z away from vortexpoints, is ultimately related to an essential singularity at | z | = ∞ . The argument is thatthe set of vortex points is complete, hence one can analytically continue the gap functionfrom this set to the upper half-plane of frequency in the same way as ∆( z ) is obtainedfrom a discrete set of Matsubara points ω m = πT (2 m + 1) in a standard diagrammaticcalculations for interacting fermions. If this analytical continuation was unique, we wouldobtain ∆( z ) ≡
0, because ∆( z ) = 0 at the vortex points. To provide a non-zero ∆ ξ ( z ), theextension must be multi-valued. A rigorous mathematical argument is that this is the casewhen the end point of the set, | z | = ∞ , goes outside the domain of analyticity. This isexactly what we have because of an essential singularity at | z | = ∞ .We conjecture that the multi-value nature of the extension is the reason why the set of∆( ω ) is a continuous one at γ = 2. This is plausible, particularly if the vortices are at thesame z i for all ξ , as Fig.12 seems to indicate. However, at the moment, we cannot provethis. VI. FINITE ω D A. Gap equation at a finite ω D We now consider the case when the bosonic mass is small but finite. By analogy with thephonon case we call this mass ω D . On the Matsubara axis, ∆ ( ω m ) changes little comparedto the case ω D = 0. The set of ∆ n ( ω m ) still exists at small ω D , but becomes discrete andholds up to a finite n max . In particular, there is no solution of the linearized gap equation at T = 0 ( n = ∞ ) for any non-zero ω D . The value of n max can be estimated by noticing that ifwe, e.g., depart from the solution on the Matsubara axis at ω D = 0 and compute correctionsdue to finite ω D , these corrections increase at small ω m and become O (1) at ω m ∼ ω D . Asimple experimentation shows that this sets n max at n max ∼ ¯ gω D (61)31n the real axis, the gap equation still has the form D ( ω ) ωB ( ω ) = A ( ω ) + C ( ω ), and A ( ω )and B ( ω ) remain the same as in (6), up to irrelevant small corrections. However, C ( ω )changes to C ( ω ) = − i π ¯ g ω D D ( ω − ω D ) − D ( ω ) (cid:112) − D ( ω − ω D ) sign ω (62)Expanding to first order in ω D and introducing, as before, D ( ω ) = 1 / sin φ ( ω ), we obtainafter straightforward algebra that the gap equation reduces to˙ φ − ω D (cid:16) ( ˙ φ ) tan φ ( ω ) + ¨ φ (cid:17) = 2 π ¯ g [ ωB ( ω ) − A ( ω ) sin φ ( ω )] + ... (63)where dots stand for the terms with higher powers of ω D . A similar equation at a finite T instead of finite ω D has been obtained by Combescot [9]For definiteness, let’s consider the case ξ = 0. The initial conditions are φ (0+) = 0, dφdω (0+) = 1 / ∆ (0). At ω ≥ ¯ g , B ( ω ) and A ( ω ) from (6) can be approximated by B ( ω ) ≈ g /ω and A ( ω ) ≈ − . g /ω . To understand the effect of ω D we use as an input thesolution at ω D =0, φ ( ω ) ≈ ω / ( π ¯ g ) + iδ . Substituting this input into (63), expanding near ω = π ¯ g/ √
2, where φ ( ω ) = π/
2, and solving for φ (cid:48)(cid:48) ( ω ), we find that it jumps to O ( ω D ) once ω exceeds π ¯ g/ √
2. The same happens at all ω n = π ¯ g/ √ n + 1) / , where tan φ (cid:48) ( ω ) = 0.After n jumps, φ (cid:48)(cid:48) ( ω ) becomes φ (cid:48)(cid:48) ( ω ) = πω D √ g n (cid:88) m =1 √ m + 1 ≈ πω D g n / ≈ ω D π ω ¯ g . (64)A more accurate, non-perturbative analysis of (63) shows that φ (cid:48)(cid:48) ( ω ) appears slightly before φ (cid:48) ( ω ) reaches π/
2. This smoothes up the jumps, but the functional form of φ (cid:48)(cid:48) ( ω ) in (64)remains intact. When both φ (cid:48) and φ (cid:48)(cid:48) , D ( ω ) is a complex function of frequency: D (cid:48) ( ω ) = ω sin φ (cid:48) ( ω ) cosh φ (cid:48)(cid:48) ( ω )sin φ (cid:48) ( ω ) + sinh φ (cid:48)(cid:48) ( ω ) sign ω (65)and D (cid:48)(cid:48) ( ω ) = − ω cos φ (cid:48) ( ω ) sinh φ (cid:48)(cid:48) ( ω )sin φ (cid:48) ( ω ) + sinh φ (cid:48)(cid:48) ( ω ) (66)At ω > π ¯ g / (2 ω D ), φ (cid:48)(cid:48) ( ω ) becomes larger than one. At such frequencies, both D (cid:48) ( ω ) and D (cid:48)(cid:48) ( ω ) oscillate with progressively decreasing magnitudes, approximately as the real and theimaginary parts of − ie iω / ( π ¯ g ) e − ωD π g ( ω/ ¯ g ) , (67)and the phase η ( ω ) gradually winds up as ω increases. We show this in Fig.13. This32 ω / g - π - π π πη ( ω ) - - - ω / g Δ ( ω ) Δ′ ( ω ) Δ′ ′ ( ω ) FIG. 13. The gap function ∆( ω ) (a) and the variation of its phase η ( ω ) at ω D = 0 . g . From (67). behavior holds as long as | A ( ω ) | (cid:28) ω , i.e., ω < ω max , where ω max ∼ ¯ g (cid:18) ¯ gω D log ¯ gω D (cid:19) / (68)At even larger frequencies, the A ( ω ) term cannot be neglected, and the forms of φ (cid:48) ( ω ) and φ (cid:48)(cid:48) ( ω ) change. We show the full numerical solution of Eq. (63) in Fig.14. We see that at ω > ω max , φ (cid:48)(cid:48) ( ω ) keeps increasing, while φ (cid:48) ( ω ) saturates. A simple analysis shows that Eq.(63) is satisfied, up to corrections of order ω D , if φ (cid:48)(cid:48) ( ω ) = 3 log ω ¯ g + 0 . , φ (cid:48) ( ω ) = − π mπ (69)where m is an integer. Substituting this complex φ ( ω ) into ∆ ( ω ) = ω/ sin φ ( ω ) ≈− iωe iφ ( ω ) , we see that at ω > ω max , the real part of the gap function gradually decreasesas ∆ (cid:48) ( ω ) = 1 . g ω (70)To obtain ∆ (cid:48)(cid:48) ( ω ) at these frequencies, we need to keep the ω D term in the l.h.s. of (63) andobtain the correction to (69), which we label as ˜ φ . Solving perturbatively for ˜ φ ( ω ) we obtain˜ φ ( ω ) = f (cid:18) ωω max (cid:19) e iω / ( π ¯ g ) (71)where f ( ... ) is a decreasing function of the argument. The ω oscillations of ˜ φ ( ω ) are clearlyvisible in the numerical results for φ (cid:48) and φ (cid:48)(cid:48) in Fig. 14. Substituting (71) into (66), we33
10 15 20 25 30 35 ω / g ϕ ( ω ) ϕ ' ( ω ) ϕ '' ( ω ) ω / g ϕ ( ω ) ∝ω ∝ω ω / g ϕ ( ω ) ∝ω ∝ω ω / g ϕ ( ω ) ∝ω ∝ω
20 25 30 35 ω / g ϕ ( ω ) π
20 25 30 35 ω / g ϕ ( ω ) ( ω / g )
20 25 30 35 ω / g ϕ ( ω ) π
20 25 30 35 ω / g ϕ ( ω ) ( ω / g ) (a) (b)(c) ω / g ϕ ( ω ) ϕ ' ( ω ) ϕ '' ( ω )
20 25 30 35 ω / g ϕ ( ω ) γ = ( γ + ) log ( ω / g )+ log ( / α )
20 25 30 35 ω / g ϕ ( ω ) γ = ( γ + ) log ( ω / g )+ log ( / α ) FIG. 14. (a) Numerical solution of Eq. (63) for a complex φ ( ω ) at ω D = 0 . g . (b) and (c):Asymptotic forms of φ (cid:48) ( ω ) and φ (cid:48)(cid:48) ( ω ). At ω < ω max , φ (cid:48) ( ω ) increases as ω , while φ (cid:48)(cid:48) ( ω ) firstdisplays a step-like behavior and then increases as ω , with ω D in the prefactor. At ω > ω max , φ (cid:48) ( ω )saturates at(2 m − / π , where m = 19 for our chosen ω D , and φ (cid:48)(cid:48) ( ω ) increases logarithmically. Thecorrections to asymptotic values oscillate with the period set by ω /π . In the numerical solution,we neglected ¨ φ term in (63) compared to ( ˙ φ ) and verified that this is a valid approximation. obtain ∆ (cid:48)(cid:48) ( ω ) ∼ ¯ g ω f (cid:18) ωω max (cid:19) cos ω π ¯ g (72)One can verify that an integer m in (69) determines the number of 2 π variations of η ( ω )on the real axis and, equivalently, the number of vortices at complex z i . The value of m decreases one-by-one as ω D increases and ω max decreases. That m is finite implies that thereis no essential singularity at | z | = ∞ . Indeed, at the largest frequencies ∆( ω ) ∝ /ω .For completeness, we verified that higher-order terms in ω D , which we neglected in thel.h.s. of (63), become important at frequencies ω ∼ ¯ g /ω D , which well exceed ω max and aretherefore irrelevant to our purposes. 34 II. DRESSED SUPERFLUID STIFFNESS
In this section we analyze superfluid stiffness and thermal corrections to a superconduct-ing order parameter. As we discussed in the Introduction, we consider the γ = 2-model asthe double limit ω D → E F → ∞ , such that Migdal-Eliashberg parameter λ E = ¯ g N /ω D remains small ( N ∼ /E F is the DoS per unit volume). Accordingly, in the analysis belowwe keep ω D small, but finite.
1. Bare stiffness
A superfluid stiffness is the ratio of the excess energy E η due to inhomogeneous variationof the phase of a superconducting order parameter ∆( r ) = ∆ e iη ( r ) and ( ∇ η ( r )) : E η = ρ s (cid:82) dr ( ∇ η ( r )) . In the momentum space, E η = ρ s (cid:88) q q η q (73)A way to compute ρ s is to choose η q = δ q,q and extract ρ s as the prefactor for q term inthe particle-particle bubble (the sum of GG and FF terms) at zero frequency and finite q (see Refs. [31, 33, 42]At ω D / ¯ g > T = 0 isa fraction of the Fermi energy, ρ s ( T = 0) = E F / (4 π ) (Refs. [31, 33]). This stiffness is muchlarger than T c [43]. At T > ρ s ( T ) drops and vanishes at T c , but at weak coupling a dropof ρ s occurs only in the immediate vicinity of T c .At small ω D / ¯ g , strong mass renormalization m ∗ /m = 1 + ¯ g /ω D changes the stiffness to ρ s ( T = 0) ∼ E F ω D ∆(0)¯ g ∼ T p λ E (74)where T p ∼ ∆(0) is the onset temperature of the pairing. As long as λ E ≤ ρ s ( T = 0) > T p .We now relate the stiffness to the strength of thermal phase fluctuations of ∆( r ) = ∆ e iη ( r ) .For this, consider the correlator (cid:104) η ( r ) η (0) (cid:105) = (cid:82) D [ η ] η ( r ) η (0) e − ρ s (cid:82) dr ( ∇ η ( r )) /T (cid:82) D [ η ] e − ρ s (cid:82) dr ( ∇ η ( r )) /T (75)We assume that in equilibrium η ( r ) = 0 and expand (cid:104) e iη ( r ) (cid:105) as 1 − (cid:104) η ( r ) (cid:105) /
2. Transforming3575) to the momentum space, we obtain (cid:104) ∆( r ) (cid:105) = ∆(1 − (cid:104) η (cid:105) ), where (cid:104) η (cid:105) = 1 N (cid:88) q (cid:81) q (cid:48) (cid:82) dη q (cid:48) η q e − ρ s q η q (cid:48) /T (cid:81) q (cid:48) (cid:82) dη q (cid:48) e − ρ s q η q (cid:48) /T (76)where N is the number of particles in the system. Evaluating the integrals, we obtain theconventional result [44] (cid:104) η (cid:105) = Tρ s N (cid:88) q q . (77)We assume for simplicity that spatial dimension D is larger than 2, in which case the sumconverges. By order of magnitude we then have (cid:104) η (cid:105) ∼ Tρ s ( T ) (78)As long as ρ s ( T ) > T , fluctuation corrections to the order parameter are small. This doesnot hold in the immediate vicinity of the onset temperature of the pairing, T p , but as longas ρ s (0) (cid:29) T c , the T range, where fluctuations are strong and destroy phase coherence, isquite narrow, i.e., superconducting T c remains close to T p . We see that this holds even when ω D is small and the reduction of ρ s by mass renormalization is strong. A. Dressing of ρ s by soft longitudinal fluctuations We now argue that in our case the expression for (cid:104) η (cid:105) is different due to the presence ofa continuum gapless spectrum of condensation energy, E c ( ξ ), where, we remind, ξ , whichruns between 0 and ∞ and the bottom of the spectrum is at ξ = 0. We will need states nearthe bottom of the continuum, at ξ (cid:28)
1. For such states, we assume E c,ξ = E c, + b N N ¯ g ξ , (79)where b = O (1) and N is the total number of particles. We will also need superfluidstiffness ρ s,ξ for the states near the bottom of the continuum. Evaluating the particle-particle susceptibility for a generic ∆ ξ ( ω m ) and extracting the q term we obtain ρ s,ξ ∼ E F ω D ¯ g (cid:90) dω m D ξ ( ω m )1 + D ξ ( ω m ) (80)For ξ = 0, the integral is determined by ω m ∼ ¯ g , where D ( ω m ) ∼
1. This yields ρ s, ∼ E F ω D / ¯ g ∼ T p /λ E , as in (74). For states with ξ >
0, the magnitude of ∆ ξ ( ω m ) is36educed, and the stiffness gets smaller. We assume that for the states near the bottom ofthe continuum, the stiffness is obtained by expanding to first order in ξ : ρ s,ξ = ρ s, (1 − b ξ ) (81)where b = O (1) is positive. The extra energy of a given state ξ due to phase variation is E η,ξ = ρ s,ξ (cid:88) q q η q (82)We assume (see reasoning below) that all states near the bottom of a continuum con-tribute to the variation of the phase, i.e., the averaging in (cid:104) ( η q ) (cid:105) is over both η q and ξ withthe weight factor e − E ξ /T , where E ξ = E c,ξ + E η,ξ = E + δE ξ , (83)and E = E c, + ρ s, (cid:88) q q η q δE (cid:15) = b N g ξ − b ρ s, ξ (cid:88) q q η q (84)If we neglected δE ξ , we would obtain the same result as before: (cid:104) η (cid:105) = Tρ s, N (cid:88) q q (85)Keeping δE ξ we find that (cid:104) η (cid:105) has an additional overall factor, which we label as I T . Drop-ping for simplicity numerical prefactors b and b , we obtain after integrating over η q (cid:48) I T = (cid:82) dξe − Nf ( ξ ) − ξ (cid:82) dξe − Nf ( ξ ) (86)where f ( ξ ) = N ¯ g T ξ − ξ ξ is non-singular. The linear in ξ term in f ( ξ ) comes from integration over η q (cid:48) with q (cid:48) (cid:54) = q (see Eq. 76). Each integration over q (cid:48) yields1 / √ − ξ , and the product of the integrals over all q (cid:48) yields 1 / (1 − ξ ) N/ = e − ( N/
2) log (1 − ξ ) ≈ e ( N/ ξ . 37t small T , the function f ( ξ ) in (87) has a minimum at ξ = T / (4 N ¯ g ) ∼ T / (4 ω D λ E ).Then I T = 1 / (1 − T / (4 ω D λ E )) and (cid:104) η (cid:105) ∼ T λ E T p − T ω D λ E (88)We see that the renormalizations coming from the low-energy states of the continuum spec-trum of the condensation energy hold in powers of T /ω D . With these renormalizations, thefully dressed stiffness is ρ s ( T ) = T p λ E (cid:18) − T ω D λ E (cid:19) (89)We see from that the value of ρ s ( T ) at T → ω D → T = 0, ρ s (0) = T p /λ E is finite and exceeds T p . At ω D → T , and ρ s ( T ) becomes comparable to T at T ∼ ω D λ E . For the largest λ E ∼
1, at which our theory is valid, this holds at T ∼ ω D . It is tempting to associate thistemperature with the actual T c above which the system looses long-range phase coherence.Further, there is an analogy between finite ω D and finite 2 − γ , as the two have similareffect on the gap function (see Paper IV). Replacing ω D by ¯ g (2 − γ ), we find that at ω D = 0and γ <
2, superconducting T c ∼ ¯ g (2 − γ ).Before concluding this Section, we elaborate on our assumption that the averaging overphase fluctuations should include low-energy states from the continuum spectra of the con-densation energy. Consider the case γ ≤
2, when the spectrum is still discrete and the n = 0solution has the lowest condensation energy E c, . The energies E c,n are close to E c, , yetthe solutions with different n are topologically distinct as ∆ n ( ω m ) has n vortices. Theseother states contribute to the renormalization of the phase of ∆ ( ω m ) only if the tunnel-ing amplitude between the states n = 0 and n > E c, and E c,n to be small. The height of the barrier depends on the path alongwhich a state without a vortex transforms into a state with a vortex at some small ω m . Avortex can either come from ω m = ∞ , in which case the barrier is high, or via a creationof a vortex-antivortex pair at ω m = 0, in which case it is low. For a generic γ <
2, ∆ n ( z )are regular at small z in the complex plane, hence one should not expect an anti-vortexnearby. However, for γ = 2, our candidate ∆ ξ ( z ), Eq. (36), possess anti-vortices at small z in the lower frequency half-plane. In this situation, it is natural to expect that for γ ≤ E c, and E c,n with n > III. PHASE DIAGRAM OF THE γ MODELA. γ = 2 , finite ω D In Fig. 15 we the phase diagram for γ = 2 for non-zero ω D and T . At ω D = 0, thetrue transition temperature into a SC state is zero, although the onset temperature for thepairing, T p is finite. At finite ω D , T c is finite but much smaller than T p , at least for small ω D .In between T c and T p the system displays pseudogap behavior: the spectral function andthe DoS display a peak at a finite frequency, but the spectral weight below the peak remainsfinite. Close to T p , the pairing is mainly induced by fermions with the two lowest Matsubarafrequencies ± πT . In this situation, the position of the peak in the spectral function and theDoS increases linearly with T , and the gap fills in as T approaches T p . In the T range near T c , fermions with all Matsubara frequencies contribute to the pairing, and the positions ofthe maxima in the spectral function and the DoS move to smaller frequencies as T increases(gap closing behavior). We show the DoS in the two regimes in Fig.18 below. B. ω D = 0 , < γ ≤ In Fig. 16 we show the phase diagram for the γ model with 0 < γ ≤
2, at ω D = 0 andfinite T . This phase diagram is based on the results of this work and previous works (PapersI-IV). For γ <
2, we found earlier the largest condensation energy is for sign-preservingsolution of the gap equation ( n = 0 in our classification). Still, for any γ >
0, there existsan infinite set of topologically distinct solutions for the gap (all with the same symmetry),labeled by integer n . This generates a discrete spectrum of the condensation energy E c,n .The spectrum is sparse near the bottom at small γ , but becomes dense and flattens up at thebottom as γ approaches 2. At γ ≤
2, the corrections to superconducting order parameterfrom the states with n (cid:54) = 0 are small at low T , but rapidly increase with increasing T anddestroy phase coherence at T c ∼ ¯ g (2 − γ ) /λ E . For γ ≤ T c (cid:28) T p , and there exists awide intermediate temperature range where the system displays a pseudogap behavior. Bycontinuity, we expect that the pseudogap region to exist for all γ > / ¯ g ω D / ¯ g T p T c Normal StatePseudogap SC γ = 2 FIG. 15. The phase diagram of the γ model for γ = 2 in variables ( T / ¯ g, ω D / ¯ g ), where ω D isthe mass of a pairing boson. T p is the onset temperature of the pairing, and T c is the actualsuperconducting transition temperature, below which the system establishes phase coherence. Inbetween the system displays pseudogap behavior, in which fermionic pairs are formed, but there isno macroscopic phase coherence. The dashed line separates the two regimes within the pseudogapphase – the one at higher T , where the system behavior is chiefly determined by fermions with thetwo lowest Matsubara frequencies ± πT , and the one at lower T , when fermions with all Matsubarafrequencies contribute to the pairing. In these two regimes the system displays gap filling and gapclosing behavior, respectively. C. Properties of the pseudogap phase
1. toy model for γ = 2 Let’s start with γ = 2. At T = 0 the DoS is the set of δ − functions, (Fig.17(a)) At afinite T , two new features appear. First, ∆ ( ω ) decreases with increasing ω and displaysno oscillations above ω max , similar to the case with finite ω D discussed in VI. As a result, δ -functional peaks in the DoS at larger frequencies get broadened and eventually disappear.Second, other D ξ ( ω ) from the continuum spectrum of condensation energies contribute to theDoS with Boltzmann factors. For all these solutions, Im∆ ξ ( ω ) remains finite down to ω = 0.As a result, the DoS also becomes non-zero at the smallest ω (this phenomenon is often40 ormal State P s e u dog a p SC T p T c T / ¯ g γ FIG. 16. The phase diagram of the γ model for a generic γ < T and vanishing ω D . For any γ <
2, the true SC transition temperature T c is finite, but is smaller than the onsettemperature for the paring, T p . In between T c and T p , the system displays a pseudogap behavior.There are two distinct behaviors in the pseudogap regime, like in Fig. 15: close to T p , the spectralfunction and the DoS display gap filling behavior, while close to T c , the behavior becomes moreconventional and the gap frequency shifts to a smaller value as T increases. called a gapless superconductivity [34, 35, 45, 46]) We model both effects by introducing aphenomenological ∆( ω ) = ω/ sin( ia + ( ω/ ¯ g ) (1 + ib )), where a and b increase with T . Weshow the corresponding DoS in Fig. 17(b) and (c).
2. gap filling vs gap closing
We argue, based on earlier works[47, 48]), that there are two different regimes of systembehavior within the pseudogap phase. At low T , the position of the peak in the DoS scaleswith ∆ (0) and decreases as T increases (the gap “closes” with increasing T ). At higher T , the peak in the DOS shifts to higher frequencies and the spectral weight below the peakincreases (gap “fills in” with increasing T ). We illustrate this in Fig.18. This last behavioris at least partly related to the fact that in some finite range of T below T p , the gap functionon the Matsubara axis is strongly peaked at the first Matsubara frequencies ± πT (Refs.41 ω / g D o S ω / g D o S ω / g D o S N N N (a) (b) (c) FIG. 17. The density of states, N ( ω ), at different temperatures, for a toy model with ∆ ( ω ) = ω/ sin( ia + ω (1+ ib )), where a and b are two parameters, which increase with T . (a) The T = 0 limit, a = b = 10 − . The DoS has a set of δ -functional peaks. (b) A finite but small T , a = b = 0 .
05. Thefirst few peaks are well defined, but the peaks at large frequencies get overdamped and disappear.(c) A higher temperature, a = b = 0 .
25. The peak at the smallest frequency is still present, atabout the same frequency as at T = 0, but other peaks are washed out, and the spectral weightbelow the peak increases, i.e., the DoS at low frequencies fulls in. [47, 48]). On the real axis the corresponding ∆( ω ) displays ω/T scaling. For such ∆( ω ), thepeak frequency in the DoS increases linearly with T .At a finite ω D and/or 2 − γ , the “gap filling” behavior holds in some range between theonset temperature of the pairing T p and a finite superconducting T c (Fig.18). To estimatethe crossover temperature between the two regimes, we then compare the actual T p withthe one obtained by neglecting the contributions from fermions with ω m = ± πT . We showthe results in Fig. 19. We see that for γ = 2 the onset temperature without ± πT fermionsis strongly reduced – it is about 1 / T p ∼ . g . This implies that the “gapfilling” behavior holds in a wide range below T p and crosses over to “gap closing” behavioronly near T c . IX. CONCLUSIONS
In this paper, we extended our earlier analysis of the γ − model to γ = 2. The γ = 2model describes, among other cases, the pairing, mediated by an Einstein boson, in the limitwhen the bosonic mass ω D tends to zero. On a real axis, the effective interaction in this42 = 0 N ω = 0 N gap closes gap fills in T c T p T = 0 SC Pseudogap Normal State γ < 2 : γ = 2 : T = 0 T p Pseudogap Normal State T low T high T low T high N ω = 0 N ω = 0 T low T high T low T high FIG. 18. The temperature evolution of the DoS N ( ω ). For γ < T < T c . In this regime and in the pseudogap state at T ≥ T c , the temperature variation of N ( ω )resembles that in a conventional BCS superconductor, i.e. when T increases, the position of themaximum in N ( ω ) moves to a smaller frequency. At larger T within the pseudogap phase, N ( ω )displays gap filling behavior when the peak position increases with increasing T and N ( ω = 0)increases towards its normal state value. For γ = 2 (lower panel), T c = 0, but the two differentregimes of pseudogap behavior are present. limit V (Ω) = − ¯ g / Ω is repulsive, and, at a first glance, should not give rise to pairing.However, the same interaction on the Matsubara axis, V (Ω m ) = ¯ g / Ω m , is attractive, andearlier calculations on the Matsubara axis found that the onset temperature of the pairing, T p , tends to a finite value T p = 0 . g at ω D → T p = 0 . ω D √ λ in terminology of43 p , / ¯ g w it hou t ± π T FIG. 19. The onset temperature of the pairing, obtained without including Matsubara frequencies ω m = ± πT . For γ = 2, this temperature is roughly 1 / T p, . Ref.[15], which is the same expression because λ = ¯ g /ω D ). The issue we discussed in thispaper is whether this T p is close to the actual superconducting T c , or T c is smaller, and thereis a range of pseudogap behavior between T c and T p . We argued that the actual T c scaleswith ω D and is much smaller than T p when ω D / ¯ g is small.To prove this, we solved the non-linear gap equation at T = 0 and ω D = 0 and founda continuum of solutions, governed by a single parameter ξ (0 ≤ ξ ≤ ∞ ). This in turngives rise to a continuum spectrum of condensation energy, E c,ξ , which can be viewed as acontinuum gapless spectrum of “longitudinal” gap fluctuations. An infinite set of the gapfunctions and the condensation energy exists already for γ <
2, but is a discrete one. For γ = 2, this spectrum becomes continuous in a manner similar how a discrete set of energylevels in a finite size crystal becomes a continuous vibration spectrum when system sizebecomes infinite. In our case, 1 / (2 − γ ) plays the role of a system size.Without the contribution from the gapless longitudinal branch, superfluid stiffness ρ s ( T =0) is larger than T p , and thermal corrections to superconducting order parameter scale ap-proximately as T /ρ s (0) and remain small at all T < T p . However, upon including contri-44utions from the longitudinal branch, we found that thermal corrections become of orderone already at much smaller T ∼ ω D . We identified this temperature with the actual super-conducting T c . We emphasize that T c vanishes at ω D = 0, and the behavior of the stiffnessdepends on the order in which the double limit ω D → T → γ = 2 model is critical at T = 0. At smaller γ , the ground state is notcritical at ω D = 0, and T c ∼ ¯ g (2 − γ ). It is finite but at γ ≤ T p ∼ ¯ g .We presented collaborative evidence that the γ = 2 model is critical, from the analysisof the continuum set of gap functions along real frequency axis and in the upper half-planeof frequency. We found that for each solution, there is an infinite array of 2 π vortices inthe upper frequency half-plane. The array of vortices stretches up to an infinite frequency,where each gap function from the continuous set has an essential singularity. We speculatedthat different gap functions from the continuous set are different extensions from the arrayof vortices, ending at an essential singularity, onto the upper half-plane of frequency.At a finite ω D , the set of gap functions becomes discrete and contains only a finite numberof solutions, all of which behave regularly in the high-frequency limit. The the number ofvortices also becomes finite. Still, at small ω D / ¯ g , the system behavior over a wide frequencyrange mimics that at ω D = 0.We showed the phase diagram of the γ = 2 model in variables T and ω D in Fig. 15 andthe phase diagram of the γ model at ω D = 0 in in variables T and γ in Fig. 16 In bothcases, there is range of pseudogap behavior between the onset temperature of the pairing T p and the actual T c . In the pseudogap region, the bound pairs are formed, but there is nomacroscopic phase coherence. We argued that in most of the pseudogap regime, the DoSand other observables display “gap filling” behavior, in which the peak position remains ata finite frequency up to T p , while the states below the peak gradually fill in.In the next (last) paper in the series we consider the behavior of the γ model for γ > T = 0, which gives rise to a reduction and eventualvanishing of the superfluid stiffness in the ground state.45 CKNOWLEDGMENTS
We thank I. Aleiner, B. Altshuler, E. Berg, D. Chowdhury, L. Classen, K. Efetov, R.Fernandes, A. Finkelstein, E. Fradkin, A. Georges, S. Hartnol, S. Karchu, S. Kivelson, I.Klebanov, A. Klein, R. Laughlin, S-S. Lee, G. Lonzarich, D. Maslov, F. Marsiglio, I. Mazin,M. Metlitski, W. Metzner, A. Millis, D. Mozyrsky, C. Pepan, V. Pokrovsky, N. Prokofiev,S. Raghu, S. Sachdev, T. Senthil, D. Scalapino, Y. Schattner, J. Schmalian, D. Son, G.Tarnopolsky, A-M Tremblay, A. Tsvelik, G. Torroba, E. Yuzbashyan, J. Zaanen, and partic-ularly R. Combescot and Y. Wang for useful discussions. The work by A.V.C. and Y.M.W.was supported by the NSF DMR-1834856. Y.-M.W, S.-S.Z.,and A.V.C acknowledge thehospitality of KITP at UCSB, where part of the work has been conducted. The research atKITP is supported by the National Science Foundation under Grant No. NSF PHY-1748958.A.V.C. also acknowledges the hospitality of Stanford University, where some results of thiswork have been obtained. His stay at Stanford has been supported through the Gordon andBetty Moore Foundation’s EPiQS Initiative, Grants GBMF4302 and GBMF8686.
Appendix A: Expansion in D ( ω m ) for γ = 2 and γ < In this Appendix we present some details of the analysis of the non-linear gap equationfor γ = 2 and elaborate on the claim in the main text that a continuous set of gap functionsexists only for γ = 2, while for smaller γ , the set is a discrete one. γ = 2 We begin with γ = 2. Consider first the limit of small frequencies ω (cid:28) ¯ g . For such ω m , ∆( ω m ) in the l.h.s. of the gap equation (2) can be neglected, as its inclusion leads toterms with extra ( ω m / ¯ g ) . This approximation is equivalent to neglecting ω m compared tothe self-energy Σ( ω m ) and is similar to the “no ω m ” approximation, used in the studies ofSYK-type models [49–51]. The non-linear gap equation at T = 0 without ∆( ω m ) in the l.h.sreduces to (cid:90) dω m (cid:48) D ( ω m (cid:48) ) − D ( ω m ) (cid:112) D ( ω m (cid:48) ) sign ω m | ω m − ω m (cid:48) | = 0 (A1)where, we recall, D ( ω m ) = ∆( ω m ) /ω m . 46he linearized gap equation is obtained from (A1) by neglecting D ( ω m (cid:48) ) in the denomi-nator. The exact solution of the linearized gap equation is Eq. (11): D ( ω m ) = 2 (cid:15) cos (cid:32) β log (cid:18) | ω m | ¯ g (cid:19) + φ (cid:15) (cid:33) sign ω (A2)where (cid:15) is an arbitrary overall factor, φ (cid:15) is yet undetermined constant, and β = 0 . πβ tanh( πβ ) = 1.We now expand Eq. (A1) in powers of D . We will be searching for the solution in theform D (cid:15) ( ω m ) = 2 ∞ (cid:88) n =0 (cid:15) n +1 Q n +1 cos (cid:32) (2 n + 1) (cid:32) β (cid:15) log (cid:18) | ω m | ¯ g (cid:19) (cid:33) + φ (cid:15) (cid:33) (A3)Substituting into (A1) and collecting contributions at each order in (cid:15) n +1 , we find that D (cid:15) ( ω m ) given by (A3) does satisfy Eq. (A1), and that all integrals are ultra-violet convergent,i.e., there is no need for regularization. The calculations are lengthy, but straightforward.We checked explicitly that β (cid:15) is the same in all terms in (A3) and is related to the original β by β (cid:15) = β (cid:18) − (cid:15) . (cid:15) + ... (cid:19) (A4)The numerical coefficients are Q = 0 .
222 + O ( (cid:15) ), Q = 0 .
043 + O ( (cid:15) ). We cited this resultand Eqs. (A3) and (A4) in Sec. III C.At larger frequencies, we need to keep ∆( ω m ) in the l.h.s. of (2). In the opposite limit ω m (cid:29) ¯ g , the leading term in D (cid:15) ( ω m ) is obtained by pulling 1 /ω from the integrand in ther.h.s. Then we obtain D (cid:15) ( ω m ) = a (cid:15) ω m (A5)where a (cid:15) = ¯ g (cid:90) dω m (cid:48) D ( ω m (cid:48) ) (cid:112) D ( ω m (cid:48) ) sign ω m (A6)Substituting this form of D ( ω m ) into the integrand in the r.h.s, we find that the integral isultra-violet convergent, i.e., the solution (A6) is self-consistent.The two solutions have to merge at ω m ∼ ¯ g . For the linearized gap equation (the limit (cid:15) →
0) we verified that this does happen for a certain value of φ (cid:15) in (A3). We conjecture thatthe same holds for other (cid:15) , i.e., for a certain φ (cid:15) , D (cid:15) ( ω m ) smoothly evolves between (A3) and(A6). We did similar analysis in Paper I. There, we demonstrated that for arbitrary φ (cid:15) , D (cid:15) ( ω )of Eq. (A3) approaches a constant at ω m → ∞ , while the desired term ( D (cid:15) ( ω m ) ∝ / | ω m | γ +1 γ ) is the subleading one. For a particular φ (cid:15) , a constant vanishes, and the high-frequency behavior becomes the expected one.We see from (A4) that β (cid:15) decreases with increasing (cid:15) , while the overall magnitude of∆( ω m ) increases. It is natural to expect that β (cid:15) = 0 at some critical (cid:15) = (cid:15) cr . We exploredthis in the main text.Another way to argue for the existence of (cid:15) cr is to depart from the opposite limit (cid:15) (cid:29) D ( ω m ) is supposed to be large. In this case, we introduce Ξ( ω m ) = 1 /D ( ω m ) andre-express the gap equation as (cid:90) dω m (cid:48) Ξ( ω m (cid:48) ) − Ξ( ω m ) (cid:112) ( ω m (cid:48) ) 1 | ω m − ω m (cid:48) | = 0 (A7)Note the absence of sign ω m in the integrand. At small Ξ, we neglect the Ξ ( ω m (cid:48) ) in thedenominator and search for the solution in the form Ξ( ω m ) = sign ω m | ω m / ¯ g | b . Substitutinginto (A7) we find b = ±
1, i.e., Ξ( ω m ) = (cid:18) A ¯ g | ω m | + A | ω m | ¯ g (cid:19) (A8)The first term does not satisfy the normalization condition and has to be discarded (seePaper I for the details on this). This leaves no parameter to adjust in order to match withthe behavior at high frequencies. This implies that there is no solution for the gap at large (cid:15) . There is a similarity between this analysis and the analysis in Paper I, where we consideredthe γ model for γ < N to make interactions inthe particle-hole and particle-particle channels non-equivalent. There, we found that thereexists N cr , which separates oscillating and non-oscillating solutions, and only oscillatingsolutions are compatible with high-frequency behavior. Here, (cid:15) plays the same role as N .We illustrate this in Fig. 20 γ < We now extend this approach to γ <
2. The gap equation for D ( ω m ) at ω m (cid:28) ¯ g has thesame form as in (A1), only | ω m − ω m (cid:48) | in the denominator is replaced by | ω m − ω m (cid:48) | γ . Thesolution of the linearized equation for D ( ω m ) is D ( ω m ) = 2 (cid:15) (cid:18) ¯ g | ω m | (cid:19) − γ/ cos (cid:18) β γ log (cid:18) | ω m | ¯ g (cid:19) γ + φ (cid:19) sign ω (A9)48 cr C β Nonzero gap NFL ground state N cr N β Nonzero gap NFL ground state (a) (b) β N β ϵ ϵ ϵ cr FIG. 20. The comparison between the behavior of β (cid:15) in the γ = 2 model and β N in the modelwith γ <
1, extended to
N > where β γ is some regular function of γ . As before, we search for the solutions in the form D (cid:15) ( ω m ) = 2 ∞ (cid:88) n =1 (cid:32) (cid:15) (cid:18) ¯ g | ω m | (cid:19) − γ/ (cid:33) n +1 Q n +1 cos (2 n + 1) (cid:18) β γ,(cid:15) log (cid:18) | ω m | ¯ g (cid:19) γ + φ (cid:15) (cid:19) (A10)Substituting into (A9), we find that the integrals that determine Q n +1 now contain infra-red divergencies. The only way to eliminate the divergencies is to assume that ∆ (cid:15) tends toa finite value at ω m →
0. But this is only possible for a discrete set of finite (cid:15) . We also notein passing that because the actual expansion parameter is (cid:15) (¯ g/ | ω m | ) − γ/ , the expansion of β γ,(cid:15) in powers of (cid:15) yields β γ,(cid:15) = β γ (1 + a ( (cid:15) (¯ g/ | ω m | ) − γ/ ) ). For a (cid:54) = 0, this gives rise toadditional terms, which are not matched by the terms in Eq. (A10). The only option thenis to set a = 0, i.e., leave β γ,(cid:15) equal to its bare value β γ .The outcome is that the continuous set of gap functions exists only for γ = 2. For smaller γ , this set is discrete. We also emphasize that the distinction between γ = 2 and γ < T = 0. At any finite T , the set of gap functions is a discrete one for all γ ≤
2, andthe solutions with different n from the set vanish at different temperatures T p, . In Fig. 21we show the results of high-accuracy numerical calculation of T p, for γ = 2. We see that T p, (cid:39) . × − ¯ g is finite. 49 IG. 21. Determination of the temperature T p, , at which the n = 1 solution develops in the γ = 2model. M c ∼ is the largest value of the Matsubara number, used in this numerical calculation.Extrapolating M c to ∞ yields a finite value T p, (cid:39) . × − ¯ g . Appendix B: The gap function ∆ ∞ The exact solution of the linearized gap equation at zero temperature has been derivedfor 0 < γ < < γ < γ = 2.We first solve for the gap function ∆ ∞ ( ω m ) along the Matsubara axis. Following thesame computational steps as in the analysis for γ <
2, we obtain D ∞ ( ω m ) = ∆ ∞ ( ω m ) /ω m in the form D ∞ ( ω m ) = (cid:15) ¯ gω m (cid:90) ∞−∞ dkb k e − ik log ( ω m / ¯ g ) , (B1)where (cid:15) is an infinitesimal number, b k = e − iI k [cosh( π ( k − β )) cosh( π ( k + β ))] / , (B2)and the function I k is I k = 12 (cid:90) ∞−∞ dk (cid:48) log | (cid:15) k (cid:48) − | tanh π ( k (cid:48) − k + iδ ) , (B3)50ere (cid:15) k (cid:48) = πk (cid:48) tanh( πk (cid:48) ) and β (cid:39) . πβ tanh ( πβ ) = 1 We cited theseresults in Eqs. (13), (14) and (15) in the main text.The integrals (B3) and (B1) can be computed numerically. We showed the result for D ∞ ( ω m ) in Fig.1 in the main text. The function D ∞ ( ω m ) oscillates at ω m < ¯ g and decaysas 1 / | ω m | at ω m > ¯ g .
1. Series expansion
The integral in Eq. (B3) can be evaluated by closing the integration contour along aninfinite arc in the complex plane of frequency. For | ω m | < ¯ g , the arc must be in the upperhalf-plane, and for ω m > ¯ g , in the lower half-plane. The integral is equal to the sum ofthe contributions from each pole of the function b k in the upper or lower half-plane. Theposition of these poles are obtained from the representation of b k as an infinite product ofthe Gamma-functions: (see Papers I and IV for details on this) b k = Γ(1 − ik )Γ(1 + ik ) Γ (cid:18)
12 + i ( k + β ) (cid:19) Γ (cid:18)
12 + i ( k − β ) (cid:19) ∞ (cid:89) m =1 Γ (cid:0) + i ( k − iβ m ) (cid:1) Γ (1 + m − ik )Γ (cid:0) − i ( k + iβ m ) (cid:1) Γ (1 + m + ik )(B4)Here β m > πβ m tan( πβ m ) = −
1. There is an infinite set of such β m ,specified by an integer m = 0 , , .. . Each β m is located within the interval 1 / m <β m < / m . Viewed as a function of complex k , b k has poles from individual Γ-functionsin the upper frequency half-plane at k = ± β + i ( n + 1 / n = 0 , , .. and at k = iβ m + i ( n + 1 / k = − i ( n + 1) and k = − i (1 + m + n ),where n = 0 , , , ... and m = 1 , , ... . | ω m | < ¯ g For | ω m | < ¯ g , the relevant poles are at k = ± β + i ( n + 1 /
2) and at k = iβ m + i ( n + 1 / n = 0 , , , ... . This yields series expansion for D ∞ ( y ) with y = ( | ω m | / ¯ g ) in the form D ∞ ( y ) = R e ∞ (cid:88) n =0 e ( iβ log y + φ ) C 2. At γ = 2, the result takes a very simple form C 1, i.e., | ω m | (cid:28) ¯ g , which is the rightbehavior of the gap function at small frequencies, see (B6).We note, however, that the first subleading term in (B8) scales as y sin( β log y ). Thiscontribution is smaller than the actual subleading term, which scales as y . and does notoscillate. This implies that, besides the leading term, the form of ∆ ∞ ( y ) is determined bynon-local corrections. | ω m | > ¯ g and logarithmic correction For y > 1, i.e., | ω m | > ¯ g , relevant poles are in the lower half-plane. These are not simplepoles, however. According to Eq. (B4), a pole at k = − i ( n + 1) ( n = 0 , , , ... ) is of order n + 1, namely a simple pole at n = 0, a double pole at n = 1, etc. The leading term inthe limit of | ω m | → ∞ is the contribution from a simple pole at k = − i ( n = 0), Thiscontribution accounts for 1 /y behavior of ∆ ∞ ( y ) at large y . However, the subleading termsfrom the rest poles contain extra logarithms on top of powers of 1 /y :∆ ∞ ( y ) = ∞ (cid:88) n =0 ˜ C >n y − n ) (log y ) n , (B9)To demonstrate the presence of the logs, consider as an example the contribution fromthe double pole at k = − i . We shift γ to 2 − δ , δ > δ → b k for γ ≤ z = − i and z = − (2 + δ/ i . In the neighborhood ofthe two poles, the function b k takes the form ∼ / ( z − z ) / ( z − z ). The contribution fromthe these two poles is obtained by circling out a loop C enclosing z and z . Evaluating heintegral and talking the limit δ → 0, we obtain √ y lim δ → + ¯ gω m (cid:73) C dz z + 2 i )( z + (2 + δ/ i ) e − iz log( y − δ/ ) = 2 π log yy . (B10)Similarly, the triple pole at k = − i gives rise to (log y ) /y , etc. Collecting the contributionsfrom every pole on the lower-half plane, we obtain (B9). 4. The universal oscillating term at large y We now show that the high-frequency form of D ∞ ( ω m ) contains an additional oscillatingcontribution. This contribution is exponentially small on the Matsubara axis, but, as we willsee, it becomes the dominant one on the real axis. To extract this contribution, we note thatfor large | ω m | / ¯ g , the argument of the cosine function, I k + k log y , passes through extremeat k ∼ k ∗ = y/π . Expanding around this point and evaluating the Gaussian integral, weobtain the universal piece D ∞ ; u ( y ) in the form D ∞ ; u ( y ) = 2 √ (cid:15)e − y cos (cid:20) ( π − π y + π (cid:21) . (B11)We see that D ∞ ; u ( y ) is exponentially small, yet this oscillating term is present. The total D ∞ ( y ) is the sum of (B9) and (B11). D ∞ ( y ) along real axis Let’s now transform from Matsubara to real axis. We use ω instead of y for bettertransparency. By construction, the gap function along the real axis is obtained by replacing iω m → ω + i + in the integrand in the r.h.s. of Eq. (B1). Under this transformation,log( | ω m | / ¯ g ) transforms into log( | ω | / ¯ g ) − iπ . The integral in Eq. (B1) splits into two parts: D ∞ ( ω ) = (cid:15) ¯ gω (cid:90) ∞ dk e − πk e − iI k − ik log( | ω | / ¯ g ) + e πk e iI k + ik log( | ω | / ¯ g ) (cid:112) cosh( π ( k − β )) cosh( π ( k + β )) . (B12)Evaluating each integral by expanding near the point where I k ± k log( | ω | / ¯ g ) passes throughextremum and approximating the denominator in (B12) by its form at large k , we find that53he first term is small in e − πk , while in the second term the exponential factor cancels out.Ignoring the first term, we obtain D ∞ ; u ( ω ) ≈ √ (cid:15)e iπ (cid:104) ( ω ¯ g ) +log ( ω ¯ g ) (cid:105) . (B13)Other contributions to D ∞ ( ω ) contain powers of ¯ g/ω and are smaller. As a result, on thereal axis, D ∞ ( ω ) ≈ D ∞ ; u ( ω ) at ω (cid:29) ¯ g .The same calculation can be carried out for an arbitrary complex frequency z = ω (cid:48) + iω (cid:48)(cid:48) in the upper frequency plane. For this, one has to replace iω m by z ≡| z | e iψ (0 < ψ < π ) inthe integrand in the r.h.s. of (B1). This changes log( | ω m | / ¯ g ) to log( | z | / ¯ g ) + i (2 ψ − π ) andgives ∆ ∞ ( z ), which we presented in Eq. (60) in the main text. Appendix C: Extended γ − model In Papers I-III and other works [47, 48, 52–58], we and others extended the γ modelto in-equal interactions in the particle-particle and particle-hole channel. This was doneby adding a factor 1 /N to the interaction in the particle-particle channel and leaving theinteraction in the particle-hole channel intact. The advantage of extending the model to N (cid:54) = 1 is that superconducting order in the ground state exists for N < N cr , while for larger N the ground state is a non-Fermi liquid. By analyzing the gap equation near this point,one can obtain useful information about how a discrete set of solutions emerges. In PaperIII we argued that the extension to N (cid:54) = 1 makes sense for γ < 1, while for γ ≥ 1, the modelwith N (cid:54) = 1 possess singularities, not present in the original γ model. We proposed anotherway to extend the model with γ > 1, which is free from singularities. The idea is to firstexplicitly cancel out singularities in the original γ model with γ > 1, and only then extendthe model to M (cid:54) = 1 by making interactions in the particle-particle and particle-hole channelin-equivalent. The extended model is then free from singularities, and one can obtain critical M cr , where superconducting order disappears at T = 0 (by our construction, it exists at M > M cr ). In this Appendix, we analyze the extended γ model for γ = 2. We show that acontinuous set of solutions for the gap equation emerges at M cr + 0.We first briefly describe the extension procedure. The two coupled Eliashberg equations54re for the pairing vertex Φ( ω m ) and the self-energy Σ( ω m ). For γ = 2, the equations areΦ( ω m ) = ¯ g πT (cid:88) m (cid:48) (cid:54) = m Φ( ω m (cid:48) ) (cid:113) ˜Σ ( ω m (cid:48) ) + Φ ( ω m (cid:48) ) 1 | ω m − ω m (cid:48) | , ˜Σ( ω m ) = ω m + ¯ g πT (cid:88) m (cid:48) (cid:54) = m ˜Σ( ω m (cid:48) ) (cid:113) ˜Σ ( ω m (cid:48) ) + Φ ( ω m (cid:48) ) 1 | ω m − ω m (cid:48) | (C1)where ˜Σ( ω m ) = ω m + Σ( ω m ). At T = 0, the r.h.s. of each of the two equations contains adivergent integral (cid:82) dx/x . To regularize the divergencies, we keep the temperature smallbut finite and set T = 0 at the end of calculations. At a finite T , the sum over m (cid:48) isnon-singular as singular self-action term with m (cid:48) = m cancels out by the same reason as thecontributions from non-magnetic impurities.We then introduce¯Φ( ω m ) = Φ( ω m ) − ¯ g ζ (2)(2 πT ) 1 (cid:113) ˜Σ ( ω m ) + Φ ( ω m ) ¯˜Σ( ω m ) = ˜Σ( ω m ) − ¯ g ζ (2)(2 πT ) 1 (cid:113) ˜Σ ( ω m ) + Φ ( ω m ) (C2)where ζ (2) = π / (cid:80) ∞ n =1 /n . Because Φ( ω m ) / ˜Σ( ω m ) = ¯Φ( ω m ) / ¯˜Σ( ω m ), Eqs. (C1) canbe re-expressed solely in terms of ¯Φ( ω m ) and ¯˜Σ( ω m ):¯Φ( ω m ) = ¯ g πT (cid:88) m (cid:48) (cid:54) = m ¯Φ( ω m (cid:48) ) (cid:113) ¯˜Σ ( ω m (cid:48) ) + ¯Φ ( ω m (cid:48) ) − ¯Φ( ω m ) (cid:113) ¯˜Σ ( ω m ) + ¯Φ ( ω m ) | ω m − ω m (cid:48) | , ¯˜Σ( ω m ) = ω m +¯ g πT (cid:88) m (cid:48) (cid:54) = m ¯˜Σ( ω m (cid:48) ) (cid:113) ¯˜Σ ( ω m (cid:48) ) + ¯Φ ( ω m (cid:48) ) − ¯˜Σ( ω m ) (cid:113) ¯˜Σ ( ω m ) + ¯Φ ( ω m ) | ω m − ω m (cid:48) | (C3)These equations are now free from singularities at T = 0, when the summation over Mat-subara numbers is replaced by the integration over ω m .We now extend the modified Eliashberg equations (C3) by multiplying the interaction inthe particle-particle channel by a factor 1 /M :¯Φ( ω m ) = ¯ g M (cid:90) dω (cid:48) m ¯Φ( ω (cid:48) m ) (cid:113) ¯˜Σ ( ω (cid:48) m ) + ¯Φ ( ω (cid:48) m ) − ¯Φ( ω m ) (cid:113) ¯˜Σ ( ω m ) + ¯Φ ( ω m ) | ω m − ω (cid:48) m | (C4)55he gap function ∆( ω m ) is expressed via ¯Φ( ω m ) and ¯˜Σ( ω m ) in the same way as via theoriginal Φ( ω m ) and ˜Σ( ω m ): ∆( ω m ) = ω m ¯Φ( ω m ) / ¯˜Σ( ω m ). The equation on ∆( ω m ) is∆( ω m ) = ¯ g M (cid:90) dω (cid:48) m | ω m − ω (cid:48) m | (cid:32) ∆( ω (cid:48) m ) − M ∆( ω m ) ω m ω (cid:48) m (cid:112) ∆ ( ω (cid:48) m ) + ( ω (cid:48) m ) − ∆( ω m )(1 − M ) (cid:112) ∆ ( ω m ) + ω m (cid:33) (C5)Absorbing 1 /M into ¯ g M = ¯ g /M , introducing a dimensionless ¯ ω m = ω m / ¯ g M , D (¯ ω m ) =∆(¯ ω m ) / ¯ ω m and re-arranging, we obtain from (C5) D (¯ ω m ) (cid:32) ¯ ω + 1 − M (cid:90) d ¯ ω (cid:48) m | ¯ ω m − ¯ ω (cid:48) m | (cid:32) (cid:112) D (¯ ω m ) − sign¯ ω (cid:48) m (cid:112) D (¯ ω (cid:48) m ) (cid:33)(cid:33) =12 (cid:90) d ¯ ω (cid:48) m | ¯ ω m − ¯ ω (cid:48) m | D (¯ ω (cid:48) m ) − D (¯ ω m ) (cid:112) D (¯ ω (cid:48) m ) sign¯ ω (cid:48) m (C6)Both integrals in (C6) are free from singularities and infra-red convergent.For infinitesimally small D (¯ ω m ), Eq. (C6) becomes D (¯ ω m ) (cid:18) ¯ ω m + 1 − M ¯ ω m (cid:19) = 12 (cid:90) d ¯ ω (cid:48) m D (¯ ω (cid:48) m ) − D (¯ ω m ) | ¯ ω m − ¯ ω (cid:48) m | sign ω (cid:48) m , (C7)At small ¯ ω m , the solution of the gap equation is D (¯ ω m ) = 2 (cid:15) cos (cid:0) β M log ¯ ω m + φ (cid:1) sign ω m , (C8)It has the same form as Eq. (11), but now β M = M/π . This form implies that M cr = 0.We now assume that M is small and solve the non-linear gap equation. Our key intensionis to check whether we still have a continuum of solutions. For this purpose, it is sufficientto focus on small ¯ ω m , when we can neglect bare ¯ ω m in the l.h.s. of (C6).As in Sec. III C, we search for the solution of (C6) in the series in (cid:15) for both D (¯ ω m ) and β . To leading order in M , we obtain D (¯ ω m ) = 2 (cid:15)π M / log ¯ ω m (cid:0) − (cid:15) + .... (cid:1) / (C9)where dots stand for (cid:15) and higher order terms. The M / dependence (same as ( M − M cr ) / as M cr = 0) is an expected one. The logarithmic dependence on frequency is consistent withthe result in Paper I, where we obtained log ¯ ω γm dependence at N = N cr . However, theresuch dependence exists only for N = N cr , while here we have an infinite set of solutions withthe same frequency dependence, but different amplitudes, parametrized by (cid:15) . All solutionsappear simultaneously at M = 0+ 56 complimentary piece of evidence for multiple solutions comes about if we simplify thel.h.s. of (C6) by dropping D terms in the denominator of the (1 − M ) term. The gapequation then reduces to D (¯ ω m ) (cid:18) ¯ ω + 1 − M ω (cid:19) = 12 (cid:90) d ¯ ω (cid:48) m | ¯ ω m − ¯ ω (cid:48) m | D (¯ ω (cid:48) m ) − D (¯ ω m ) (cid:112) D (¯ ω (cid:48) m ) sign¯ ω (cid:48) m (C10)The full gap equation for the original γ model is reproduced if we set M = 1, so Eq. (C10)can be viewed as another extension of the original model. The solution of the linearized gapequation is the same as before, with β M = M/π , hence M cr = 0. at M → 0. The expansionin (cid:15) now yields, to leading order in M ≡ M − M cr : D (¯ ω m ) = 2 (cid:15) (cid:18) cos f M (¯ ω m ) − (cid:15) M cos 3 f M (¯ ω m ) + ... (cid:19) (C11)where f M (¯ ω m ) = β M,(cid:15) log ¯ ω m + φ (C12)and β M,(cid:15) = β M (cid:18) − (cid:15) M + ... (cid:19) (C13)We see that the expansion holds in powers of (cid:15) /M and is valid up to (cid:15) ∼ ( M ) / , at which β M,(cid:15) vanishes. At larger (cid:15) , β M,(cid:15) becomes negative, and the solution disappears (there is nonormalized solution of the linearized gap equation). We see that there is again an infiniteset of solutions, specified by (cid:15) , which runs between 0 and (cid:15) cr = O ( √ M ). [1] A. Abanov and A. V. Chubukov, Interplay between superconductivity and non-fermi liquidat a quantum critical point in a metal. i. the γ model and its phase diagram at T = 0: Thecase 0 < γ < 1, Phys. Rev. B , 024524 (2020).[2] Y.-M. Wu, A. Abanov, Y. Wang, and A. V. Chubukov, Interplay between superconductivityand non-fermi liquid at a quantum critical point in a metal. ii. the γ model at a finite T for0 < γ < 1, Phys. Rev. B , 024525 (2020).[3] Y.-M. Wu, A. Abanov, and A. V. Chubukov, Interplay between superconductivity and non-fermi liquid behavior at a quantum critical point in a metal. iii. the γ model and its phasediagram across γ = 1, Phys. Rev. B , 094516 (2020). 4] Y.-M. Wu, S.-S. Zhang, A. Abanov, and A. V. Chubukov, Interplay between superconductivityand non-fermi liquid at a quantum critical point in a metal. iv. the γ model and its phasediagram at 1 < γ < 2, Phys. Rev. B , 024522 (2021).[5] T.-H. Lee, A. Chubukov, H. Miao, and G. Kotliar, Pairing mechanism in hund’s metal su-perconductors and the universality of the superconducting gap to critical temperature ratio,Phys. Rev. Lett. , 187003 (2018).[6] Y.-M. Wu, A. Abanov, and A. V. Chubukov, Pairing in quantum critical systems: Transitiontemperature, pairing gap, and their ratio, Phys. Rev. B , 014502 (2019).[7] A. Karakozov, E. Maksimov, and A. Mikhailovsky, The investigation of eliashberg equationsfor superconductors with strong electron-phonon interaction, Solid State Communications ,329 (1991).[8] F. Marsiglio and J. P. Carbotte, Gap function and density of states in the strong-couplinglimit for an electron-boson system, Phys. Rev. B , 5355 (1991), for more recent results see F.Marsiglio and J.P. Carbotte, “Electron-Phonon Superconductivity”, in “The Physics of Con-ventional and Unconventional Superconductors”, Bennemann and Ketterson eds., Springer-Verlag, (2006) and references therein; F. Marsiglio, Annals of Physics 417, 168102-1-23 (2020).[9] R. Combescot, Strong-coupling limit of eliashberg theory, Phys. Rev. B , 11625 (1995).[10] Y. Wang, Solvable strong-coupling quantum-dot model with a non-fermi-liquid pairing tran-sition, Phys. Rev. Lett. , 017002 (2020).[11] I. Esterlis and J. Schmalian, Cooper pairing of incoherent electrons: An electron-phononversion of the sachdev-ye-kitaev model, Phys. Rev. B , 115132 (2019).[12] D. Hauck, M. J. Klug, I. Esterlis, and J. Schmalian, Eliashberg equations for an electron-phonon version of the sachdev-ye-kitaev model: Pair breaking in non-fermi liquid supercon-ductors, Annals of Physics , 168120 (2020).[13] The model also describes strong coupling limit of the interaction between dispersion-lessfermions and phonons (SYK-Yukawa model) [10–12, 26, 59].[14] D. J. Scalapino, J. R. Schrieffer, and J. W. Wilkins, Strong-coupling superconductivity. i,Phys. Rev. , 263 (1966); D. J. Scalapino, The electron-phonon interaction and strong-coupling superconductors, Superconductivity, Ed. R.D. Parks (Dekker, New York, 1969) p (1969).[15] P. B. Allen and R. C. Dynes, Transition temperature of strong-coupled superconductors re- nalyzed, Phys. Rev. B , 905 (1975).[16] G. Bergmann and D. Rainer, The sensitivity of the transition temperature to changes in α f ( ω )., Z. Physik , 59 (1973).[17] D. Rainer, Principles of ab initio calculations of superconducting transition temperatures. In progress in low temperature physics , Chapter 4 (Elsevier, 1986) pp. 371–424.[18] P. Hertel, Die sprungtemperatur stark koppelnder supraleiter, Zeitschrift f¨ur Physik , 272(1971).[19] V. T. Geilikman and N. F. Masharov, Influence of the form of the phonon spectrum on thetransition temperature of superconductors, Journal of Low Temperature Physics , 131 (1972).[20] S. M. A.E. Karakozov, E.G. Maksimov, Effect of the frequency dependence of the electron-phonon interaction spectral function on the thermodynamic properties of superconductors,JETP , 971 (1975).[21] O. V. Dolgov, I. I. Mazin, A. A. Golubov, S. Y. Savrasov, and E. G. Maksimov, Criticaltemperature and enhanced isotope effect in the presence of paramagnons in phonon-mediatedsuperconductors, Phys. Rev. Lett. , 257003 (2005).[22] Y. Wang and A. Chubukov, Quantum-critical pairing in electron-doped cuprates, Phys. Rev.B , 024516 (2013).[23] F. Marsiglio, Eliashberg theory in the weak-coupling limit, Phys. Rev. B , 024523 (2018).[24] S. Mirabi, R. Boyack, and F. Marsiglio, Thermodynamics of eliashberg theory in the weak-coupling limit, Phys. Rev. B , 214505 (2020).[25] P. B. Allen and R. C. Dynes, Transition temperature of strong-coupled superconductors re-analyzed, Phys. Rev. B , 905 (1975).[26] A. V. Chubukov, A. Abanov, I. Esterlis, and S. A. Kivelson, Eliashberg theory of phonon-mediated superconductivity – when it is valid and how it breaks down, Annals of Physics ,168190 (2020).[27] This formula was originally obtained semi-analytically by Allen and Dynes [25]. They ex-pressed it as T p ∼ ω D √ λ to emphasize that at strong coupling T p becomes larger than ω D .Given that λ = (¯ g/ω D ) , their formula reduces to T p ∼ ¯ g .[28] C. A. R. S´a de Melo, M. Randeria, and J. R. Engelbrecht, Crossover from bcs to bose su-perconductivity: Transition temperature and time-dependent ginzburg-landau theory, Phys.Rev. Lett. , 3202 (1993). 29] D. J. Scalapino, S. R. White, and S. Zhang, Insulator, metal, or superconductor: The criteria,Phys. Rev. B , 7995 (1993).[30] L. Benfatto, A. Toschi, S. Caprara, and C. Castellani, Phase fluctuations in superconductors:From galilean invariant to quantum XY models, Phys. Rev. B , 140506 (2001).[31] L. Benfatto, S. Caprara, C. Castellani, A. Paramekanti, and M. Randeria, Phase fluctua-tions, dissipation, and superfluid stiffness in d-wave superconductors, Phys. Rev. B , 174513(2001).[32] S. G. Sharapov, V. P. Gusynin, and H. Beck, Effective action approach to the leggett’s modein two-band superconductors, The European Physical Journal B - Condensed Matter andComplex Systems , 45 (2002).[33] A. V. Chubukov, I. Eremin, and D. V. Efremov, Superconductivity versus bound-state for-mation in a two-band superconductor with small fermi energy: Applications to fe pnic-tides/chalcogenides and doped srtio , Phys. Rev. B , 174516 (2016).[34] D. J. Scalapino, Y. Wada, and J. C. Swihart, Strong-coupling superconductor at nonzerotemperature, Phys. Rev. Lett. , 102 (1965).[35] Y. Wada, The effect of quasiparticle damping on the ratio between the energy gap and thetransition temperature of lead, Rev. Mod. Phys. , 253 (1964).[36] F. Marsiglio, M. Schossmann, and J. P. Carbotte, Iterative analytic continuation of the electronself-energy to the real axis, Phys. Rev. B , 4965 (1988).[37] J. Bardeen and M. Stephen, Free-energy difference between normal and superconductingstates, Phys. Rev. , A1485 (1964).[38] R. Haslinger and A. V. Chubukov, Condensation energy in strongly coupled superconductors,Phys. Rev. B , 214508 (2003).[39] E. A. Yuzbashyan, A. V. Chubukov, A. Abanov, and B. L. Altshuler, Non BCS pairing neara quantum phase transition – effective mapping to a spin chain. in preparation, .[40] The second term in (37) diverges logarithmically at ξ = 0 if we use D ( ω ) ≈ ∆(0) /ω atsmall frequencies. This divergence comes from the putative normal state energy E ∆=0 whilethe ground state energy E ∆ remains finite. For ξ > 0, both E ∆=0 and E ∆ have logarithmicsingularities, which cancel out in E c = E ∆ − E ∆=0 .[41] For a generic γ , ∆ ( ω m ) and ∆ (cid:48)(cid:48) ( ω ) are related by Cauchy formula: ∆ ( ω m ) =(1 /π ) (cid:82) ∞ dω ∆ (cid:48)(cid:48) ( ω ) ω/ ( ω + ω m ). For γ < 2, typical ω are of order ω m , and to reproduce ( ω m ) ∝ / | ω m | γ one needs ∆ (cid:48)(cid:48) ( ω ) = sin( πγ/ 2) sign ω/ | ω | γ . The case γ = 2 is an exceptionhere because 1 /ω m dependence of ∆ ( ω m ) is obtained by pulling 1 /ω m out of the denominatorin the Cauchy formula. The remaining integral is determined by ω = O (¯ g ) rather than O ( ω m ).Because of this, the fact that ∆ ( ω m ) ∝ /ω m at large frequencies does not imply that ∆ (cid:48)(cid:48) ( ω )must behave as 1 /ω .[42] D. Mozyrsky and A. V. Chubukov, Dynamic properties of superconductors: Anderson-bogoliubov mode and berry phase in the bcs and bec regimes, Phys. Rev. B , 174510(2019).[43] J. Schrieffer, in Theory of Superconductivity (Benjamin, New York, 1964); V. J. Emery andS. A. Kivelson, Importance of phase fluctuations in superconductors with small superfluiddensity, Nature , 434 (1995).[44] V. Pokrovsky, Properties of ordered, continuously degenerate systems, Advances in Physics , 595 (1979).[45] K. Maki, The magnetic properties of superconducting alloys i, Physics Physique Fizika , 21(1964).[46] K. Maki, Theory of electron-spin resonance in gapless superconductors, Phys. Rev. B , 191(1973).[47] A. Abanov, Y.-M. Wu, Y. Wang, and A. V. Chubukov, Superconductivity above a quantumcritical point in a metal: Gap closing versus gap filling, fermi arcs, and pseudogap behavior,Phys. Rev. B , 180506 (2019).[48] Y.-M. Wu, A. Abanov, Y. Wang, and A. V. Chubukov, Special role of the first matsubarafrequency for superconductivity near a quantum critical point: Nonlinear gap equation below T c and spectral properties in real frequencies, Phys. Rev. B , 144512 (2019).[49] Y. Gu, A. Kitaev, S. Sachdev, and G. Tarnopolsky, Notes on the complex sachdev-ye-kitaevmodel, Journal of High Energy Physics , 157 (2020).[50] A. A. Patel and S. Sachdev, Theory of a planckian metal, Phys. Rev. Lett. , 066601 (2019).[51] Y. Wang and A. V. Chubukov, Quantum phase transition in the yukawa-syk model, Phys.Rev. Research , 033084 (2020).[52] S. Raghu, G. Torroba, and H. Wang, Metallic quantum critical points with finite bcs couplings,Phys. Rev. B , 205104 (2015).[53] Y. Wang, A. Abanov, B. L. Altshuler, E. A. Yuzbashyan, and A. V. Chubukov, Superconduc- ivity near a quantum-critical point: The special role of the first matsubara frequency, Phys.Rev. Lett. , 157001 (2016).[54] H. Wang, S. Raghu, and G. Torroba, Non-fermi-liquid superconductivity: Eliashberg approachversus the renormalization group, Phys. Rev. B , 165137 (2017).[55] H. Wang, Y. Wang, and G. Torroba, Superconductivity versus quantum criticality: Effects ofthermal fluctuations, Phys. Rev. B , 054502 (2018).[56] J. A. Damia, S. Kachru, S. Raghu, and G. Torroba, Two-dimensional non-fermi-liquid metals:A solvable large- n limit, Phys. Rev. Lett. , 096402 (2019).[57] A. V. Chubukov, A. Abanov, Y. Wang, and Y.-M. Wu, The interplay between superconduc-tivity and non-fermi liquid at a quantum-critical point in a metal, Annals of Physics ,168142 (2020).[58] J. A. Damia, M. Sol´ıs, and G. Torroba, How non-fermi liquids cure their infrared divergences,Phys. Rev. B , 045147 (2020).[59] L. Classen and A. V. Chubukov, Superconductivity of incoherent electrons in yukawa-sykmodel (2021)., 045147 (2020).[59] L. Classen and A. V. Chubukov, Superconductivity of incoherent electrons in yukawa-sykmodel (2021).