Interplay between superconductivity and non-Fermi liquid at a quantum critical point in a metal. II. The γ-model at a finite T for 0<γ<1
IInterplay between superconductivity and non-Fermi liquid at aquantum critical point in a metal. II. The γ -model at a finite T for < γ < . Yi-Ming Wu, Artem Abanov, Yuxuan Wang, 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 Department of Physics, University of Florida, Gainesville, FL 32611, USA (Dated: August 4, 2020)
Abstract
In this paper we continue the analysis of the interplay between non-Fermi liquid and supercon-ductivity for quantum-critical systems, the low-energy physics of which is described by an effectivemodel with dynamical electron-electron interaction V (Ω m ) ∝ / | Ω m | γ (the γ model). In paper I [A.Abanov and A. V. Chubukov, arXiv:2004.13220 (2020)] two of us analyzed the γ model at T = 0for 0 < γ < n ( ω m ) forthe n th solution changes sign n times as the function of Matsubara frequency. In this paper weanalyze the linearized gap equation at a finite T . We show that there exist an infinite set of pairinginstability temperatures, T p,n , and the eigenfunction ∆ n ( ω m ) changes sign n times as a functionof a Matsubara number m . We argue that ∆ n ( ω m ) retains its functional form below T p,n and at T = 0 coincides with the n th solution of the nonlinear gap equation. Like in paper I, we extendthe model to the case when the interaction in the pairing channel has an additional factor 1 /N compared to that in the particle-hole channel. We show that T p, remains finite at large N dueto special properties of fermions with Matsubara frequencies ± πT , but all other T p,n terminate at N cr = O (1). The gap function vanishes at T → N > N cr and remains finite for N < N cr .This is consistent with the T = 0 analysis. a r X i v : . [ c ond - m a t . s up r- c on ] A ug . INTRODUCTION In this paper we continue our analysis of the competition between non-Fermi-liquid (NFL)physics and superconductivity near a quantum-critical point (QCP) in a metal for a classof models with dynamical four-fermion interaction V ( q, Ω m ), mediated by a critical softboson. The interaction V ( q, Ω m ) gives rise to strong fermionic self-energy, which destroysFermi-liquid behavior, in most cases in dimensions D ≤
3, and also gives rise to an attrac-tion in at least one pairing channel. We consider the class of models in which bosons areslow excitations compared to fermions (like phonons in the case when the Debye frequencyis much smaller than the Fermi energy). In this situation, the momentum integration inthe expressions for the fermionic self-energy and the pairing vertex in the proper attractivechannel can be carried out, and both the self-energy and the pairing vertex are determinedby the effective, purely dynamical local interaction V (Ω m ) ∝ R d q V ( q, Ω), with q connect-ing points on the Fermi surface. The coupled equations for the self-energy Σ( k F , ω m ) andthe pairing vertex Φ( ω m ) have the same structure as Eliashberg equations for a phononsuperconductor. We consider quantum-critical models for which V (Ω m ) = ¯ g γ / | Ω m | γ (the γ model). In the previous paper, Ref. , hereafter referred to as paper I, we listed a number ofquantum-critical systems, the low-energy physics of which is described by the γ model withdifferent γ . This paper also contains an extensive list of references to earlier publications onthis subject.The interaction V (Ω m ) is singular, and gives rise to two opposite tendencies: NFL behav-ior in the normal state, with Σ( ω m ) ∝ ω − γm , and the pairing. The two tendencies competewith each other as a NFL self-energy reduces the magnitude of the pairing kernel, while thefeedback from the pairing reduces fermionic self-energy.In paper I we analyzed zero-temperature behavior of the γ model for 0 < γ <
1. We foundthat the system does become unstable towards pairing, i.e., in the ground state the pairingvertex Φ( ω m ) and the pairing gap ∆( ω m ) = Φ( ω m ) / (1 + Σ( ω m ) /ω m ) are finite. However,in qualitative distinction with BCS/Eliashberg theory of superconductivity, in which thereis a single solution of the gap equation, here we found an infinite discrete set of solutions∆ n ( ω m ). The solutions have the same spatial gap symmetry, but ∆ n ( ω m ) changes sign n times as a function of Matsubara frequency. The gap functions ∆ n ( ω m ) tend to finite valuesat zero frequency, but the magnitude of ∆ n (0) decreases with n and at large enough n scales2s ∆ n (0) ∝ e − An , where A is a function of γ . The solutions with different n are topologicallydifferent. To see this, one has to consider positive ω m and analytically continue the function∆ n ( ω m ) into the upper half-plane of frequency ( ω m → − iz ). The gap function ∆ n ( z ) is acomplex function for a generic z , ∆( z ) = ∆ ( z ) + i ∆ ( z ), and one can define z − dependentphase as ∆( z ) = | ∆( z ) | e iψ ( z ) . One can verify that ψ ( z ) winds up by 2 π upon circling everynodal point on the Matsubara axis, i.e., ∆ n ( z ) is a solution with n vortices. This samephysics can be detected along the real axis ( z = ω ): the total change of the phase between ω = −∞ and ∞ for ∆ n ( ω ) has extra 2 πn compared to the case n = 0. In this paper we analyze the γ model for 0 < γ < T . We solve the linearizedgap equation and show that there is an infinite, discrete set of the onset temperatures forthe pairing T p,n . The magnitude of T p,n decreases with n and at large n , T p,n ∝ e − An withthe same A as for ∆ n (0) from the T = 0 analysis. The gap function ∆ n ( ω m ) at T = T p,n − n times as a function of discrete Matsubara frequency ω m = πT (2 m + 1). Weargue that ∆ n ( ω m ), which develops at T p,n , retains its functional form with n sign changesalso at T < T p,n , and at T → n th solution of the full nonlinear gapequation, which we obtained in paper I.In paper I we also added an extra parameter N to the γ model to control the ratio ofinteractions in the particle-particle and the particle-hole channels. Specifically, we added thefactor 1 /N to the interaction in the pairing channel and left the interaction in the particle-hole channel intact. We found that at T = 0 there exists a critical N cr >
1, which separatesthe NFL ground state for
N > N cr and the state with a nonzero ∆ for N < N cr . We foundthat all ∆ n ( ω m ) emerge simultaneously at N = N cr −
0. Here we analyze the extended γ model at a finite T . We show that the onset temperature T p, for sign-preserving gap function∆ ( ω m ) remains finite for all N and at large N scales as T p, ∝ /N /γ . This T p, has beenstudied before , and its existence even for very large N was attributed to special propertiesof fermions with ω m = ± πT , for which nonthermal part of the self-energy vanishes. Namely,it was argued that the pairing at T p, at large N is predominantly between these fermions,and ∆ n =0 ( ± πT ) is much larger than at other ω m .Here, we show that there exists an infinite set of other T p,n with n >
0, which terminate at N = N cr and are related to the solutions of the nonlinear gap equation at T = 0, ∆ n> ( ω m ).We show that the structure of the eigenfunctions ∆ n> ( ω m ) at T p,n is opposite to that of∆ ( ω m ) in the sense that ∆ n> ( ω m ) has additional smallness at ω m = ± πT . We show that3or any N < N cr , T p,n gradually decrease with n , and the set become more dense as n increases. As the consequence, in the limit T →
0, the density of eigenvalues (DoE) remainsfinite for all
N < N cr . We show that it is the same at T = 0 and at T →
0, i.e., the systemevolves smoothly between T = 0 and T >
0. For completeness, we obtain the DoE awayfrom a QCP, when a pairing boson has a finite mass. We show that in this case the solutionsat T → N .The structure of the paper is the following. In Sec. II we present Eliashberg equationsfor the γ model and discuss in some detail the extension to N = 1 at a finite T . In Sec. IIIwe analyze analytically and numerically the linearized gap equation at a finite T . We showthat there are two possibilities to get a solution at small T : one can either set N ∝ /T γ (the n = 0 solution), or keep N finite. We consider these two cases separately in Secs. III Aand III B. In Sec. III A we find the analytic solution for T p, , which reproduces earlier resultsin Refs. . In Sec. III B we first solve numerically the gap equation at T →
0, obtain thedensity of eigenvalues, and show that it is nonzero for all N between N = N cr and 0. Thisis consistent with the T = 0 result in paper I. We then extend calculations to finite T andshow that for any N < N cr (and any γ from 0 < γ < T p,n . We argue that all T p,n with n > N = N cr and that T p,n at large n should decay exponentially with n as T p,n ∝ e − An .We relate A to the form of ∆ n ( ω m ) at T = 0, express A via γ and N , and argue that T p,n is proportional to ∆ n (0) at T = 0. In Sec. III C we present numerical results for the gapfunction ∆ n ( ω m ) at T = T p,n −
0. We show that ∆ n ( ω m ) changes sign n times as a functionof discrete Matsubara frequency. We extend the numerical analysis to a range of T ≤ T p,n ,where ∆ n ( ω m ) is small but finite, and show that ∆ n ( ω m ) still changes sign n times, i.e., itstopological structure does not change below T p,n . We argue that at T → n ( ω m ) approaches the n th solution of the nonlinear gapequation at T = 0, which has the same topological structure. In Sec. IV we discuss howthe density of eigenvalues at T → T p,n evolve away from aQCP, when the pairing boson acquires a finite mass M . Finally, in Sec. V we present ourconclusions and briefly outline what we will do next.Before we proceed, we briefly outline our reasoning to look for solutions with different ∆ n .For γ <
1, which we study here, the set of temperatures T p,n is infinite, but discrete, and T p, for n = 0 is the largest onset temperature for the pairing. there is just one phase transition,4t the largest T p, ( N ). The solutions with n > T , but in the Eliashberg approximation, there is justone transition temperature, and the gap function has no nodes along the Matsubara axis.However, as γ increases, condensation energies of the solutions for n > n = 0, and, finally, the condensation energies for all finite n become thesame as for n = 0, while the ones for n → ∞ form a continuous spectrum. This happens for γ ≥
1, when N = 1 and for γ = 2 for the original γ model with N = 1. Once the spectrumbecomes continuous, a new channel for low-energy longitudinal fluctuations emerges, with anumber of consequences. The key goal of the study of the gamma model is to demonstratethat it has this new physics. However, to show this, we first have to establish the presenceof the discrete set of solutions for smaller γ , and, most important, to prove that the set isinfinite and extends to n = ∞ (the solution of the linearized gap equation). This is whatwe do in this paper (and in paper I), where we study the gamma model for smaller γ < II. THE γ -MODEL, ELIASHBERG EQUATIONS The γ model was introduced in paper I and in earlier papers, and we refer the readerto these works for the justification of the model and its relation to various quantum-criticalsystems, e.g., fermions at the verge of a nematic of spin-density wave/charge-density wavetransition, or Sachdev-Ye-Kitaev-type models of dispersion-less fermions randomly inter-acting with optical phonons. The model describes low-energy fermions, interacting by ex-changing soft bosonic excitations. The effective dynamical four-fermion interaction, whichcontributes to both the fermionic self-energy Σ( k F , ω m ) = Σ( ω m ) and the pairing vertexΦ( ω m ), is V (Ω m ) = ¯ g γ / | Ω m | γ with the exponent γ . [For Φ( ω ), this holds in an attractivechannel with a particular spatial symmetry, determined by the underlying model with fullmomentum dependence of the interaction]. We assume, like in earlier works, that soft bosonsare slow compared to dressed fermions. In this situation, the coupled equations for Σ( ω m )and Φ( ω m ) are similar to Eliashberg equations for the case when the effective dynamicalinteraction is mediated by phonons and we will use the term “Eliashberg equations” for ourcase. 5t a finite T the coupled Eliashberg equations for Φ( ω m ) and Σ( ω m ) are, in Matsubaraformalism Φ( ω m ) = ¯ g γ πT X m Φ( ω m ) q ˜Σ ( ω m ) + Φ ( ω m ) 1 | ω m − ω m | γ , ˜Σ( ω m ) = ω m + ¯ g γ πT X m ˜Σ( ω m ) q ˜Σ ( ω m ) + Φ ( ω m ) 1 | ω m − ω m | γ (1)where ˜Σ( ω m ) = ω m + Σ( ω m ). Observe that we define Σ( ω m ) with the overall plus sign andwithout the overall factor of i , that is G − ( k, ω m ) = i ˜Σ( ω m ) − (cid:15) k . In these notations, Σ( ω m )is a real function, odd in frequency.The superconducting gap function ∆( ω m ) is defined as∆( ω m ) = ω m Φ( ω m )˜Σ( ω m ) = Φ( ω m )1 + Σ( ω m ) /ω m (2)The equation for ∆( ω m ) is readily obtained from (1):∆( ω m ) = ¯ g γ πT X m ∆( ω m ) − ∆( ω m ) ω m ω m q ( ω m ) + ∆ ( ω m ) 1 | ω m − ω m | γ . (3)This equation contains a single function ∆( ω m ), but at the cost that ∆( ω m ) appears also inthe right-hand-side (rhs). Both Φ( ω m ) and ∆( ω m ) are defined up to an overall U (1) phasefactor, which we set to zero.The full set of Eliashberg equations for electron-mediated pairing contains the additionalequation, which describes the feedback from the pairing on V (Ω). In this paper we do notinclude this feedback into consideration, to avoid additional complications. In general terms,the feedback from the pairing makes bosons less incoherent and can be modeled by takingthe exponent γ to be T dependent and by moving it towards larger values as T decreases.Both Φ( ω m ) and Σ( ω m ) contain divergent contributions from the m = m terms in thesummation over internal Matsubara frequencies. The divergencies, however, cancel outin the gap equation (3) by the Anderson theorem , because scattering with zero-frequencytransfer mimics the effect of scattering by nonmagnetic impurities. To get rid of divergenciesin (1), one can use the same trick as in the derivation of the Anderson theorem: pull outthe terms with m = m from the frequency sums, move them to the left-hand-side (lhs). ofthe equations, and introduce new variables Φ ∗ ( ω m ) and Σ ∗ ( ω m ) asΦ ∗ ( ω m ) = Φ( ω m ) (1 − Q ( ω m )) , ˜Σ ∗ ( ω m ) = ˜Σ( ω m ) (1 − Q ( ω m )) (4)6here Q ( ω m ) = πT V (0) q ˜Σ ( ω m ) + Φ ( ω m ) (5)The ratio Φ( ω m ) / ˜Σ( ω m ) = Φ ∗ ( ω m ) / ˜Σ ∗ ( ω m ), hence ∆( ω m ), defined in (2), is invariant underthis transformation. One can easily verify that the equations on Φ ∗ ( ω m ) and ˜Σ ∗ ( ω m ) are thesame as in (1), but the summation over m now excludes the divergent term with m = m .To simplify the formulas, we re-define Φ ∗ and ˜Σ ∗ back as Φ and ˜Σ. We haveΦ( ω m ) = ¯ g γ πT X m = m Φ( ω m ) q ˜Σ ( ω m ) + Φ ( ω m ) 1 | ω m − ω m | γ , ˜Σ( ω m ) = ω m + ¯ g γ πT X m = m ˜Σ( ω m ) q ˜Σ ( ω m ) + Φ ( ω m ) 1 | ω m − ω m | γ (6)In the normal state (Φ ≡ m > ω m ) = ¯ g γ (2 πT ) − γ m X m =1 | m | γ = ¯ g γ (2 πT ) − γ H m,γ (7)where H m,γ is the harmonic number. This expression holds for ω m = ± πT . For the twolowest Matsubara frequencies, Σ( ± πT ) = 0. This will be essential for our paper. We remindthe reader that Σ( ω m ) in (7) is not the full self-energy, as the summation is only over m = m The full self-energy includes also the term with m = m and is nonzero for all Matsubarafrequencies.We now extend the model and introduce a parameter N to control the relative strengthof the interactions in the particle-hole and particle-particle channels. In essence we justassume that the effective, momentum-integrated interactions in the two channels have thesame frequency dependence, but different amplitudes. We treat N as a continuous variable,but keep the same notation as in the original work, Ref. ), where the original model wasextended to a matrix SU ( N ) model, and N was treated as an integer. With this extensionΦ( ω m ) = ¯ g γ N πT X m = m Φ( ω m ) q ˜Σ ( ω m ) + Φ ( ω m ) 1 | ω m − ω m | γ , ˜Σ( ω m ) = ω m + ¯ g γ πT X m = m ˜Σ( ω m ) q ˜Σ ( ω m ) + Φ ( ω m ) 1 | ω m − ω m | γ (8)and ∆( ω m ) = ¯ g γ N πT X m = m ∆( ω m ) − N ∆( ω m ) ω m ω m q ( ω m ) + ∆ ( ω m ) 1 | ω m − ω m | γ . (9)7e emphasize that for our paper the extension to N = 1 is just a convenient way to betterunderstand the interplay between the tendencies towards NFL and pairing. Our ultimategoal is to understand the physics in the physical case of N = 1. This is why we firsteliminated the divergent terms, which cancel out in the gap equation for N = 1, and onlythen extended the model to N = 1. An alternative approach, suggested in Ref. , is toextend to N = 1 without first subtracting the m = m terms in the equations for Φ( ω m ) andΣ( ω m ). In this case, one has to deal with the actual divergencies in the rhs of these equationsand also in the gap equation. Using the same analogy as before, once the interactions inthe particle-hole and the particle-particle channel become nonequivalent (the case N = 1),thermal fluctuations effectively act as singular magnetic impurities, which do not cancel outin the gap equation .We will chiefly analyze the linearized equation for the pairing vertex Φ( ω m ) N Φ( ω m ) = ¯ g γ πT X m = m Φ( ω m ) | ˜Σ( ω m ) | | ω m − ω m | γ (10)and the gap function∆( ω m ) = ¯ g γ N πT X m = m ∆( ω m ) − N ∆( ω m ) ω m ω m | ω m | | ω m − ω m | γ . (11)The two equations are equivalent in the sense that both have a nontrivial solution at the onsettemperature for the pairing. We label this temperature as T p rather than superconducting T c . because in general T c is smaller than T p due to gap fluctuations. In this paper we focuson the solution of the Eliashberg equations and neglect gap fluctuations, in which case, T p = T c . The comparative analysis of T p and T c will be presented in a separate paper.To distinguish between finite T and T = 0, below we label the gap function at a finite T as ∆ n ( m ), where m is a discrete Matsubara number and ω m = πT (2 m + 1), and at T = 0as ∆( ω m ), where ω m is a continuous Matsubara frequency. We use the same notations forother quantities. III. THE LINEARIZED GAP EQUATION, ANALYTIC CONSIDERATION
We re-write the self-energy (7) asΣ( m ) = πT KA ( m )sgn(2 m + 1) (12)8here A ( m ) = 2 m X n γ , m > A ( m ) = A ( − m − , m < − A (0) = A ( −
1) = 0 (13)and K = (cid:18) ¯ g πT (cid:19) γ (14)Using these notations and the fact that Φ( m ) is even under m → − m −
1, we re-write Eq.(10) for Φ( m ) as the set of two coupled equations, by singling out Φ(0) = Φ( −
1) (Ref. ):( N − K )Φ(0) = ∞ X n =1 Φ( n ) A ( n ) + n +1 K n γ + 1( n + 1) γ ! (15)Φ( m >
0) = 1 N ∞ X n =1 ,n = m Φ( n ) A ( n ) + n +1 K | n − m | γ + 1( n + m + 1) γ ! + Φ(0) KN n γ + 1( n + 1) γ ! (16)Eliminating Φ(0) from this set we obtain N Φ( m >
0) = ∞ X n =1 ,n = m Φ( n ) A ( n ) + n +1 K | n − m | γ + ∞ X n =1 Φ( n ) A ( n ) + n +1 K n + m + 1) γ + KN − K ∞ X n =1 Φ( n ) A ( n ) + n +1 K n γ + 1( n + 1) γ ! m γ + 1( m + 1) γ ! (17)A quick inspection of Eq. (17) shows that at low T , when K (cid:29)
1, there are two regions of N , in which the solution with Φ( m ) = 0 may appear: a) N ≈ K , such that N − K = const,and b) N = O (1), such that K (cid:29) N . We consider these two regions separately. A. The region of large N , N ≈ K We introduce b γ via N as N = K + b γ , substitute this relation into (17), and take thelimit K → ∞ . The divergent K cancels out, and we obtain b γ Φ( m >
0) = ∞ X n =1 Φ( n ) A ( n ) n γ + 1( n + 1) γ ! m γ + 1( m + 1) γ ! (18)The form of the kernel implies that the solution should be in the formΦ( m >
0) = C m γ + 1( m + 1) γ ! (19)Substituting this form into (18) we obtain that b γ is given by b γ = ∞ X n =1 A ( n ) n γ + 1( n + 1) γ ! (20)9 .2 0.4 0.6 0.8 1.024681012 γ b γ FIG. 1. The plot of b γ vs γ , Eq.(20). At large n , A ( n ) ≈ n − γ / (1 − γ ). The sum in the rhs of (20) then converges for all 0 < γ < b γ is indeed of order 1. We plot b γ in Fig.1 At small γ , b γ ∝ /γ . At γ = 1, b = 1 . T p = (cid:18) ¯ g π (cid:19) N /γ b γ N γ ! (21)We emphasize that T p remains finite for arbitrary large N , i.e., for arbitrary weak pairinginteraction. Eq. (21) has been earlier obtained in Ref. using somewhat different computa-tional procedure.Using (15), (19), and (20), we find Φ(0) = C . Comparing with (19), we see that themagnitude of the pairing vertex at Matsubara frequencies ω m = ± πT is of the same orderas at other frequencies [a more detailed analysis shows that Φ( m ) is the largest at m = 2,see Fig. 8]. At the same time, the pairing gap ∆( m ) is strongly peaked at m = 0 , − −
1) = C , while ∆( m >
0) = ∆(0) /K is much smaller. In fact, the leading termin the expression for T p in Eq. (21) can be obtained from the gap equation (9) if we keepthere only the terms with m = 0 , − P m sign(2 m + 1) / | m | γ = 0. Inthis case, the gap equation simplifies to∆(0) = ¯ g γ N ∆( − πT ) γ ∆( −
1) = ¯ g γ N ∆(0)(2 πT ) γ (22)Solving this set we obtain T p = ¯ g/ (2 πN /γ ). The outcome is that at large N , the pairingpredominantly involves fermions with Matsubara frequencies ω m = ± πT (Matsubara num-10 a) (b) FIG. 2. The gap function and the density of states below T p, , for N > N cr . We set γ = 0 .
9, inwhich case N cr = 1 . m = 0) ≡ ∆( πT ) as a function of T . Panel (b): the densityof states ν ( ω ) in real frequency, normalized to its value in the normal state. The density of statesshows a gap-filling behavior: as T increases towards T p, , states at low frequencies fill in, and theposition of the maximum shifts to higher ω bers m = 0 , − ω m ) in (12) vanishes. The gap function at other Matsubarafrequencies is induced by a nonzero ∆(0) = ∆( − n , which is thenumber of times Φ m changes sign as a function of Matsubara number m . In this classification,the sign-preserving pairing vertex Φ( m ) in (19) corresponds to n = 0. Accordingly, T p in(21) is T p, and Φ( m ) = Φ ( m ), ∆( m ) = ∆ ( m ).At a first glance, the result that T p, is nonzero at large N contradicts the T = 0 analysisin paper I and earlier studies that at T = 0 there exists critical N = N cr , whichseparates the state with a finite ∆( ω m ) for N < N cr and the NFL normal state for N > N cr .We, however, showed in Refs. that Φ( m ) and ∆( m ) evolve nonmonotonically below T p, and for N > N cr vanish at T →
0. We show the behavior of ∆( m = 0) ≡ ∆( πT ) vstemperature in Fig.2a. Such a state has a number of other unusual properties, e.g., thegap-filling behavior of the density of states and the spectral function (see Fig.2b). We referthe reader to Refs. for the detailed analysis of the Eliashberg equations below T p, and ofthe role of gap fluctuations. In Appendix A we compute the free energy F ,p of a state witha finite ∆ ( m ) for N > N cr and analyze δF = F p, − F n , where F n is the free energy at∆ ( m ) = 0. We use δF = δE − T δS and study the two terms separately. We find that at T ≤ T p, , negative δF comes primarily from negative δE , while − T δS is positive. However,below a certain T , close to where ∆( m = 0) changes its behavior in Fig. 2, − T δS becomes11 cr N = 1 N cr N = 1 T p ,0 T p ,1 T p ,2 T p ,3 N T N T T p ,0 n = 0 n = 1 n = 2 n = 3 . . . T p ,4 FIG. 3. Two potential phase diagrams of the γ model. (a) There exists only one onset temperature T p, for any given N , like in BCS theory. In this case, N cr is an isolated T = 0 QCP, withno transition line attached to it. (b) There exists an infinite set of T p,n , which all terminate at N = N cr . In this case, a QCP at N = N cr is critical point of an infinite order. We show that thecorrect phase diagram is the one in panel (b). negative, and δE almost vanishes. In this T range, a negative δF comes primarily from theentropy. This is consistent with the vanishing of ∆ ( m ) at T = 0.Now, if T p, was the only solution of the linearized gap equation, the phase diagram withinthe Eliashberg theory would be as in Fig.3(a), i.e., there would be an isolated QCP at T = 0and N = N cr , with no transition line coming out of it. We show below that this is not thecase, and the actual phase diagram is the one in Fig.3(b), with an infinite number of linesterminating at N cr . For this to hold, the linearized gap equation must have an infinite setof onset temperatures T p,n for any N < N cr . This is what we analyze next. B. The region N = O (1) . Our reasoning to search for additional solutions of the gap equation in the region N = O (1) comes from the T = 0 analysis in 1aper I. The key result of paper I is that for any N < N cr [which is O (1) for a generic γ < n ( ω m ). If system properties evolve smoothly between T = 0and T >
0, each gap function evolves with T and has to vanish at some finite T p,n . By this12ogic, there must be an infinite set of T p,n for any given N < N cr . Because all ∆ n ( ω m ) at T = 0 vanish at N = N cr , all T p,n ( N ) with finite n must approach 0 at N = N cr , as inFig.3(b).We first verify the very assumption that the system behavior evolves smoothly between T = 0 and T >
0. This is not a’priori guaranteed, as at T = 0 the infinite set of solutionsfor all N < N cr emerges due to a fine balance between the tendencies towards NFL andpairing. A finite T perturbs this balance, and it is possible in principle that an infinite setof solutions exists only at T = 0, while at a finite T only a single solution with a finite ∆survives.To address this issue, we recall that the origin for the appearance of the set of ∆ n ( ω m ) at T = 0 is the existence of the solution of the linearized gap equation for all N < N cr ratherthan only for N = N cr . In paper I we obtained the exact solution of the linearized gapequation for all N < N cr and 0 < γ <
1. At small ω m (cid:28) ¯ g , the exact solution reduces to∆( ω m ) = E | ω m | γ/ cos β N log | ω m | ¯ g + φ N ! (23)where E is an infinitesimally small overall magnitude, φ N = O (1) is the phase factor, and β N is the real solution of (cid:15) β N = N , where (cid:15) β = 1 − γ | Γ( γ/ iβ )) | Γ( γ ) πγβ )cos( πγ/ ! (24)A solution of this equation with a real β = ± β N exists for N < N cr , where N cr = 1 − γ ( γ/ γ ) πγ/ ! (25)One can verify that N cr > < γ <
1. At γ (cid:28) N cr ≈ /γ . At γ → N cr → (cid:15) β N = N as the dispersion relation and identify β N with theeffective momentum and N with the effective energy. Then one can define the DoE as ν ( N ) ∝ dβd(cid:15) β (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) β = β N . (26)We plot this function in Fig.4. As expected, it is nonzero for all N < N cr . It is singular near N = 0, where ν ( N ) ∝ (1 /N ) (2 − γ ) / (1 − γ ) and near N cr , where ν ( N ) ∝ / √ N cr − N . This lastsingularity, however, affects ν ( N ) only in the immediate vicinity of N = N cr , as one cansee from Fig. 4. The DoE in (26) is defined up to an overall factor. In Fig.4 we plot ν ( N )without normalizing it. 13 N ν ( N ) N cr γ = 0.5 FIG. 4. The DoE ν ( N ) (in arbitrary units) at T = 0, for γ = 0 .
5. The DoE is non-zero for all
N < N cr . It has a strong singularity near N = 0 and a weaker singularity at N = N cr . Now, for a smooth evolution of system properties between T = 0 and T > ν ( N ) must emerge if we solve the linearized gap equation by approaching the T = 0 limitfrom a finite T . To verify whether this is the case, we keep N = O (1) and set K ∝ /T γ toinfinity in (17).
1. The limit T → Keeping N finite and setting K to infinity, we obtain from (17), after symmetrizationand rescaling ˜Φ( m >
0) = ∞ X n =1 ,n = m ˜Φ( n ) K n,m ( B m B n ) / , (27)where ˜Φ( m ) = Φ( m )( S m /A ( m )) / , B m = S m A ( m ), and S m and K n,m are given by S m = N − A ( m ) m + 1) γ − m γ + 1( m + 1) γ ! , (28) K n,m = 1 | n − m | γ + 1( n + m + 1) γ − m γ + 1( m + 1) γ ! n γ + 1( n + 1) γ ! (29)In explicit form, B m = 2 N ∞ X n =1 n γ − m + 1) γ + m γ + 1( m + 1) γ ! (30)14t large n, m (cid:29) B m ≈ N m − γ / (1 − γ ), K n,m ≈ | n − m | γ + n + m ) γ , ˜Φ( m ) ≈ Φ( m ), andthe equation for the pairing vertex can be approximated by an integral equationΦ( m ) = 1 − γ N Z ∞ dn Φ( n ) n − γ | n − m | γ + 1( n + m ) γ ! (31)This equation has been analyzed in paper I and in earlier works . The solution isΦ( m >
0) = Em γ/ cos ( β N log m + φ ) (32)where β N is the same as in (23) and E is an infinitesimally small overall factor. For these m , Σ m ∝ m − γ , hence ∆( m >
0) = Em γ/ cos ( β N log m + φ ) (33)The functional form of ∆( m ) is the same as for the T = 0 solution, but at this stage φ isa free parameter. To be the solution of the gap equation for all m , this ∆( m ) (or Φ( m ))has to match with the solution at small m , when the discreteness of Matsubara numbersbecomes relevant. If this can be achieved by fixing the value of φ , then ν ( N ) at T → T = 0.This is similar to the T = 0 case, where the oscillating ∆( ω m ) from (23) has to matchwith the nonoscillating solution ∆( ω m ) ∝ / | ω m | γ for | ω m | (cid:29) ¯ g . For T = 0, the exactsolution and the approximate, but highly accurate, analytical solution show that this is thecase. Namely, for some particular φ = φ N , the exact solution has the form of Eq. (23)for | ω m | (cid:28) ¯ g and decays as 1 / | ω m | γ for | ω m | (cid:29) ¯ g . For T →
0, we solve the gap equationnumerically. We present the results for the DoE in Fig.5. While for any finite number M ofMatsubara points in numerical calculations the DoE consists of a discrete set of points, wesee from the histogram of the eigenvalues in panel (a) that at larger M , more eigenvaluesmove to larger N and the number of eigenvalues in any fixed interval of N increases. The“smoothened” ν ( N ) in panel (b) weakly depends on M and is quite similar to the DoE at T = 0 in Fig. 4. This strongly indicates that Eq. (27) has a solution for any N < N cr , likeat T = 0, i.e. the system evolves continuously between T = 0 and T →
2. The case of small but finite T We next consider the case when T is small, but finite. Like we said at the beginningof Sec. III B, at T = 0 and N < N cr , there exists an infinite, discrete set of solutions of15 a) (b) γ = 0.5 γ = 0.5 FIG. 5. Numerical results for the DoE at T →
0. (a) The histogram of the eigenvalues for differentnumbers of sampled Matsubara points M . As M increases, more eigenvalues are shifted towardslarger N , and the number of eigenvalues in any given interval of N increases. (b) the “smoothened”DoE for different M , normalized such that ν ( N = 0 .
5) is kept fixed. As M increases, ν ( N ) doesnot change much. The curves are similar to that in Fig. 4 except for the narrow peak near N cr , asa much larger M is required to get enough eigenvalues near N cr . the nonlinear gap equation, ∆ n ( ω m ). It is natural to expect that each ∆ n smoothly evolveswith T and ends at a finite T p,n . We then expect that Eq. (17) should have an infinite setof solutions at large but finite K . The equation for the pairing vertex in this case has thesame form as in Eq. (27), but one should keep the bare ω term along with the self-energy,i.e., replace A ( m ) by A ( m ) + 2 m + 1 K (34)The implication here is that for large enough m ∼ K /γ , the second term becomes compa-rable to the first one, and this will modify the functional form of ∆( m ) compared to that inEq. (32). The other dependence on K , from K/ ( N − K ) ≈ − N/K , is irrelevant as the1 /K term never becomes comparable to other terms in the rhs of (17).A qualitative argument for the existence of the infinite set of T p,n is the following: • at m (cid:29)
1, the difference between summation over Matsubara numbers n and inte-gration is negligible. Hence, the solution of the gap equation, expressed in terms of ω m ≈ πmT , should be the same as the exact solution at T = 0, Eq. (23). In terms16f m we have ∆( m ) = Em γ/ cos (cid:16) β N log m − β N log K /γ + φ N (cid:17) (35)The phase φ N is fixed by the requirement that ∆( m ) ∝ / | m | γ for m > K /γ , butthere is another phase factor − β N log K /γ , which depends on K . • From the analysis at T → N < N cr .The form of a solution is the same as in (35), ∆( m ) = Em γ/ cos (cid:16) β N log m + ˜ φ N (cid:17) ,with some particular ˜ φ N which depends on N (for a fixed γ ). • We can match the two forms of ∆( m ) by relating the phases:˜ φ N = φ N − β N log K /γ + πn, (36)where n is an integer. Solving (36), we obtain the set of K n , for which this identityholds. Note that the additional factor is πn , not 2 πn , because we can independentlyflip the sign of ∆( m ) at T → m ) at T = 0 intact.Solving Eq. (36) for the critical temperature, we obtain T p,n ∼ ¯ ge − An , A = πβ N γ (37)This is consistent with the result in paper I that at T = 0, ∆ n ( ω m = 0) ∼ ¯ ge − An . As N increases towards N cr , β N ∝ ( N cr − N ) / gets smaller, and all T p,n with n > N cr − N : T p,n ∼ ¯ ge − bn/ ( N cr − N ) / , (38)where b = O (1). Eq. (38) shows that all T p,n ( N ) terminate simultaneously at N = N cr .A remark is in order here. Eq. (37) and the T = 0 result for ∆( m ) are based onthe assumption that the solution of the linearized gap equation oscillates up to ω m ∼ ¯ g and decays as 1 / | ω m | γ at larger ω m . This holds for most of γ <
1, but for very small γ ,oscillations extend to larger scale. In this case, Eq. (36) has to be modified. We discuss thisin Appendix B.In Fig. 6 we show the numerical results for T p,n for n = 1 −
4, along with the result for T p, from Eq. (21). We set representative γ to be γ = 0 . γ = 0 .
5. We verified that17 g / T ( n ) c ) γ = 0.3 g / T ( n ) c ) N γ = 0.5 N γ = 0.3 γ = 0.51/log( ¯ g / T p , n ) 1/log( ¯ g / T p , n ) n=0n=1n=2n=3n=4 n=0n=1n=2n=3n=4 FIG. 6. The numerical solution of the gap equation for small but finite T for γ = 0 . .
5. Thetemperatures T p,n are the onset temperatures for the pairing in different topological sectors (thecorresponding eigenfunctions change sign n times as functions of discrete Matsubara frequency).The highest T p, ∝ /N /γ terminates at N = ∞ . Analytical reasoning shows that all other T p,n vanish at N = N cr (big red dot). Numerical results show that T p,n with finite n > N = N cr . T p, behaves as 1 /N /γ , as expected. At other T p,n lines exponentially approach zero as N tends to N cr , and the slope becomes larger as n increases. To obtain this behavior we useda “hybrid frequency scale” method, which allowed us to numerically cover an exponentiallylarge frequency range and reach very low T , while keeping track of the Matsubara summationat the lowest frequencies. This is achieved by adopting a frequency mesh that overlaps withMatsubara frequencies ω m = πT, πT, · · · at small values and crosses over to a logarithmicalspacing beyond a certain scale, above which the discreteness of the Matsubara sum becomesunimportant. We discuss this method in Appendix C. It was also used in Ref. 9 in addressingthe interplay between the first-Matsubara physics and thermal fluctuations.In Fig. 7 we plot T p,n with n up to 17 for particular N = 1. We clearly see that T p,n scale as e − An , as in Eq. (37). We extracted β N =1 from the fit to the exponential form andobtained β N =1 = 1 .
62 for γ = 0 . β N =1 = 1 .
12 for γ = 0 .
5. These values are quite closeto the exact values, extracted from Eq. (24): β N =1 = 1 .
71 for γ = 0 . β N =1 = 1 .
27 for γ = 0 .
5. The small difference comes from the numerical error of the hybrid frequency scale,which effectively shifts N up by roughly 0 .
07 for γ = 0 . .
13 for γ = 0 . n -40-30-20-10 l og T =0.3, N=1 numerical resultslinear fittinglog T=-2.805n+0.4508 n -40-35-30-25-20-15-10 l og T =0.5, N=1 numerical resultslinear fittinglog T=-2.453n+0.5721 log T p , n = − 2.805 n + 0.4508 log T p , n = − 2.453 n + 0.5721 FIG. 7. T p,n as functions of n for N = 1 and γ = 0 . .
5. The black dots are the data, obtainedby varying temperature to find the eigenvalue N = 1, and the red lines are linear fitting of log T p,n ,with fitting parameters given in the boxes. We clearly see that T p,n scales as e − An , as we obtainedanalytically. C. The structure of the gap function, ∆ n ( m ) Another result of the T = 0 analysis is that the solutions of the nonlinear gap equationwith different n are topologically distinct — the gap function ∆ n ( ω m ) changes sign n timesas a function of Matsubara frequency. Because we expect that ∆ n ( m ), which develops below T p,n , becomes ∆ n ( ω m ) at T = 0, it should change sign n times as a function of Matsubaranumber m . The same should be true for the pairing vertex Φ n ( m ). In Fig.8 we show Φ n ( m )for a few smallest n and for n = 16 ,
17. We see that Φ n ( m ) indeed changes sign n times. Atlarge n , Φ n ( m ) oscillates at large m as a function of log m , with the amplitude proportional to1 / | m | γ/ . This is exactly the same behavior as in Eq. (35), given that ∆ n ( m ) ∼ Φ n ( m ) | m | γ .For comparison, in the last panel of Fig. 8 we plot the exact Φ( ω m ) at T = 0. We see thatat T = T p,n , the form of Φ n ( m ) for m (cid:29) T = 0.That Φ n ( m ) has to change sign at least once follows from the relation between Φ(0) andΦ( m > K (cid:29) N :Φ(0) ≈ − K ∞ X m =1 Φ( m ) A ( m ) + m +1 K m γ + 1( m + 1) γ ! (39)This relation shows that even if Φ( m ) has the same sign for all m >
0, Φ(0) would still beof opposite sign. This is consistent with Fig. 8, which shows that Φ ( m ) changes sign at19 -40 -30 -20 -10 -40 -30 -20 -10 -40 -30 -20 -10 -40 -30 -20 -10 -11 -40 -30 -20 -10 -11 ω γ / m Φ ( ) m ω γ / m Φ ( ) m ω γ / m Φ ( ) m ω γ / m Φ ( ) m ω γ / m Φ ( ) m ω m / ¯ g T p ,17 / ¯ g = 5 × 10 −42 T p ,16 / ¯ g = 1 × 10 −39 T p ,0 / ¯ g = 0.8 T p ,1 / ¯ g = 3.4 × 10 −3 T p ,2 / ¯ g = 0.2 × 10 −5 ω / ¯ g ω γ / Φ ( ω ) (a)(b)(c)(d)(e) (f) ω = πT p ,0 ω = πT p ,1 ω = πT p ,2 ω = πT p ,16 ω = πT p ,17 n=1, one sign changen=2, two sign changes n=16n=17n=0, sign preserving ω γ / m Φ ( m ) ω γ / m Φ ( m ) ω γ / m Φ ( m ) ω γ / m Φ ( m ) ω γ / m Φ ( m ) ω γ / m Φ ( ω m ) ω m / ¯ g FIG. 8. Panels (a)-(e) – The pairing vertex, Φ n ( m ) at T = T p,n , as a function of the Matsubarafrequency ω m = πT (2 m + 1) for representative parameters γ = 0 . N = 1. We show Φ n ( m )for n = 0 , , , ,
17. The corresponding T p,n are shown in the figures. For n = 0, Φ ( m ) ∝ (1 | m | γ + 1 / | m + 1 | γ ) does not change sign. Other Φ n ( m ) change sign n times, and T p,n ∝ e − An .The results for n = 16 ,
17 show that Φ n ( m ) oscillates at m (cid:29) m , with theamplitude proportional to 1 / | m | γ/ . Panel(f) - Φ( ω m ) from the exact solution of the linearizedequation for the pairing vertex at T = 0. The positions of zeros of Φ n ( m ) are marked by crosses.The smallest frequencies ω = πT p,n (different for different n ) are shown by arrows. -21 nonlinear eq. at 0.9T c(0) linear eq. at T c(0) -22 nonlinear eq. at 0.9T c(1) linear eq. at T c(1) -22 nonlinear eq. at 0.9T c(2) linear eq. at T c(2) ω γ / m Φ ( ) m ω γ / m Φ ( ) m ω γ / m Φ ( ) m ω m / g T p ,0 T p ,0 T p ,1 T p ,1 T p ,2 T p ,2 n = 0 n = 1 n = 2 FIG. 9. The solution of the nonlinear equation for the pairing vertex for T = 0 . T p,n , along withthe solution at T = T p,n −
0. The three panels are the solutions for n = 0 , ,
2. The number of signchanges remains the same at T p,n and 0 . T p,n , as indicated by the blue arrows, and the frequencies,at which Φ n ( m ) changes sign, do not shift with T . m = O (1) and keeps the same sign at larger m . The same holds for larger n - the first signchange occurs at m = O (1), i.e., at ω m ∼ T p,n . In Fig. 9 we show the results for Φ n ( m ) for n = 0 , , T ≤ T p,n .We expanded to order Φ m and used the solution at T = T p,n as the source. We see that thenumber of sign changes remains the same, and the frequencies, at which the sign of Φ n ( m )changes, remain essentially independent on T . This is consistent with the result in paper Ithat at T = 0, ∆ n ( ω m ) changes sign n times at finite ω m .Finally, in Fig.10 we follow ∆ n ( m ) along the line T p,n ( N ) and show the evolution of thefrequency ω max , at which ∆ n ( m ) changes sign for the last time (i.e., all sign changes are at ω m < ω max ). Because all T p,n terminate at N = N cr , ω max must shrink, as T p,n decreases,and vanish at T p,n →
0, because right at N = N cr and T = 0, ∆( ω m ) is sign-preserving (seepaper I). The data in Fig.10 show that ω max indeed decreases with decreasing T p,n for allvalues of n , as long as n remains finite. 21 a) (b) n=15 n=13 n=11n=9 n=15n=13 n=11n=9n=7 n=7 T p , n / ¯ g T p , n / ¯ g FIG. 10. The highest frequency ω max , at which ∆ n ( m ) changes sign, plotted vs T p,n for various n .The arrows indicate the direction towards smaller T p,n . We see that ω max decreases together with T p,n , i.e., n sign changes of ∆ n ( m ) occur at progressively smaller Matsubara frequencies. IV. AWAY FROM A QCP
We now analyze how T p,n and the DoE at T → M b . We argue that a finite M b introduces qualitativechanges in the system behavior, i.e., there is a qualitative difference between the structure ofthe DoE at a finite M b and at M b = 0. Specifically, we argue that a finite M b (i) makes thenumber of T p,n for a given N finite and (ii) splits T p,n at the smallest T such that different T p,n terminate at different N cr,n , all smaller than N cr . The temperature T p, is still nonzerofor any N , but at a finite M b it acquires a conventional form T p, ∼ M b e − /λ , where λ ∝ /N [see Eq. (43) below].We begin with the analytical analysis. We model the bosonic propagator away from aQCP by V (Ω m ) = ¯ g γ (Ω m + M b ) γ/ (40)The linearized gap equation becomes∆( ω m ) = ¯ g γ N πT X m = m ∆( ω m ) − N ∆( ω m ) ω m ω m | ω m | | ω m − ω m | + M b ) γ/ . (41)Consider first the limit T →
0. Replacing the summation over frequency by integration with T as the lower limit, we immediately find that at a finite M b there is a simple difference22etween sign-preserving and sign-changing solutions for ∆( ω m ). For the sign-preservingsolution, ∆(0) is finite. Taking properly the limit ω m →
0, we obtain T p, from the self-consistent equation on ∆(0): (cid:18) (cid:18) ¯ gM b (cid:19) γ (cid:19) = 1 N (cid:18) ¯ gM b (cid:19) γ log M b T p, . (42)The second term in the lhs is the contribution from the self-energy, which away from a QCPhas a Fermi-liquid form at frequencies below M b . Solving for T p, , we find T p, ∼ M b e − N (1+ λ ) /λ , λ = 1 + (cid:18) ¯ gM b (cid:19) γ (43)We see that T p, is still nonzero for any N , however its dependence on N is exponential.This is similar to the case of a BCS superconductor, where T c is finite for arbitrary weakcoupling, albeit exponentially small. For N = 1, Eq. (43) has the same structure asMcMillan formula T c ∼ ω D e − (1+ λ ) /λ (Ref. ). In qualitative distinction to the behavior at aQCP, now the existence of a nonzero T p, for arbitrary large N is due to ordinary Cooperlogarithm in a Fermi-liquid rather than to special properties of fermions with ω m = ± πT ina NFL regime. As a consequence, ∆ ( ω m ), emerging below T p, , does not vanish at T = 0,i.e, at any N the ground state is a superconductor. In this respect, there is no critical N cr ,separating normal and superconducting states at T = 0.For solutions with n >
0, ∆ n ( ω m ) must vanish at ω = 0 because there is just a singlesolution with a finite ∆(0). This sets the condition∆(0) (cid:18) (cid:18) ¯ gM b (cid:19) γ (cid:19) = ¯ g γ N Z ∞ ∆ n ( ω m ) ω m dω m (( ω m ) + M b ) γ/ = 0 (44)We show below that ∆ n ( ω m ) scale as ( ω m ) at small ω m , hence the integral in (44) isinfra-red convergent. Eq. (44) then implies that ∆ n ( ω m ) must change sign n times atsome finite ω m . This is qualitatively different from the situation at a QCP. There, all T p,n with finite n terminate at T = 0 at N = N cr . The gap function at N = N cr vanishes at ω m = 0, yet it remains sign-preserving (see paper I for details). This holds at a QCP because∆( ω m ) ∝ | ω m | γ/ , in which case ∆( ω m ) /ω m is singular at ω m →
0, and one cannot just set ω m = 0 in the gap equation, as it is done in Eq. (44). At small M b , T p,n approaches zero at N only slightly below N cr , and ∆ n ( ω m ) must recover the gap function at a QCP at N = N cr at frequencies above some scale, which vanishes when M b →
0. Because ∆( ω m ) at N = N cr is sign-preserving, the n sign changes of ∆ n ( ω m ) have to occur below this scale.23 =0n=1n=2n=3n=4 n=0n=1n=3n=4n=0n=1n=3n=4 n=0n=2n=3n=4 T p , n / ¯ g T p , n / ¯ gT p , n / ¯ g T p , n / ¯ g NNN N (a) M b / ¯ g = 0 (b) M b / ¯ g = 10 −4 (c) M b / ¯ g = 10 −3 (d) M b / ¯ g = 10 −1 n=2 n=2 n=1 FIG. 11. The solutions of the linearized gap equation for a finite boson mass M b . Different panelsare for different M b / ¯ g , shown in the figures. We set γ = 0 .
5. The critical temperatures T p,n nowterminate at different N cr,n . This is qualitatively different from the behavior at a QCP, where all T p,n with n > N = N cr . We now expand in ω m the rhs of the gap equation (41) for ∆ n> ( ω m ). Expanding andusing (44) to cancel out the leading term, we obtain ∆ n ( ω m ) = A n ( ω m /M b ) , where A n isgiven by A n = (cid:18) ¯ gM b (cid:19) γ γN Z ∞ dx ∆ n ( x ) x x (1 + γ ) − x + 1) γ/ (45)where x = ω m /M b . The integral in (45) converges at x = O (1), hence by order of magnitude A n ∼ ¯ g γ / ( M γb N ) with n − dependent prefactor. For an estimate, we assume that at large n the integral is determined by x n ∼ /n , before oscillations begin, and that ∆ n ( x 5. We see thatthe histogram is heavily shifted towards N = 0 because now there is a finite number of points in agiven interval around a particular N even when the total number of sampling Matsubara points M tends to infinity. This is qualitatively different from the case M b = 0 in Fig. 5, where the numberof points in a given interval around any N < N cr scales with M . M b . We see that T p,n indeed terminate at different N . We verified that T p, is exponentialin N , like in Eq. (43). The number of solutions, for which T p,n crosses N = 1, is finite forany nonzero M b . It decreases one by one with increasing M b and vanishes once M b exceedssome critical value. At larger M b , there is only one onset temperature T p, for the physicalcase N = 1, and the behavior below T p, is qualitatively the same as in BCS theory.The same behavior shows up in the analysis of the DoE at T → 0. Because there isonly a finite number of termination points of T p,n in any finite interval of N , the normalizedDoE ν ( N ) vanishes for all N = 0 in the formal limit M → ∞ , where M is the numberof Matsubara points, probed in a numerical calculation. Because termination lines clusteraround N = 0, and the total R dN ν ( N ) = 1, the DoE becomes ν ( N ) = δ ( N ) in this limit.In practice, this implies that the DoE will decrease for any finite N when M increases, andthe histogram of ν ( N ) will shift with increasing M towards N = 0. In Fig... we show thenumerical results for ν ( N ). We see precisely this behavior. Note also that the largest N cr, shifts down from N cr as M b increases. 25 . CONCLUSIONS In this paper, the second in a series, we continued our analysis of the interplay betweenpairing and NFL in the model of fermions with singular dynamical interaction V (Ω m ) ∝ / | Ω m | γ (the γ model). We introduced a knob (the parameter N ) to vary the relativestrength of the dynamical interaction in the particle-hole and particle-particle channels. Inpaper I we studied the γ model at T = 0 and 0 < γ < 1. We found that superconductingorder exists for N < N cr ( γ ) and that for any such N there is an infinite number of solutionsof the nonlinear gap equation, ∆ n ( ω m ), specified by the number of times ( n ) the gap functionchanges sign as a function of ω m . At large n , ∆ n (0) = e − an is exponentially small in n .In this paper, we analyzed the same model at a finite T . We showed the DoE remainscontinuous at T → T = 0, i.e., the existence of an infinite setof solutions of the gap equation is not a special property of T = 0. We further showed thateach ∆ n ( ω m ) vanishes at a critical temperature T p,n ∼ ∆ n (0). Viewed as functions of N , all T p,n ( N ) with n > N = N cr , which turns out to be a multicriticalpoint (see Fig.3(b) and Fig.6). The eigen-function ∆ n ( m ) at T = T p,n − n times as a function of Matsubara number m , and retains this property down to T = 0,where it becomes a continuous function ∆ n ( ω m ). The temperature T p, behaves differently– it remains finite for all N and scales at large N as T p, ∼ /N /γ . In this limit, the pairingat T p, predominantly involves fermions with the two lowest Matsubara frequencies ± πT ,because for these fermions nonthermal self-energy vanishes, which implies that pairing doesnot compete with NFL. At N > N cr , ∆ ( m ) initially grows below T p, , but then changestrend and vanishes at T → T p,n with different n terminate atdifferent N cr,n , and for any N > 0, the number of T p,n is finite. It drops one by one as thebosonic mass M b gets larger, and disappears above a certain M b . The onset temperature T p, is still finite for any N , but the gap ∆ ( m ) does not vanish at T → 0, i.e., the groundstate is a superconductor for any value of N .The next paper will be devoted to analysis of how the physics of the γ model at a QCPchanges between γ < γ > 1. 26 CKNOWLEDGMENTS We thank I. Aleiner, B. Altshuler, E. Berg, D. Chowdhury, L. Classen, R. Combescot,K. Efetov, R. Fernandes, A. Finkelstein, E. Fradkin, A. Georges, S. Hartnol, S. Karchu,S. Kivelson, I. Klebanov, A. Klein, R. Laughlin, S-S. Lee, G. Lonzarich, D. Maslov, F.Marsiglio, M. Metlitski, W. Metzner, A. Millis, D. Mozyrsky, C. Pepin, V. Pokrovsky, N.Prokofiev, S. Raghu, S. Sachdev, T. Senthil, D. Scalapino, Y. Schattner, J. Schmalian, D.Son, G. Tarnopolsky, A-M Tremblay, A. Tsvelik, G. Torroba, E. Yuzbashyan, and J. Zaanenfor useful discussions. The work by AVC was supported by the NSF DMR-1834856. Appendix A: Free energy below T p, at large N In this Appendix we address one issue about the system behavior at large N , below T p, .As we said in the text and in earlier papers , T p, remains finite for arbitrary large N because the pairing in the large N limit predominantly involves fermions with Matsubarafrequencies ω m = ± πT (Matsubara numbers m = 0 , − T p, involves fermions with ω m = ± πT , whichdevelop ∆ ( m = 0) = ∆ ( − 1) below T p, . A much smaller gap ∆ ( m ) at other Matsubarafrequencies is then induced by proximity. The gap ∆ (0) initially increases as T decreases,but never exceeds T . At the smallest T , ∆ (0) decreases proportional to T and vanishes at T = 0, where a Matsubara frequency becomes a continuous function, and the special role offrequencies ± πT is lost.The very existence of the solution of the Eliashberg equations with a nonzero ∆ ( m )implies that the free energy of such state, F p, , is smaller than that of the normal state (seePaper I and Refs. ). In explicit form, in this case δF = F p, − F n = − ¯ g ν N /γ t − γ ( t γ − (A1)where t = T /T p, and ν is the density of states at the Fermi level in the normal state.The issue we address here is whether a negative δF = δE − T δS comes from the changeof the energy, like in BCS superconductor, or from the change of the entropy. In our case,27 .0 0.2 0.4 0.6 0.8 1.0 - - - × - - × - - × - - × - - × - - × - - × - - × - - × - - - T / T (0) c δF = F sc − F n δS = S sc − S n δE = E sc − E n in units ofin units ofin units of g N gN g N ¯ g ν ¯ g ν ¯ gν T / T p ,0 FIG. 13. The difference of free energy, entropy and energy between the SC state and normal state,as a function of temperature. We choose γ = 0 . N = 5 > N cr . we can compute both terms explicitly at large N . We obtain δS = S p, − S n = ¯ gν πN /γ t − γ ( t γ − 1) ( t γ ( γ + 2) + γ − δE = E p, − E n = ¯ g ν N /γ t − γ ( t γ − t γ ( γ + 1) + γ − 1) (A2)We plot δF, δS , and δE in Fig.13 for a particular case of γ = 0 . N = 5 N cr ( γ = 0 . 5) =4 . < 5. We clearly see that all three quantities vanish at T = 0, as expected. We also seethat there is a clear change of behavior at larger t , when ∆ ( m ) increases with decreasing t , and at smaller t , where ∆ ( m ) scales with t . In the first case, δE is negative and δS ispositive, i.e., negative δF is due to negative δE , as in a BCS superconductor. However, atsmaller t , δE is reduced and actually becomes positive at the smallest t . In this t range,a negative δF is entirely due to positive T δS , i.e smaller free energy for the state with afinite ∆ ( m ) is an entropic effect. This is also consistent with the fact that the gap functionvanishes at T = 0. Appendix B: Onset temperatures T p,n at small γ The limit γ → γ model attracted a lot of attention from various sub-communitiesin physics and has been analyzed in both Eliashberg-type and renormalization28roup approaches. In this Appendix, we present two expressions for T p,n at small γ and N = O (1), which is much smaller than N cr = 4 /γ . One expression is for the pure γ modelwith V (Ω) = (¯ g/ | Ω | ) γ . Another is for the case when we introduce an upper cutoff for V (Ω m )at Ω m = Λ, i.e., modify the interaction to V (Ω m ) = ¯ g | Ω m | ! γ " − | Ω m | Λ ! γ (B1)for | Ω m | < Λ, and V (Ω m ) = 0 for | Ω m | > Λ. The distinction is important for small γ becauseat γ → V (Ω) ∝ | Ω | − γ becomes frequency independent, and the pairing problem withinthe γ model becomes equivalent to BCS, but without the upper cutoff for the interaction.In this situation, the cutoff at Λ is necessary to avoid ultra-violet singularities. At the sametime, at a finite γ , the gap equation is free from ultra-violet singularities already without Λ,and if T p,n are smaller than Λ, the cutoff can be safely neglected.We consider the cases when Λ is set to infinity (the pure γ model) and when Λ is finiteseparately, and then obtain an interpolation formula linking the two cases. T p,n in the pure γ -model The point of departure for the calculation of T p,n at small γ is the observation in PaperI and in earlier works that for N = O (1), T p, is determined by fermions withfrequencies much larger than ¯ g . For these fermions, the normal state self-energy at T = 0,Σ( ω m ) ≈ | ω m | − γ ¯ g γ sgn ω m is smaller than the bare ω m and can be safely neglected. Withoutthe self-energy, Φ( m ) = ∆( m ) and the gap equation reduces to∆( m ) = KN X m = m ∆( m ) | m + 1 | | m − m | γ . (B2)where, we remind the reader, K = (¯ g/ (2 πT )) γ . For m, m (cid:29) 1, the summation over m canbe replaced by the integration. Also, for small γ , 1 / | m − m | γ is reasonably well approximatedby 1 / | m | γ for | m | > | m | and by 1 / | m | for | m | > | m | . As a result, the gap equation for m (cid:29) m ) = KN | m | γ Z m dm ∆( m ) m + Z ∞ m dm ∆( m )( m ) γ ! (B3)Introducing z = m γ and differentiating twice over z , we obtain the differential equation( z ∆( z )) = − KN γ ∆( z ) z (B4)29he solution of this equation must satisfy two boundary conditions. One follows from B2:∆( z → ∞ ) ∝ /z . Another comes from the boundary at m = O (1), i.e., at z ≈ 1, wherethe discretization in Eq. B2 becomes relevant. Sign-preserving ∆ ( m ) is flat at m = O (1)and deviates from ∆ (0) only by O ( γ ). Then, to leading order in γ , the second boundarycondition is ∆ ( z = 1) = 0. For sign-changing ∆ n ( m ), we can use the fact that the gapfunction changes sign at m = O (1) and use ∆( z = 1) = 0.The generic solution of (B4) is∆( z ) = 1 √ z " A J s KN γz ! + A Y s KN γz ! , (B5)where J and Y are Bessel and Neumann functions.At large z , J ∝ / √ z and Y ∝ √ z . To satisfy the boundary condition at z → ∞ wemust set A = 0. Then ∆( z ) = A √ z J s KN γz ! . (B6)A similar result has been obtained in Ref. . Using the second boundary condition for sign-preserving solution ( n = 0 in our nomenclature) and using ( xJ ( x )) = xJ ( x ), we obtain J s KN γ ! = 0 , (B7)where, we remind the reader, K = (¯ g/ (2 πT )) γ . The highest T at which this equation hasthe solution is T = T p, , where T p, = c ¯ g (1 . N γ ) − /γ (B8)where c = O (1). For N = 1 this agrees with Ref. (see also earlier results ). We see that T p, is indeed much larger than ¯ g and tends to infinity when γ → 0. As we said, this is theconsequence of the fact that at γ → γ model reduces to the BCS model without thecutoff on the interaction.The prefactor c has been computed in Ref. , and the result is c = 0 . T =0, ∆ ( ω m = 0) has been obtained in Paper I and in Ref. : ∆ (0) = 0 . g (1 . N γ ) − /γ .The ratio 2∆ /T c = 3 . 53, as in BCS theory.For sign-changing solutions, the boundary condition J (2 q K/N γ ) = 0 yields T p,n = c n ¯ g ( a n N γ ) − /γ (B9)30here c = O (1) and a = 3 . a = 12 . n , forwhich T p,n = O (¯ g ). For larger n , T p,n ∝ e − An , like in Eq. (37). The ratio 2∆ n (0) /T p,n is oforder 1 for n = O (1), but we did not compute it explicitly.Note that although T p,n with n = O (1) are much larger than ¯ g , the ratio T p, /T p, ∼ (0 . /γ is parametrically small and vanishes when γ → T p,n in the γ model with the infra-red cutoff We now consider the opposite limit when γ → 0, but the upper cutoff for V (Ω m ) is afinite Λ. In this case, we use Eq. (B1) for the interaction.The derivation of the differential equation proceeds along the same lines as before, andthe result is the same Eq. (B4) as before. However, now the solution, ∆( z ), has to satisfythe boundary condition ∆( K ∗ ) = 0, where K ∗ = (Λ / πT ) γ . Such a solution is∆( z ) ∝ √ z × " J s KN γz ! Y s KN γK ∗ ! − Y s KN γz ! J s KN γK ∗ ! . (B10)In the limit K/ ( N γK ∗ ) = (¯ g/ Λ) γ / ( N γ ) (cid:29) 1, valid when γ → J ( x ) ≈ s πx cos ( x − π/ ,Y ( x ) ≈ s πx sin ( x − π/ , (B11)and obtain ∆( z ) ∝ z / sin s KN γK ∗ − s KN γz ! (B12)To the leading order in γ , this reduces to∆( m ) ∝ sin r γN log 2 πT | m | Λ ! (B13)Using the second boundary condition for sign-preserving solution, we find T p, ∼ Λ e − π/ √ λ γ , λ γ = γN (B14)This agrees with Refs. . For n > 0, using ∆( z = 1) = 0, we obtain T p,n ∼ Λ e − nπ/ √ λ γ , n = 1 , ... (B15)31gain, T p, is parametrically smaller than T p, , the ratio T p, /T p, vanishes at γ → T p,n with larger n are even smaller. Note that T p,n is independent on γ .The zero-temperature gap ∆ n ( ω m = 0) in the same limit has been computed in paper I.Comparing it with T p,n , Eqs. (B14) and (B15), we find that T p,n and ∆ n (0) are of the sameorder. A higher accuracy of calculations is required to compute 2∆ /T p ratios. 3. Intermediate regime We now argue that Eq. (B10), describes the gap functions in both, the γ model withthe cutoff at Λ, and the pure γ model. Indeed, in the opposite limit of large Λ and finite γ , K/ ( N γK ∗ ) = (¯ g/ Λ) γ / ( N γ ) (cid:28) 1. In this case Y (cid:16) q KNγK ∗ (cid:17) (cid:29) J (cid:16) q KNγK ∗ (cid:17) . Then only Y (cid:16) q KNγK ∗ (cid:17) should be kept in (B10), and we recover Eq. (B6). In this limit, the functionalform of ∆( z ) does not depend on Λ and we recover the result for T p,n for the pure γ model.The crossover between the two regimes occurs at (¯ g/ Λ) γ / ( N γ ) = O (1). Appendix C: The hybrid frequency scale Within our computation capacity, the maximum size of frequency mesh we can take,is approximately 10 . If we directly sum over Matsubara frequencies ω m = (2 m + 1) πT ,the maximum frequency we can reach is around 10 T . Because we are interested in fre-quencies ω m < ¯ g , the lowest temperature one can reach is T ∼ − ¯ g . To access expo-nentially lower temperatures, one can in principle use a logarithmical frequency mesh inwhich log( ω m /πT ) ∝ m , but in this scheme one loses the special role of the first Matsubarafrequency and, more generally, of Matsubara frequencies with m = O (1), which, e.g., set thephase in the expression for ∆( ω m ) at T → 0, Eq. (33) To overcome this difficulty we adopta hybrid frequency scaling. Namely, we set ω m = (2 m + 1) πT, m < m L (2 m + 1) πT + e km − b πT, m ≥ m L (C1)Here n L is some number, below which we adopt the original formula for a Matsubara fre-quency, and beyond which we add an exponentially growing term. In practice, we havetaken m L ∼ . M , where M is the total number of frequency points. One should also32roperly choose the parameters k and b such that when m → m L , this exponential term canbe neglected compared to πT and when m → M , the exponential term dominates over thelinear term and can reach our upper limit of frequency. The change of the frequency formalso induces a corresponding change to the measure of summation through a Jacobian, i.e.when the hybrid frequency is used ( m ≥ m L ), the following adjustment should be applied, πT X ... → πT X k e km − b ! ... (C2)On the other hand, the Matsubara summation for the self-energy is well approximatedusing the Euler-Maclaurin formula. In the normal state and at zero bosonic mass M b = 0,Σ( ω m ) = ¯ g γ (2 πT ) − γ H m,γ , Eq. (7), where H m,γ is a generalized harmonic number. Al-though H m,γ is a tabulated function in some computation libraries, we note that it is wellapproximated by H m,γ = m X p =1 p γ ≈ − γ (cid:16) m − γ − (cid:17) + 12 (cid:18) m γ (cid:19) + γ (cid:18) − m γ +1 (cid:19) − γ ( γ + 1)( γ + 2)720 (cid:18) − m γ +3 (cid:19) + γ ( γ + 1)( γ + 2)( γ + 3)( γ + 4)30240 (cid:18) − m γ +5 (cid:19) (C3)In Fig. 14 we compare our numerical Σ( ω m ) at about the smallest T that we used, to thezero-temperature expression Σ( ω m ) = g γ − γ ω − γm . We see that the numerical and the T = 0expressions expectedly coincide at Matsubara numbers m (cid:29) 1, but differ at m = O (1), andthe difference increases as γ gets larger.In the case when the bosonic mass M b is finite, the self-energy in the normal state is givenby Σ( ω m ) = ¯ g γ (2 πT ) − γ H m,γ ( δ ), where δ = M b / (2 πT ) and H m,γ ( δ ) = P mp =1 1( p + δ ) γ/ . This H m,γ ( δ ) is not a tabulated function, but based on what we know about the case M b = 0, we33 -50 n /g -5 ( n ) / g =0.9 finite T resultT=0 result -50 n /g -10 ( n ) / g =0.7 finite T resultT=0 result -50 n /g -40 -20 ( n ) / g =0.5 finite T resultT=0 result -50 n /g -40 -20 ( n ) / g =0.3 finite T resultT=0 result FIG. 14. Comparison between self energy obtained though finite temperature( T /g = 10 − ) anal-ysis and zero-temperature analysis. In all subplots the red curves represent results obtained viahybrid scaling frequency technique and though the Euler-Maclaurin formula, while the dark bluecurves represent T = 0 analysis Σ( ω ) = g γ − γ ω − γ . expect the Euler-Maclaurin formula to work well. Using this formula, we obtain H m,γ ( δ ) = m X p =1 p + δ ) γ/ ≈ (cid:18) δ (cid:19) γ m F [ 12 , γ , , − m δ ] − F [ 12 , γ , , − δ ] ! + 12 δ ) γ/ + 1( m + δ ) γ/ ! + γ − m ( m + δ ) γ/ + 1(1 + δ ) γ/ ! − γ ( γ + 2) − m ( γ + 4)( m + δ ) γ/ + 3 m ( m + δ ) γ/ + γ + 4(1 + δ ) γ/ − δ ) γ/ ! + 130240 γ ( γ + 2)( γ + 4) − m ( γ + 6)( γ + 8)( m + δ ) γ/ + 10 m ( γ + 6)( m + δ ) γ/ − m ( m + δ ) γ/ + ( γ + 6)( γ + 8)(1 + δ ) γ/ − γ + 6)(1 + δ ) γ/ + 15(1 + δ ) γ/ ! (C4)We used this form to compute Σ( ω m ) for a finite M b .34ll numerical results, reported in this paper, are obtained using Matlab2017a with thedefault double precision and Wolfram Mathematica 11.1. A. Abanov and A. V. Chubukov, arXiv:2004.13220 (2020). This issue will be discussed in more detail for the case γ > 1, where additional vortices appearin the complex plane for all ∆ n ( z ). Topological aspects of the existence of zeros of the gap function in momentum space have beenextensively discussed in the literature (see e.g., Refs. ). Our case is a dynamical versionof a nodal topological superconductor. The classification and consequences of dynamical nodaltopology in frequency space remain an open question. Y. Wang, A. Abanov, B. L. Altshuler, E. A. Yuzbashyan, and A. V. Chubukov, Phys. Rev.Lett. , 157001 (2016). Y.-M. Wu, A. Abanov, Y. Wang, and A. V. Chubukov, Phys. Rev. B , 144512 (2019);A. Abanov, Y.-M. Wu, Y. Wang, and A. V. Chubukov, Phys. Rev. B , 180506 (2019); A. V.Chubukov, A. Abanov, Y. Wang, and Y.-M. Wu, Annals of Physics , 168142 (2020). Although ∆ n ( ω m ) at T > n still exists because, e.g., on a realaxis the total change of the phase of a complex ∆ n ( ω ) = | ∆ n ( ω ) | e iψ n ( ω ) between ω = −∞ and ω = ∞ has extra 2 πn compared to the case n = 0. A. A. Abrikosov, L. P. Gorkov, and I. E. Dzyaloshinski,Methods of Quantum Feld Theory in Statistical Physics (Pergamon Press, Oxford, 1965). S. Raghu, G. Torroba, and H. Wang, Phys. Rev. B , 205104 (2015). H. Wang, Y. Wang, and G. Torroba, Phys. Rev. B , 054502 (2018). A. A. Abrikosov and L. P. Gorkov, Soviet Phys. JETP , 1243 (1961). H. Wang, S. Raghu, and G. Torroba, Phys. Rev. B , 165137 (2017). A. Abanov, A. V. Chubukov, and A. M. Finkel’stein, EPL (Europhysics Letters) , 488 (2001);A. Abanov, A. V. Chubukov, and J. Schmalian, Advances in Physics , 119 (2003); E.-G.Moon and A. Chubukov, Journal of Low Temperature Physics , 263 (2010). W. L. McMillan, Phys. Rev. , 331 (1968). G. M. Eliashberg, JETP , 696 (1960). Y. Wada, Phys. Rev. , A1481 (1964); D. J. Scalapino, Y. Wada, and J. C. Swihart,Phys. Rev. Lett. , 102 (1965); J. Bardeen and M. Stephen, Phys. Rev. , A1485 (1964);R. Haslinger and A. V. Chubukov, Phys. Rev. B , 214508 (2003). E. A. Yuzbashyan, A. V. Chubukov, A. Abanov, and B. L. Altshuler, In preparation. Y.-M. Wu, A. Abanov, and A. V. Chubukov, Phys. Rev. B , 014502 (2019). A. L. Fitzpatrick, S. Kachru, J. Kaplan, S. Raghu, G. Torroba, and H. Wang, Phys. Rev. B , 045118 (2015). M. A. Metlitski, D. F. Mross, S. Sachdev, and T. Senthil, Phys. Rev. B , 115111 (2015). D. F. Mross, J. McGreevy, H. Liu, and T. Senthil, Phys. Rev. B , 045121 (2010). D. T. Son, Phys. Rev. D , 094019 (1999); A. V. Chubukov and J. Schmalian, Phys. Rev. B , 174520 (2005). M. Sato and Y. Ando, Reports on Progress in Physics , 076501 (2017). R.-X. Zhang, W. S. Cole, X. Wu, and S. Das Sarma, Phys. Rev. Lett. , 167001 (2019). S. Matsuura, P.-Y. Chang, A. P. Schnyder, and S. Ryu, New Journal of Physics , 065001(2013)., 065001(2013).