The special role of the first Matsubara frequency for superconductivity near a quantum-critical point -- the non-linear gap equation below $T_c$ and spectral properties in real frequencies

TL;DR

This paper investigates how the special role of the first Matsubara frequency influences superconductivity near a quantum-critical point, revealing unique spectral properties and gap behaviors below the critical temperature.

Contribution

It provides a detailed analysis of the non-linear gap equations in real frequencies, highlighting the distinct spectral features when the first Matsubara frequency is special, and compares these with experimental data in cuprates.

Findings

Spectral functions show gap filling rather than closing at $T_c$.

Peak positions in the density of states can remain finite at $T_c$.

Different behaviors (gap filling vs. Fermi arc) depend on thermal contributions.)

Abstract

Near a quantum-critical point in a metal a strong fermion-fermion interaction, mediated by a soft boson, destroys fermionic coherence and also gives rise to superconductivity. In a class of large $N$ models, one would naively expect an incoherent (non-Fermi liquid) normal state behavior to persist down to $T=0$. However, this is not the case for quantum-critical systems described by Eliashberg theory, where non-Fermi liquid part of the self-energy $\Sigma (\omega_m)$ is large for a generic Matsubara frequency $\omega_m = \pi T(2m+1)$, but vanishes at $\omega_m = \pm \pi T$, while the pairing interaction between fermions with these two frequencies remains strong. This peculiarity gives rise to a non-zero $T_c$ even at large $N$[Y. Wang et al PRL 117, 157001 (2016)]. We consider the system behavior below $T_c$ and contrast the two cases when $\omega_m = \pm \pi T$ are either special or not. We obtain the solution of the non-linear gap equations in real frequency axis and obtain the spectral function $A(\omega)$ and the density of states $N(\omega)$. In a conventional BCS-type superconductor $A(\omega)$ and $N(\omega)$ are peaked at the gap value $\Delta (T)$, and the peak position shifts to a smaller $\omega$ as temperature increases, i.e. the gap "closes in". However when superconductivity is induced by fermions with $\omega_m = \pm \pi T$, the peak in $N(\omega)$ remains at a finite frequency even at $T =T_c-0$, the gap just "fills in". The spectral function $A(\omega)$ either shows almost the same "gap filling" behavior, or its peak position shifts to zero frequency already at a finite $\Delta$ ("emergent Fermi arc" behavior), depending on the strength of the thermal contribution. We compare our results with the data for the cuprates and argue that "gap filling" behavior holds in the antinodal region, while the "emergent Fermi arc" behavior holds in the nodal region.

Figures

Equations

$$N(\omega )=N_0\mathrm{Re}\left[\frac{1}{\sqrt{1-{\left(\Phi (\omega )/(\omega +\Sigma (\omega ))\right)}^2}}\right]$$

$$A(\omega )$$

$$L(\omega )$$

$$N(\omega )=N_0\mathrm{Re}\left[\frac{1}{\sqrt{1-{\left(\frac{\Delta (T)}{\omega +i\Gamma (T)}\right)}^2}}\right]$$

$$A(\omega )$$

$$\Phi ({\omega }_m)$$

$$\tilde{\Sigma }({\omega }_m)$$

$$\Delta ({\omega }_m)={\omega }_m\frac{\Phi ({\omega }_m)}{\tilde{\Sigma }({\omega }_m)}$$

$$\Delta ({\omega }_m)=\pi T{g}^{\gamma }\sum _{m^{\prime }}\frac{\Delta ({\omega }_{m^{\prime }})-\Delta ({\omega }_m)\frac{{\omega }_{m^{\prime }}}{{\omega }_m}}{\sqrt{{\omega }_{m^{\prime }}^2+{\Delta }^2({\omega }_{m^{\prime }})}}\frac{1}{|{\omega }_m-{\omega }_{m^{\prime }}{|}^{\gamma }}.$$

$${\Phi }^{\ast }({\omega }_m)$$

$${\tilde{\Sigma }}^{\ast }({\omega }_m)$$

$$Q({\omega }_m)$$

$$\Phi ({\omega }_m)$$

$$\tilde{\Sigma }({\omega }_m)$$

$$Q^{\ast }({\omega }_m)$$

$${\Phi }^{\ast }({\omega }_m)=\frac{\pi T}{N}{g}^{\gamma }\sum _{m^{\prime }\ne n}\frac{{\Phi }^{\ast }({\omega }_{m^{\prime }})}{\sqrt{({\tilde{\Sigma }}^{\ast }({\omega }_{m^{\prime }}){)}^2+({\Phi }^{\ast }({\omega }_{m^{\prime }}){)}^2}}\frac{1}{|{\omega }_m-{\omega }_{m^{\prime }}{|}^{\gamma }},$$

$${\tilde{\Sigma }}^{\ast }({\omega }_m)={\omega }_m+g^{\gamma }\pi T\sum _{m^{\prime }\ne m}\frac{{\tilde{\Sigma }}^{\ast }({\omega }_m)}{\sqrt{({\tilde{\Sigma }}^{\ast }({\omega }_{m^{\prime }}){)}^2+({\Phi }^{\ast }({\omega }_{m^{\prime }}){)}^2}}\frac{1}{|{\omega }_m-{\omega }_{m^{\prime }}{|}^{\gamma }},$$

$$\Delta ({\omega }_m)=\frac{\pi T}{N}{g}^{\gamma }\sum _{m^{\prime }\ne m}\frac{\Delta ({\omega }_{m^{\prime }})-N\Delta ({\omega }_m)\frac{{\omega }_{m^{\prime }}}{{\omega }_m}}{\sqrt{{\omega }_{m^{\prime }}^2)+{\Delta }^2({\omega }_{m^{\prime }})}}\frac{1}{|{\omega }_m-{\omega }_{m^{\prime }}{|}^{\gamma }}.$$

$$F_{sc}=-N_0\left(2\pi T\sum _m\frac{{\omega }_m^2}{\sqrt{{\omega }_m^2+{\Delta }_m^2}}+{\pi }^2{T}^2{g}^{\gamma }\sum _{m\ne m^{\prime }}\frac{{\omega }_m{\omega }_{m^{\prime }}+\frac{1}{N}{\Delta }_m{\Delta }_{m^{\prime }}}{\sqrt{{\omega }_m^2+{\Delta }_m^2}\sqrt{{\omega }_{m^{\prime }}^2+{\Delta }_{m^{\prime }}^2}}\frac{1}{|{\omega }_m-{\omega }_{m^{\prime }}{|}^{\gamma }}\right)$$

$$F_n=-N_0\left(2\pi T\sum _m|{\omega }_m|+{\pi }^2{T}^2{g}^{\gamma }\sum _{m\ne m^{\prime }}\frac{\mathrm{sgn}{\omega }_m\mathrm{sgn}{\omega }_{m^{\prime }}}{|{\omega }_m-{\omega }_{m^{\prime }}{|}^{\gamma }}\right)$$

$$\delta F=F_{sc}-F_n$$

$$\delta F=$$

$$+\frac{3}{4}\pi T{N}_0\sum _m|{\omega }_m|D_m^4-N_0{\pi }^2{T}^2{g}^{\gamma }\sum _{m\ne m^{\prime }}\frac{\mathrm{sgn}{\omega }_m\mathrm{sgn}{\omega }_{m^{\prime }}}{|{\omega }_m-{\omega }_{m^{\prime }}{|}^{\gamma }}$$

$$\times \left(\frac{1}{4}{D}_m^2{D}_{m^{\prime }}^2+\frac{3}{8}\left(D_m^4+D_{m^{\prime }}^4\right)-\frac{1}{2N}{D}_m{D}_{m^{\prime }}\left(D_m^2+D_{m^{\prime }}^2\right)\right)$$

$${\Phi }^{\ast }({\omega }_m)$$

$${\Sigma }^{\ast }({\omega }_m)$$

$$N_{cr}=(1-\gamma )\Gamma (\gamma /2)\left[\frac{\Gamma (\gamma /2)}{2\Gamma (\gamma )}+\frac{\Gamma (1-\gamma )}{\Gamma (1-\gamma /2)}\right],$$

$${\Phi }^{\ast }(\pi T)\approx \frac{1}{N}{\left(\frac{g}{2\pi T}\right)}^{\gamma }{\Phi }^{\ast }(-\pi T)$$

$${\Phi }^{\ast }({\omega }_m)=\frac{1}{N}{\left(\frac{g}{2\pi T}\right)}^{\gamma }\left(\frac{{\Phi }^{\ast }(\pi T)}{|\frac{1}{2}-\frac{{\omega }_{m^{\prime }}}{2\pi T}{|}^{\gamma }}+\frac{{\Phi }^{\ast }(-\pi T)}{|\frac{1}{2}+\frac{{\omega }_{m^{\prime }}}{2\pi T}{|}^{\gamma }}\right)$$

$${\Phi }_0^{\ast }$$

$${\tilde{\Sigma }}_0^{\ast }$$

$${\Phi }_0^{\ast }$$

$${\tilde{\Sigma }}_0^{\ast }$$

$${\Delta }_0=\pi T{\left(\frac{2}{N}\right)}^{1/2}{\left(1-{\left(\frac{T}{T_c}\right)}^{\gamma }\right)}^{1/2}$$

$${\Phi }^{\ast }({\omega }_m)\approx {\Phi }_0^{\ast }\left[{\left(\frac{2\pi T}{|{\omega }_m-\pi T|}\right)}^{\gamma }+{\left(\frac{2\pi T}{|{\omega }_m+\pi T|}\right)}^{\gamma }\right]$$

$${\Sigma }^{\ast }({\omega }_m)\approx 2N{\tilde{\Sigma }}_0^{\ast }H\left(\frac{|{\omega }_m|-\pi T|}{2\pi T},\gamma \right)\mathrm{sgn}(m+\frac{1}{2})$$

$${\Phi }^{\ast }({\omega }_m)\approx \frac{2{\Phi }_0}{|m{|}^{\gamma }},{\tilde{\Sigma }}^{\ast }({\omega }_m)\approx 2N\frac{|{\tilde{\Sigma }}_0|}{1-\gamma }|m{|}^{1-\gamma }\mathrm{sgn}(m)$$

Full text

I Introduction.

The pairing near a quantum-critical point (QCP) in a metal is a fascinating subject due to highly non-trivial interplay between superconductivity and non-Fermi liquid (NFL) behavior  Combescot (1995); Ber ; *Bergmann2; *ad; Marsiglio et al. (1988); *Marsiglio_91; Karakozov et al. (1991); Bonesteel et al. (1996); Abanov et al. (2001a, 2003); *acs2; Abanov et al. (2001b); Son (1999); *son2; Lee (2009); *sslee2; Sachdev et al. (2009); *subir2; Moon and Chubukov (2010); Metlitski and Sachdev (2010a); *max2; Mross et al. (2010); Mahajan et al. (2013); *raghu2; *raghu3; *raghu4; *raghu5; Monthoux et al. (2007); *scal2; boo ; rev ; Fratino et al. (2016); Metlitski et al. (2015); Efetov et al. (2013); *efetov; Wang and Chubukov (2013a); *wang2; Raghu et al. (2015); Wang et al. (2016); Lederer et al. (2015); Tsvelik (2017); *Rice; Vojta and Sachdev (1999); Fradkin et al. (2010); Bok et al. (2016); Shibauchi et al. (2014); Vilardi et al. (2018); *Metzner18; Gerlach et al. (2017); *berg_2; *berg_3; Haule and Kotliar (2007); *kotliar2; Fratino et al. (2016); *tremblay_2; Georges et al. (2013); *georges2; Khveshchenko and Shively (2006); Lee et al. (2018); *we_last_D. In most cases, the dominant interaction between low-energy fermions near a QCP is mediated by critical fluctuations of the order parameter. In dimensions $D\leq 3$, this interaction gives rise to a singular fermionic self-energy, and a coherent Fermi-liquid behavior get destroyed below a certain temperature $T_{coh}$, either on the full Fermi surface Rech et al. (2006); Metlitski and Sachdev (2010a); Mahajan et al. (2013); q= ; *q=01; *q=02 or in the hot regions Abanov et al. (2001a, 2003); *acs2; Abanov et al. (2001b, 2008); Metlitski and Sachdev (2010a); Wang and Chubukov (2013a); *wang2; Altshuler et al. (1995); *2kf2; *2kf3. The same interaction, however, also mediates fermion-fermion interaction in the particle-particle channel. The electron-mediated interaction is positive (repulsive), but it depends on both momentum and frequency and generally has at least one attractive component ($d-$wave for antiferromagnetic QCP, $p-$wave for a ferromagnetic QCP, $s,p,d$-wave for a nematic QCP, Ref. Scalapino (2012b); *review2; *review3; *review4) If this system becomes superconducting below some finite $T_{c}$, the range of NFL behavior shrinks to $T_{coh}>T>T_{c}$, and even vanishes when $T_{c}>T_{coh}$ Metlitski et al. (2015). A naked quantum-critical $T=0$ behavior can only be observed either if the pairing interaction is repulsive, or if fermionic incoherence prevents superconductivity to develop down to $T=0$.

In all known physical quantum-critical (QC) models of fermions, superconducting $T_{c}$ is finite Bonesteel et al. (1996); Abanov et al. (2001a); Metlitski et al. (2015); Raghu et al. (2015); Wang et al. (2016); Khveshchenko and Shively (2006); Lee et al. (2018); *we_last_D. This can be interpreted as an evidence that the tendency to pairing is stronger than towards incoherent, NFL behavior. The situation can potentially be reversed if the interaction in the pairing channel is somehow reduced compared to that in the particle-hole channel. This can be achieved by either modifying the momentum dependence of the interaction mediated by critical fluctuations to reduce the partial pairing component in the cannel, where it is attractive, or by keeping the interaction intact but extending the model to an $SU(N)$ global symmetry Raghu et al. (2015) (the original model corresponds to $N=1$). Under this extension, the pairing interaction get reduced by $1/N$, but the self-energy stays intact Raghu et al. (2015). In both cases, the functional form of equation for the (frequency dependent) pairing vertex in the attractive channel does not change, but the magnitude of the eigenvalue needed for superconductivity gets larger. The analysis of a large-$N$ QC model at $T=0$ shows Raghu et al. (2015); Wang et al. (2016) that there exists a critical $N_{cr}$, separating a superconducting region at $N<N_{cr}$ and a region of a $T=0$ NFL normal state behavior at $N>N_{cr}$ (see Fig. 1). A conventional reasoning in this situation would be that the superconducting $T_{c}(N)$ terminates at $T=0$, $N=N_{cr}$, and vanishes for $N>N_{cr}$. However, numerical studies of large-$N$ QC models yield a different result Wang et al. (2016) – $T_{c}$ remains finite at any $N$, and the critical line $T_{c}(N)$ by-passes $N=N_{cr}$, and $T_{c}(N)$ remains finite at all $N$ (see Fig. 2).

This unusual behavior was argued in Ref. Wang et al., 2016 to be the consequence of the special form of Matsubara fermionic self-energy $\Sigma(\omega_{m})$ at the two lowest Matsubara frequencies: $\omega_{m}=\pi T$ and $\omega_{m}=-\pi T$. Namely, in Eliashberg theory $\Sigma(\pm\pi T)$ only contains the self-action term (thermal contribution to $\Sigma(\omega_{m})$ from zero bosonic Matsubara frequency), all other contributions cancel out. The thermal piece in $\Sigma(\omega_{m})$ comes from scattering with zero frequency and finite momentum transfer and mimics the scattering by impurities. The same thermal scattering also contributes to the pairing vertex $\Phi(\omega_{m})$. For spin-singlet pairing, the two contributions cancel out in equation for the gap function $\Delta(\omega_{m})=\Phi(\omega_{m})/(1+\Sigma(\omega_{m})/\omega_{m})$ by Anderson’s theorem Anderson (1959); Abrikosov and Gor’kov (1959). As a consequence, fermions with $\omega_{m}=\pm\pi T$ can be treated for the pairing as free quasiparticles. Meanwhile the pairing interaction between fermions with $\omega_{m}=\pi T$ and $\omega_{m}=-\pi T$ remains strong. This strong interaction, not countered by the self-energy, gives rise to the emergence of $\Delta(\pm\pi T)$ below a certain $T_{c}(N)$, which remains finite for all values of $N$. A finite $\Delta(\pm\pi T)$ then induces non-zero $\Delta(\omega_{m})$ at other Matsubara frequencies, for which the self-energy without self-action is strong.

In this communication we extend the analysis of superconductivity induced by first fermionic Matsubara frequencies to $T<T_{c}(N)$. We argue that, although $T_{c}(N)$ by-passes $N=N_{cr}$, there is a crossover in the system behavior at $T_{cross}(N)<T_{c}(N)$. The crossover line $T_{cross}(N)$ originates at $T=0$ for $N=N_{cr}$ and ends at $T_{cross}\leq T_{c}$ for the physical case $N=1$. At $T_{cross}(N)<T<T_{c}(N)$, superconductivity can be viewed as induced by fermions with $\omega_{m}=\pm\pi T$, at smaller $T<T_{cross}(N)$ fermions with all $\omega_{m}$ contribute to superconductivity, and the ones with $\omega_{m}=\pm\pi T$ are no longer special. We show the schematic phase diagram in Fig. 3.

We analyze the evolution of the gap $\Delta(\omega_{m})$ below $T_{c}(N)$ at $N>N_{cr}$ and $N<N_{cr}$ and then convert from Matsubara to real frequencies and analyze the behavior of $\Delta(\omega)$, the spectral function at the Fermi surface $A(\omega)$, and the density of states (DOS) $N(\omega)$. We argue that the system behavior below $T_{c}(N)$ is different for $N>N_{cr}$ and $N<N_{cr}$ (see the paths (a) and (b) in Fig.3). At $N>N_{cr}$ it is qualitatively different from BCS. At $N<N_{cr}$, the system behavior is similar to a BCS superconductor for $T<T_{cross}(N)$ and to that at $N>N_{cr}$ for $T>T_{cross}(N)$.

Along the Matsubara axis, we find that at large $N>N_{cr}$, the pairing vertex $\Phi(\omega_{m})$ is smaller than $\Sigma(\omega_{m})$ for all temperatures and all Matsubara frequencies, including $\pm\pi T$. In fact, $\Sigma(\omega_{m})$ with $m\neq 0,-1$ remains essentially the same as in the normal state, i.e., the feedback effect from superconductivity on this self-energy is weak. The self-energy $\Sigma(\pm\pi T)$ becomes finite below $T_{c}(N)$, but remains smaller by $1/N$ than $\Sigma(\omega_{m})$ at other Matsubara frequencies. Still, it is larger by $\sqrt{N}$ than $\Phi(\pm\pi T)$. We show that in this situation, $\Delta(\omega_{m},T)$ is monotonic as a function of $\omega_{m}$, with the largest value at $\pm\pi T$, but non-monotonic as a function of temperature, i.e., $\Delta(\pi T)$ first increases when $T$ decreases below $T_{c}(N)$, then passes through a maximum and eventually vanishes at $T=0$. At $N<N_{cr}$, $\Delta(\pi T)$ becomes non-zero at $T=0$, and its magnitude increases as $N$ gets progressively smaller than $N_{cr}$. At $N\leq N_{cr}$, the temperature dependence of $\Delta(\pi T)$ is still non-monotonic, with the maximum at a finite $T$. At smaller $N$, the maximum becomes more shallow, and at $N\gtrsim 1$, $\Delta(\pi T)$ monotonically increases with decreasing $T$.

We use the results along the Matsubara axis as an input and obtain the behavior of $\Phi(\omega)$ and $\Sigma(\omega)$ along real frequency axis. Using these $\Phi(\omega)$ and $\Sigma(\omega)$, we obtain the DOS

[TABLE]

The thermal contributions to $\Phi(\omega)$ and to $\omega+\Sigma(\omega)$ are the same and they cancel out in the DOS, i.e., in the calculations one can replace $\Phi(\omega)$ and $\Sigma(\omega)$ by $\Phi^{*}(\omega)$ and $\Sigma^{*}(\omega)$, which are the solutions of the Eliashberg equations with thermal contributions explicitly taken out. We show that, for $N>N_{cr}$, $N(\omega)$ is finite for all frequencies, including $\omega=0$, and its dependence on $\omega$ is determined by a universal scaling function of $\omega/T$. As the consequence, the frequency at which $N(\omega)$ has a maximum, linearly increases with increasing $T$. As $T$ approaches $T_{c}$ from below, DOS “fills in”, i.e., the $N(\omega)$ approaches $N_{0}$, but the position of the maximum in $N(\omega)$ remains at a finite frequency.

At $N<N_{cr}$, the DOS $N(\omega)$ at the lowest $T<T_{cross}(N)$ displays a sharp gap, , i.e., it nearly vanishes at $\omega<\omega_{0}$, where $\omega_{0}$ is roughly equal to $\Delta(0)$ (more exactly, $\omega_{0}$ is the solution of $\Delta(\omega_{0})=\omega_{0}$). As $T$ increases, the position of the maximum in the DOS initially shifts to a lower frequency, as in a BCS superconductor, because $\Delta(0)$ gets smaller with increasing $T$, i.e., the gap in the DOS “closes in” with increasing temperature. However, once temperature exceeds $T_{cross}(N)$, this behavior changes and becomes the same as for larger $N$, i.e., at these $T$ the the position of the maximum in DOS shifts to a higher frequency with increasing $T$ and remains finite at $T=T_{c}(N)-0$, i.e., the DOS “fills in” with increasing $T$. We emphasize that these two distinct regimes of system behavior are present also in the original physical model with $N=1$. In this respect, the extension to $N>1$ is just a convenient way to understand the origin of such behavior by extending the regime in which superconductivity is generated by fermions with $\omega=\pm\pi T$. A representative of our results for the DOS is shown in Fig.4

The behavior of the spectral function is more involved because in $A(\omega)$ the thermal contribution does not cancel out. The expression for $A(\omega)=-(1/\pi)\operatorname{Im}[G(k_{F},\omega)]$ at $\omega>0$ is (see Eq.(68) below)

[TABLE]

where frequency-independent $P=P(T)$ describes the thermal contribution to self-energy, and $\Sigma^{*}(\omega)$, $\Phi^{*}(\omega)$ are obtained from $\Sigma(\omega)$, $\Phi(\omega)$ by excluding thermal contributions (see (44) and the Appendix for more details). For large $P$, $A(\omega)\propto N(\omega)$, i.e., the spectral function displays the same crossover from “gap closing” to “gap filling” as the DOS. For smaller $P$, when the term next to $P$ in (2) is larger than $P$, $A(\omega)$ at $T<T_{cross}$ shows two sharp peaks at $\omega=\pm\omega_{0}$. At temperatures above $T_{cross}$, the two peaks merge, and $A(\omega)$ develops a maximum at $\omega=0$, like in the normal state. A representative of our results for $A(\omega)$ is shown in Fig.5.

The transformation from “gap closing” to “gap filling” behavior in the DOS has been observed in several superconducting materials, most notably the cuprates Fischer et al. (2007); Reber et al. (2012); Kondo et al. (2013); *Kaminski2; Kanigel et al. (2006); *kanigel2; *kanigel3; nor ; Damascelli et al. (2003); *shen2; Johnson et al. (2001); Kordyuk and Borisenko (2006); *kordyuk2; He et al. (2014); Peng et al. (2013) The spectral function in the cuprates shows the same behavior as the DOS in the antinodal regions, where the fermionic incoherence is the strongest, and the $d-$wave gap is the largest. In the regions near the Brillouin zone diagonals, the symmetrized spectral function has peaks at a finite frequency $\pm\omega_{0}$ at low temperatures, and a single maximum at $\omega=0$ at higher temperatures. The angular range in which the system displays a single peak above a certain $T$ is termed as a Fermi arc nor .

The crossover from “gap closing” to “gap filling” in the DOS and in $A(\omega)$ in the antinodal regions, and the crossover from two peaks to a single peak in $A(\omega)$ in the nodal regions, have been phenomenologically described by assuming that the pairing vertex $\Phi^{*}(\omega)=\Delta(T)$, as in a BCS superconductor, and $\Sigma^{*}(\omega)=i\Gamma(T)$ (Refs. Balatsky et al. (2006); *imp_general2; *imp_general3; Norman et al. (1998a); *imp_cuprates2; *imp_cuprates3; *imp_cuprates4; *imp_cuprates5; Chubukov et al. (2007)) Under these approximation, the DOS becomes

[TABLE]

Without $\Gamma(T)$, the DOS vanishes at $\omega<\Delta$ and is singular at $\omega=\Delta+0$. A non-zero $\Gamma(T)$ makes $N(\omega)$ continuous and non-zero down to $\omega=0$. Furthermore, the position of the peak in $N(\omega)$ shifts to a higher frequency from $\omega=\Delta(T)$ (see Fig.6) At vanishing $\Delta(T)$ the peak in $N(\omega)\approx N_{0}\left(1+\frac{1}{2}\Delta^{2}\operatorname{Re}\left[\frac{1}{(\omega+i\Gamma)^{2}}\right]\right)$ remains at a finite $\omega=\sqrt{3}\Gamma$. In other words, the magnitude of the deviation of $N(\omega)$ from $N_{0}$ is set by $\Delta^{2}$, while its frequency dependence is set by $\operatorname{Re}\frac{1}{(\omega+i\Gamma)^{2}}$ and does not depend on $\Delta$. If one additionally sets $\Gamma=O(T)$, as in marginal FL theory, one obtains that the position of the maximum in the DOS increases linearly with $T$ near $T=T_{c}$, when $\Gamma>\Delta(T)$. The same phenomenolgical model with $\Gamma(T)\propto T$ was used Chubukov et al. (2007) to explain Fermi arcs (assuming that the thermal contribution can be neglected). Indeed, at $P=0$, we have from (2)

[TABLE]

This spectral function has two separate peaks at positive and negative $\omega$ at $\Gamma(T)<\sqrt{3}\Delta(T)$, and a single maximum at $\omega=0$ at $\Gamma(T)>\sqrt{3}\Delta(T)$ (Fig.6)

This phenomenon when $N(0)$ becomes finite is known as “gapless superconductivity”. It was originally found by Abrikosov and Gorkov in their analysis of an s-wave BCS superconductor with magnetic impurities Abrikosov and Gor’kov (1961). At $T=0$, gapless superconductivity exists in a finite parameter range before magnetic impurities destroy superconductivity. Several researchers later argued Maki (1964); *maki2; *maki3; *maki4 that any phonon-mediated $s$-wave superconductor at a finite $T$ is a gapless superconductor due to scattering on thermally excited phonons, although in practice $\Gamma$ due to such scattering is extremely small at small coupling. For electronically-mediated superconductivity in a clean metal, self-energy $\Sigma(\omega)$ in the normal state contains the imaginary part. In a superconducting state, the imaginary part of $\Sigma(\omega)$ is reduced at $\omega<\Delta$ due to the reduction of the phase space for low-energy scattering. This holds for any symmetry of the gap function and gives rise to peak-dip-hump feature of the spectral functions, studied extensively in the cuprates nor ; Fink et al. (2006); Abanov and Chubukov (1999b); Eschrig (2006) As long as $T$ is finite, $\Sigma^{\prime\prime}(\omega=0)$ remains non-zero, but at low $T$ it gets substantially reduced compared to its value at $T=T_{c}$. Numerical analysis of Eliashberg equations for several models of magnetically-induced $d$-wave superconductivity Abanov et al. (2003, 2008); Monthoux et al. (2007) and for strong coupling (small Debye frequency) limit of electron-phonon superconductivity Combescot (1995); Ber ; *Bergmann2; *ad; Marsiglio et al. (1988); *Marsiglio_91; Karakozov et al. (1991) did find that $\Sigma^{\prime\prime}(0)$ rapidly increases at $T$ near $T_{c}$, and the maximum in the theoretical DOS shifts up from $\Delta(T)$ and remains at a finite frequency at $\sim T_{c}$, where $\Delta(T)$ vanishes in the Eliashberg theory. This is roughly consistent with the phenomenology of Eq. (3), although temperature variation of the peak position in the DOS has not been explicitly verified.

We view our results as the microscopic explanation of the rapid increase of $\Sigma^{\prime\prime}(0)$ above a certain $T$ within the superconducting state and the related transformation from “gap closing” to “gap filling” behavior of the DOS and the spectral function at large $P$ (and the transformation from gap to Fermi arc behavior at smaller $P$). To reiterate – we argue that the conventional “gap closing” behavior occurs at $T<T_{c,1}$, while the “gap filling” behavior occurs at $T_{c,1}<T<T_{c}$, where in Matsubata formalism the pairing is induced by two lowest Matsubata frequencies at which self-energy vanishes, and would not happen if fermions with these two frequencies were eliminated from the gap equation.

The issue which we do not address here is the role of pairing fluctuations. We remind the reader that Eliashberg theory neglects phase and amplitude fluctuations of the pairing vertex and in this respect should be treated as effectively a “mean-field” theory. It is very likely that in some range below Eliashberg $T_{c}$ fluctuations destroy long-range superconducting order, and the actual $T_{c,act}<T_{c}$. Our results, that the DOS $N(\omega)$ is non-zero at all $\omega$ and the position of its maximum increases with $T$, should survive at $T_{c,act}<T<T_{c}$ as our reasoning only explores the fact that in this $T$ range the feedback from the pairing on the fermionic self-energy is weak. Gap fluctuations reduce this feedback even further. The same holds for the spectral function both at large $P$ and at smaller $P$. In other words, our theory describes gap filling and Fermi arcs in the pseudogap region. Still, to fully address the issue of gap fluctuations one needs to go beyond Eliashberg theory and analyze the full Luttinger-Ward functional Luttinger and Ward (1960).

The paper is organized as follows. In Sec. II we present the microscopic model of pairing mediated by a gapless boson with $\chi(\Omega_{m})=(g/|\Omega_{m}|)^{\gamma}$ (the $\gamma$-model) and its extension to $N>1$. We present the set of coupled Eliashberg equations along Matsubara axis for the pairing vertex $\Phi(\omega_{m})$ and the fermionic self-energy $\Sigma(\omega_{m})$ and summarize, in Sec. II.1 earlier results of the analysis of the linearized equation for $\Phi(\omega_{m})$. At $T=0$, these results show that there exists the critical $N_{cr}$, separating the superconducting state at $N<N_{cr}$ and the normal state at $N>N_{cr}$. At $T>0$, these calculations show that superconductivity emerges for all $N$, below a certain $T_{s}(N)$ which only vanishes at $N=\infty$. In Sec. III we discuss the system behavior at $N>N_{cr}$, first in Matsubara frequencies, in Sec. III.1, and then in real frequencies, in Sec. III.3. We present the analytical solution of the Eliashberg equations at large $N$ and discuss the behavior of the gap, the Free energy and the specific heat, the DOS, and the spectral function. In Sec. IV we discuss system behavior at $N<N_{cr}$, again first in Matsubara frequencies, in Sec. IV.1, and then in real frequencies, in Sec.IV.2. In Sec. V we summarize our results and compare them with the experimental data.

II The model.

We consider a model of itinerant fermions at the onset of a long-range order in either spin or charge channel. At the critical point the propagator of a soft boson becomes massless and mediates singular interaction between fermions. We follow earlier works Abanov et al. (2001a, 2003); Moon and Chubukov (2010); Metlitski and Sachdev (2010a); Mross et al. (2010); Monthoux et al. (2007); Efetov et al. (2013); Metlitski et al. (2015); Raghu et al. (2015); Haslinger and Chubukov (2003); Wang et al. (2016); Lee et al. (2018) and assume that this interaction as attractive in at least one pairing channel and that bosons can be treated as slow modes compared to fermions, i.e., the Eliashberg approximation is valid. Within this approximation one can explicitly integrate over the momentum component perpendicular to the Fermi surface (for a given pairing symmetry) and reduce the pairing problem to a set of coupled integral equations for frequency dependent self-energy $\Sigma(\omega_{m})$ and the pairing vertex $\Phi(\omega_{m})$ for fermions on the Fermi surface, with effective frequency-dependent dimensionless interaction $\chi(\Omega)=(g/|\Omega|)^{\gamma}$ (the $\gamma$-model, Refs. Abanov et al. (2003, 2001a); Moon and Chubukov (2010); Wang et al. (2016); Lee et al. (2018)). This interaction simultaneously gives rise to NFL form of the self-energy in the normal state and to pairing. The equations we analyze are

[TABLE]

where ${\tilde{\Sigma}}(\omega_{m})=\omega_{m}+\Sigma(\omega_{m})$. Note that we define $\Sigma(\omega_{m})$ as a real function of frequency, i.e., without the overall factor of $i$. The self-energy along Matsubara axis, related by Kramers-Krong (KK) formula to $\Sigma^{{}^{\prime\prime}}(\omega)$ along the real frequency axis, does contain the factor $i$. The superconducting gap $\Delta(\omega_{m})$ is defined as a real variable

[TABLE]

The equation for $\Delta(\omega)$ is readily obtained from (II):

[TABLE]

This equation contains a single function $\Delta(\omega)$, but for the prize that $\Delta(\omega_{m})$ appears on both sides of the equation, which makes (7) less convenient for the analysis than Eqs. (II).

The r.h.s. of the equations for $\Phi(\omega_{m})$ and $\Sigma(\omega_{m})$ contain divergent pieces from the terms with $m^{\prime}=m$, i.e., from $\chi(0)$. The divergence can be regularized by moving slightly away from a QCP, in which case $\chi(0)$ is large but finite. This term mimics the effect of non-magnetic impurities. To get rid of the thermal piece in the equations for $\Phi(\omega)$ and $\Sigma(\omega)$, we follow Millis et al. (1988); Abanov et al. (2008) and use the same trick as for the derivation of the Anderson theorem for impurity scattering agd Namely, we pull out the term with $m^{\prime}=m$ from the sum, move it to the l.h.s., and introduce

[TABLE]

The ratio $\Phi(\omega_{m})/{\tilde{\Sigma}}(\omega_{m})=\Phi^{*}(\omega_{m})/{\tilde{\Sigma}}^{*}(\omega_{m})$, hence $\Delta(\omega_{m})$, defined in (6), is invariant under $\Phi(\omega_{m})\to\Phi^{*}(\omega_{m})$ and ${\tilde{\Sigma}}(\omega_{m})\to{\tilde{\Sigma}}^{*}(\omega_{m})$. Using (8), one can easily verify that the equations on $\Phi^{*}(\omega_{m})$ and ${\tilde{\Sigma}}^{*}(\omega_{m})$ are the same as in (II), but without the thermal piece, i.e., the summation over $m^{\prime}$ now excludes the divergent term with $m^{\prime}=m$. The gap function $\Delta(\omega_{m})$, defined in (6) is equally expressed in terms of $\Phi^{*}(\omega_{m})$ and ${\tilde{\Sigma}}^{*}(\omega_{m})$, and the gap equation (7) preserves its form: the sum over $m^{\prime}$ now excludes the term with $m^{\prime}=m$, but this term vanishes anyway because the numerator in the r.h.s. of (7) vanishes at $m^{\prime}=m$. One can also solve (8) backwards and express $\Phi(\omega_{m})$ and ${\tilde{\Sigma}}(\omega_{m})$ via $\Phi^{*}(\omega_{m})$ and ${\tilde{\Sigma}}^{*}(\omega_{m})$ as

[TABLE]

Eq. (II) describes color superconductivity Son (1999) ($\gamma=0_{+}$, $\chi(\Omega_{m})\propto\log{|\omega_{m}|}$), spin- and charge-mediated pairing in $D=3-\epsilon$ dimension Mross et al. (2010); Metlitski et al. (2015); Raghu et al. (2015) ($\gamma=O(\epsilon)\ll 1$), a 2D pairing  Altshuler et al. (1995) with interaction peaked at $2k_{F}$ ($\gamma=1/4$), pairing at a 2D nematic/Ising-ferromagnetic QCP Bonesteel et al. (1996); Lederer et al. (2015); Wang et al. (2001); *triplet2; *triplet3 ($\gamma=1/3$), pairing at a 2D $(\pi,\pi)$ SDW QCP Abanov et al. (2001a, 2003); Millis (1992); Wang and Chubukov (2013a) and an incommensurate CDW QCP Castellani et al. (1995); *ital2; *ital3; Chowdhury and Sachdev (2014a); *wang_22; *wang23 ($\gamma=1/2$), a 2D pairing mediated by an undamped propagating boson ($\gamma=1$), and the strong coupling limit of phonon-mediated superconductivity Combescot (1995); Ber ; *Bergmann2; *ad; Marsiglio et al. (1988); *Marsiglio_91; Karakozov et al. (1991) ($\gamma=2$). The pairing models with parameter-dependent $\gamma$ have also been considered (Refs. Sachdev et al., 2009; Moon and Chubukov, 2010). In this communication we consider the set of $\gamma$-models with $\gamma<1$. The analysis for $\gamma>1$ requires a separate consideration because of the divergence of the normal state self-energy at $T=0$.

The full set of Eliashberg equations for electron-mediated pairing contains also the equation describing the feedback from the pairing on $\chi(\Omega)$, e.g., the emergence of a propagating mode (often called a resonance mode) in the dynamical spin susceptibility for $d-$wave pairing mediated by antiferromagnetic spin fluctuations. To avoid additional complications, we do not include this feedback into our consideration. In general terms, the feedback from the pairing makes bosons less incoherent and can be modeled by assuming that $\gamma$ moves towards $\gamma=1$ as $T$ moves down from $T_{c}$.

The two equations in (II) describe the interplay between two competing tendencies – the tendency towards superconductivity, specified by $\Phi$, and the tendency towards incoherent non-Fermi liquid behavior, specified by $\Sigma$. The competition between the two tendencies is encoded in the fact that $\Sigma$ appears in the denominator of the equation for $\Phi$ and $\Phi$ appears in the denominator of the equation for $\Sigma$. Accordingly, a large, non-FL self-energy is an obstacle to Cooper pairing, while once $\Phi$ develops, it reduces the strength of the self-energy, i.e., moves a system back into a FL regime. Like we said in the Introduction, our goal is to analyze the special role of fermions with Matsubara frequencies $\omega_{m}=\pm\pi T$ in the situation when the tendency towards pairing is reduced compared to that for NFL normal state. For this, we extend the model to matrix $SU(N)$. Under this extension, the interaction in the particle-hole channel, which gives rise to fermionic self-energy, remains intact, while the interaction in the particle-particle channel acquires an additional factor $1/N$. We emphasize that we extend to $N\neq 1$ after we invoke the analog of the Anderson theorem and eliminate the thermal contributions to $\Phi(\omega_{m})$ and $\Sigma(\omega_{m})$. In this respect our approach differs from the one in Ref. Raghu et al. (2015) There, the extension to large $N$ was done without first subtracting the thermal contributions. As a result, at a finite $N$ there appeared additional terms, singular at a QCP, which gave rise to qualitative changes in the system behavior. In our extension to $N>1$ these additional terms do not appear. Put it more simply, in our case after the extension the eigenvalues in the pairing channel get multiplied by $1/N$, i.e., a larger magnitude of the original eigenvalue is needed for superconductivity.

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

[TABLE]

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

[TABLE]

Below we will also need the expression for the Free energy $F_{sc}$ of a superconductor, described by the Eliashberg theory. The formula for $F_{sc}$ has been obtained in Refs. Luttinger and Ward (1960); Bardeen and Stephen (1964); Eliashberg (1960) in the studies of phonon-mediated superconductivity ($\gamma=2$ case at finite $\omega_{D}$ and $N=1$). Extending the results to $\gamma<1$, QC regime, and $N\neq 1$, we obtain

[TABLE]

where $\Delta_{m}=\Delta(\omega_{m})$. The gap equation (11) is obtained from the condition $\delta F_{sc}/\delta\Delta_{n}=0$ In the normal state the expression for the Free energy reduces to

[TABLE]

The difference between $F_{sc}$ and $F_{n}$ at $T=0$ is known as the condensation energy of a superconductor. At a finite $T$,

[TABLE]

where $D_{n}=D(\omega_{n})=\Delta(\omega_{n})/\omega_{n}$. Near $T=T_{c}$, one can expand $\delta F$ in powers of $\Delta_{m}$:

[TABLE]

II.1 Linearized gap equation

To obtain $T_{c}$ it is sufficient to consider the linearized gap equation. It is obtained from (II) by setting $\Phi^{*}$ to be infinitesimally small. Then $\Phi^{*}(\omega_{m^{\prime}})$ in the denominators of (II) can be ignored and the self energy $\Sigma^{*}(\omega_{m})$ is approximated by its normal state value. The resulting equations are:

[TABLE]

By power-lay counting we expect $\Sigma^{*}(\omega_{m})\propto g^{\gamma}\omega^{1-\gamma}$. Substituting this into the equation for $\Phi$ in (16) we obtain that at $|\omega_{m^{\prime}}|>|\omega_{m}|$, the pairing kernel $K=(g^{\gamma}/N)/(|\omega_{m^{\prime}}+\Sigma^{*}(\omega_{m^{\prime}})|)/|\omega_{m^{\prime}}|^{\gamma}$ is marginal at $|\omega_{m^{\prime}}|<g$: $K\propto 1/|\omega_{m^{\prime}}|$ (with prefactor independent on $g$), and decays as $K\propto g^{\gamma}/|\omega_{m^{\prime}}|^{1+\gamma}$ at $|\omega_{m^{\prime}}|>g$ This implies that $T_{c}$, if it exists, should be generally of order $g$. The marginal form of the kernel is similar to the BCS case and it gives rise to logarithmical growth of the pairing susceptibility within the perturbation theory. However, in distinction to BCS, the marginal form of $K$ holds only if $|\omega_{m^{\prime}}|>|\omega_{m}|$, i.e., at each order of perturbation the logarithm is cut by the running frequency in the next cross-section in the Cooper ladder. As the consequence, the summation of the logarithms alone does not lead to the divergence of the pairing susceptibility. In this situation, the conventional wisdom is that the pairing is the threshold phenomenon, i.e., it occurs if the pairing vertex exceeds some finite value. The pairing strength in Eq. (16) is controlled by $1/N$, hence by this logics there should be a critical $N_{cr}$ separating superconducting state at $N<N_{cr}$ and non-superconducting naked critical non-FL state at $N>N_{cr}$. At larger $N$ the tendency towards pairing is stronger than the tendency towards a non-FL behavior; at smaller $N$ the situation is the opposite. The analysis of the pairing problem at $T=0$ does yield exactly this king of behavior Raghu et al. (2015); Wang et al. (2016). Namely, there exists

[TABLE]

separating superconducting and non-superconducting states ($\Gamma(...)$ is a Gamma function). We plot $N_{cr}(\gamma)$ in Fig.7

The existence of $N_{cr}$ at $T=0$ would normally imply that this is the termination point of the line $T_{c}(N)$. However, the numerical solution of (16) yields qualitatively different result: $T_{c}$ is non-zero at any $N$, and the line $T_{c}(N)$ by-passes $N_{cr}$ and approaches zero only at $N\to\infty$ (see Fig.8). The reason for this behavior has been clarified in Ref. Wang et al. (2016). It turns out that power counting argument that $\Sigma^{*}(\omega_{m})\propto\omega^{1-\gamma}_{m}$, does not work for the first two Matsubara frequencies $\omega_{m}=\pm\pi T$, for which Eq. (16) yields $\Sigma^{*}(\pm\pi T)=0$. The reason is the presence of the sign-changing factor $\operatorname{sgn}(\omega_{m^{\prime}})$ in the r.h.s. of the formula for $\Sigma^{*}(\omega_{m})$. For $\omega_{m}=\pm\pi T$, contributions from positive and negative $\omega_{m^{\prime}}$ cancel out. To see the consequence of $\Sigma^{*}(\pm\pi T)=0$, consider the limit $N\gg 1$ and set external $\omega_{m}=\pi T$. For $\omega_{m^{\prime}}=O(T)$, but $\omega_{m^{\prime}}\neq-\pi T$, the product $\pi TK(\omega_{m^{\prime}})$ is independent of $T$ and is small in $1/N$. However, for $\omega_{m^{\prime}}=-\pi T$, this product is $\pi TK=(1/N)(g/(2\pi T))^{\gamma}$, and it becomes large at small enough $T$. A simple experimentation shows that in this situation the gap equation reduces to

[TABLE]

The last equation is for $\omega_{m}\neq\pm\pi T$. We will be searching for even-frequency solutions of the gap equation: $\Phi^{*}(\omega_{m})=\Phi^{*}(-\omega_{m})$. Then the first equation in (18) sets $T_{c}=(g/2\pi)1/N^{1/\gamma}$, and the second shows that a non-zero $\Phi^{*}(\omega_{m})$ is induced by $\Phi^{*}(\pm\pi T)$.

The functional form $T_{c}\propto 1/N^{1/\gamma}$ at large $N$ has been verified numerically in Ref. Wang et al. (2016) for a particular choice of $\gamma=0.1$. In Fig.8 we show that the same behavior holds for $\gamma=0.3$ and $0.9$. We now go beyond Ref. Wang et al. (2016) and verify that this behavior of $T_{c}$ (i.e., that $T_{c}(N)$ line by-passes $N_{cr}$) is indeed due to vanishing of the self-energy at the first two Matsubara frequencies. For this, we exclude $\omega_{m}=\pi T$ from the set of Matsubara frequencies and then solve again the linearized gap equation. The result is shown in Fig.8 We clearly see that ${\tilde{T}}_{c}$, obtained this way, tends to zero above some critical value of $N$, which numerically is close to $N_{cr}(\gamma)$ in Eq. (17). The outcome is that, without the first two Matsubara frequencies, the system would display a conventional behavior with ${\tilde{T}}_{c}(N)$ line terminating at a QCP at $N=N_{cr}$. At larger $N$, superconductivity would be absent because of stronger tendency towards a (competing) non-FL ground state. That the actual $T_{c}(N)$ by-passes $N_{cr}$ and vanishes only at $N=\infty$ is then entirely due to the vanishing of the self-energy for fermions with $\omega_{m}=\pm\pi T$.

The discrepancy between $T_{c}(N)$ and ${\tilde{T}}_{c}(N)$ suggests that physical properties below the actual onset temperature for the pairing $T_{c}(N)$ depend on whether $N$ is smaller or larger than $N_{cr}$. When $N>N_{cr}$, the pairing is induced by fermions with $\omega_{m}=\pm\pi T$ and, the order parameter $\Phi(\omega_{m})$ emerges at $T_{c}(N)$ and vanishes at $T=0$, i.e., it is is non-monotonic as a function of temperature. For $N<N_{cr}$, there are two regimes of qualitatively different behavior – in between $T_{c}(N)$ and ${\tilde{T}}_{c}(N)$, the pairing is induced by fermions with $\omega_{m}=\pi\pi T$, while at $T<{\tilde{T}}_{c}(N)$, fermions with all Matsubara frequencies contribute to the pairing. This last behavior is a conventional one, in the sense that it holds in a non-critical, BCS superconductor, while the behavior at ${\tilde{T}}_{c}(N)<T<T_{c}(N)$ is of non-BCS type as it is due to strong non-FL self-energy at all $\omega_{m}$ except for $\pm\pi T$. At small $\gamma$, $N_{cr}\approx 4/\gamma$, and the $T_{c}(N)$ and ${\tilde{T}}_{c}(N)$ lines remain close down to a very small $T\sim g(\gamma)^{1/\gamma}\ll g$. However, for $\gamma\leq 1$, the two lines separate already at $T\leq g$. We note in this regard that the range between ${\tilde{T}}_{c}(N)<T<T_{c}(N)$ exists for the physical case of $N=1$, and the lower boundary of this range rapidly decreases as $\gamma$ approaches the value equal to $1$. In other words, even for $N=1$, there exists an intermediate $T$ range where the pairing is induced by fermions with $\pm\pi T$, and would not exist if these fermions were excluded from the gap equation. The behavior of a system in this intermediate $T$ range at $N=1$ should be, at least qualitatively, the same as that at large $N$.

Below we study superconductivity induced by fermions with $\omega_{m}=\pm\pi T$ in some detail by solving non-linear gap equation at $T<T_{c}$. We first solve the gap equation in Matsubara frequencies and obtain the gap, the Free energy, and the specific heat, and then convert to real frequencies and obtain the spectral function and the DOS.

III Non-linear gap equation, $N>N_{cr}$

We begin with the case $N>N_{cr}$ when ${\tilde{T}}_{c}=0$, i.e. the pairing would be impossible if the self-energy did not vanish at $\omega_{m}=\pm\pi T$. The limit $N\gg 1$ can be treated analytically and we consider it in some detail below.

III.1 Non-linear gap equation in Matsubara frequencies.

The non-linear equation for the pairing vertex $\Phi^{*}(\omega_{m})$ along with the equation for the fermionic self-energy $\Sigma^{*}(\omega_{m})$ with the feedback from the pairing are given in (II). We recall that at large $N$ the pairing temperature $T_{c}(N)$ is obtained by solving the linearized equation for $\Phi^{*}(\omega_{m})$ for fermions with only two Matsubara frequencies $\omega_{m}=\pm\pi T$; the pairing vertex $\Phi^{*}(\omega_{m})$ for other $\omega_{m}$ is then expressed via $\Phi^{*}(\pi T)=\Phi^{*}(-\pi T)$. We assume and then verify that this holds also for $T<T_{c}$, i.e., that the non-linear gap equation can be approximated by restricting to $\omega_{m^{\prime}}=\pm\pi T$ in the r.h.s. of Eq. (II). Re-labeling $\Phi^{*}(\pi T)=\Phi^{(}-\pi T)=\Phi^{*}_{0},\Sigma^{*}(\pi T)=-\Sigma^{*}(-\pi T)=\Sigma^{*}_{0}$, and ${\tilde{\Sigma}}^{*}_{0}=\pi T+\Sigma^{(}\pi T)$ to shorten notations, we obtain from (II)

[TABLE]

The solution of (19) to leading order in $1/N$ is

[TABLE]

The superconducting gap $\Delta(\pm\pi T)=\Delta_{0}$ is

[TABLE]

The gap $\Delta_{0}$ vanishes both at $T=0$ and at $T=T_{c}$. In between, it is finite, but for any $T$, $D_{0}=\Delta_{0}/(\pi T)$ is small and at most of order $1/N^{1/2}$. In other words, the gap at $N\gg 1$ remains smaller than the temperature.

Solving next the set of Eliashberg equations for other $\omega_{m}\neq\pm\pi T$ we obtain at large $N$

[TABLE]

where $H(a,b)=\sum_{1}^{a}n^{-b}$ is a Harmonic number. We plot $\Phi^{*}(\omega_{m})$ and $\Sigma^{*}(\omega_{m})$ in Fig.9. Note that at $\omega_{m}\sim T$, $\Sigma^{*}(\omega_{m})\sim N\omega_{m}$, i.e., $\Sigma^{*}(\omega_{m})\approx{\tilde{\Sigma}}^{*}(\omega_{m})$.

At large $m$ (but still when $\Sigma^{*}(\omega_{m})\gg\omega_{m}$)

[TABLE]

Note that below $T_{c}$, the self-energy at all $\omega_{m}$, including $\omega_{m}=\pm\pi T$, behaves as $\Sigma^{*}(\omega_{m})\propto T^{1-\gamma}$, consistent with the scaling $\Sigma^{*}(\omega_{m})\propto(\omega_{m})^{1-\gamma}$. Still, the self-energy at $\pm\pi T$ is smaller in $1/N$ than $\Sigma^{*}(\omega_{m})$ at other Matsubara frequencies.

From (22) we have

[TABLE]

both at $m=O(1)$ and at $m\gg 1$. We see that at any $T<T_{c}$, $\Delta(\omega_{m})$ at any Matsubara frequency is parametrically smaller than $T$. Put it differently, $D(\omega_{m})=\Delta(\omega_{m})/\omega_{m}$ is small, of order $1/N^{3/2}$, at $m=O(1)$, and even smaller at larger $m$. We plot $\Delta(\omega_{m})$ and $D(\omega_{m})$ in Fig.10.

Taking $-iD_{0}$ as an estimate for small frequency limit of $D^{R}(\omega)\equiv\Delta^{R}((\omega)/\omega$ in real frequencies, we find that $D^{R}(\omega\to 0)$ tends to a finite imaginary value, i.e., at large $N$ is a gapless superconductivity in the sense that $\Delta^{R}(\omega)\propto i\omega$. For notational simplicity for all functions of real frequencies below we will drop the superscript “$R$”.foo Using then $N(\omega)=N_{0}\operatorname{Re}[1/\sqrt{1-D^{2}(\omega)}]$ for the DOS ($N_{0}$ is the normal state value), we find that the DOS at zero frequency $N(\omega=0)=N_{0}/\sqrt{1+D^{2}_{0}}\approx N_{0}\left(1-\frac{1}{2}D^{2}_{0}\right)$ is reduced below $T_{c}$ compared to the normal state value, but remains finite for any $T$, like it is expected in a gapless superconductor.

To verify this result and to get the full form of $N(\omega)$ we need to obtain $\Delta(\omega)$ as a function of a real frequency $\omega$. This is what we will do in Sec. III.3. Before that, we use the result for $D(\omega_{m})$ and obtain the Free energy $F_{sc}(T)$ and the specific heat $C(T)$ at $N>N_{cr}$.

III.1.1 The Free energy and the specific heat

The Free energy $F_{sc}$ and $\Delta F=F_{sc}-F_{n}$ are given by Eqs. (12)-(II) At large $N$, we keep only contributions which contain $D_{m}$, $D_{m^{\prime}}$ with $m,m^{\prime}=0,-1$. Contributions from $D_{m}$ with other $m$ are smaller in $1/N$, as we explicitly verified. Using that $\sum_{m}\frac{\operatorname{sgn}~{}m}{|\pi T\pm\omega_{m}|^{\gamma}}=0$, we obtain from (II)

[TABLE]

Varying $\delta F$ by $\Delta_{0}$, one reproduces Eq. (21). Substituting $D_{0}$ from (21) into (25), we obtain

[TABLE]

The specific heat variation between the superconducting and the normal state $\delta C_{v}=-T\partial^{2}\delta F/\partial T^{2}$ is

[TABLE]

where

[TABLE]

At $T\to 0$, $C_{\gamma}(0)\to 0$, i.e., $\delta C_{v}$ vanishes and $C_{v}$ recovers its normal state limiting behavior $C_{V}\propto T^{1-\gamma}$. At $T=T_{c}-0$, $C_{\gamma}=2\gamma^{2}$, i.e., the magnitude of the specific heat jump at $T_{c}$ is

[TABLE]

The specific heat in the normal state is obtained from (13). The first term in (13) gives the conventional free-fermion contribution to Free energy $F_{n,free}(T)=F_{n,free}(0)-N_{0}\pi^{2}T^{2}/3$. The second term gives

[TABLE]

At $T\sim T_{c}$, this second term is larger by $N$ than the free-fermion contribution. The calculation of the double sum in (30) requires care as one needs to extract the universal constant on top of formally ultra-violet divergent contribution, which actually is the factor in $F_{n,int}(0)$. To extract the universal constant, we note that the summation over $m-m^{\prime}$ can be done explicitly. The result is

[TABLE]

where, we remind, $H(m,\gamma)=\sum_{1}^{m}1/p^{\gamma}$ is the Harmonic number. For the remaining summation we use the Euler-Maclaurin formula

[TABLE]

where $B_{n}$ are Bernoulli numbers. The first term in the upper line in (32) contributes to $F_{n,int}(T=0)$, the second term determines the universal prefactor in the temperature-dependent piece in the Free energy. It is essential that the argument of the function under the sum is $m+1/2$ because this is how Matsubara frequency $\omega_{m}$ depends on $m$. Accordingly, we re-define $H(m+1/2,\gamma)=\sum_{1}^{m+1/2-1/2}1/p^{\gamma}$ and extend it to a function $H(x,\gamma)$ of a continuous variable $x$. Evaluating then the integral and the derivatives in the second line in (32) numerically, we obtain

[TABLE]

We plot $Q_{\gamma}$ in Fig.11.

Substituting the result into (30) and differentiating the free energy over $T$, we obtain

[TABLE]

The ratio of the specific heat jump to its value at $T=T_{c}+0$ is then

[TABLE]

We see that the relative jump of $C_{v}$ at $T_{c}$ is by $1/N^{2}$ smaller than in a BCS superconductor. In Fig.12 we plot $C_{v}(T)=C_{v,n}(T)+\delta C_{v}(T)$ in the full temperature range below $T_{c}$. At sufficiently small $T$, both $C_{v}$ and $C_{v,n}$ scales as $T^{1-\gamma}$.

III.2 Beyond leading order in $1/N$

We now go beyond the leading order in $1/N$. The goal here is to analyze how fermions with other $\omega_{m}$ affect the magnitudes of $\Phi(\pi T)=\Phi_{0}$ and $D(\pi T)=D_{0}$ at a small but finite temperature. We recall that at large $N$, $\Phi_{0}\approx(2/N)^{1/2}\pi T(T_{c}/T)^{\gamma}$ and $D_{0}\approx(2/N)^{1/2}$. We show that both $\Phi_{0}$ and $D_{0}$ increase as $N$ get smaller.

For the analysis to next order in $1/N$ we use the fact that $D_{0}\propto 1/N^{1/2}$, while for other Matsubara frequencies $D(\omega_{m})\propto 1/N^{3/2}$ (Eqs (21) and (24). Because $D$ appears in even powers in the equation for the self-energy in (II), the inclusion of these $D(\omega_{m})$ with $m\neq 0,-1$ would lead to corrections of at least of order $1/N^{2}$. To order $O(1/N)$ we then still have the same equation for ${\tilde{\Sigma}}^{*}_{0}$ as in (19). Expanding in this equation in two orders of $D^{2}_{0}\propto 1/N$ and setting $T\ll T_{c}$, we obtain

[TABLE]

The expansion to next order in $1/N$ in the equation for $\Phi_{0}$ requires more care as the leading term (the one kept in the first equation in (19)) is of order $1/N^{1/2}$, while other terms in the r.h.s. of the equation for $\Phi^{*}(\omega_{m})$ in (II) are of order $D(\omega_{m})\propto 1/N^{3/2}$, i.e., they contain only one additional power of $1/N$. These terms then should be kept in calculation to subleading order in $1/N$. Keeping these terms, we obtain from (II):

[TABLE]

Substituting $D(\omega_{m})$ from Eq. (22):

[TABLE]

we obtain

[TABLE]

where

[TABLE]

and, we remind, $H(m,\gamma)=\sum_{1}^{m}{1/n^{\gamma}}$ is a Harmonic number. We plot $W_{\gamma}$ in inset of Fig.14.

Solving (36) and (39) to order $1/N$ we obtain at low $T\ll T_{c}$

[TABLE]

The analysis at larger $T\leq T_{c}$ proceeds in the same way and we refrain from presenting the full formulas. In Fig. 13 we show $\Delta_{0}=\Delta(\pi T)$ as a function of $T/T_{c}$ for $\gamma=0.9$ and two different values of $N>N_{cr}$ ($N_{cr}\sim 1.3$ for $\gamma=0.9$). In both cases, $\Delta_{0}$ vanishes at $T=0$, but the slope of $\Delta_{0}(T)$ at small $T$ gets larger when $N$ decreases.

The result for $\Phi^{*}_{0}$ can be cast into $\Phi^{*}_{0}\approx(2/(N-N^{*}_{cr}))^{1/2}\pi T\left(\frac{T_{c}}{T}\right)^{\gamma}$ where $N^{*}_{cr}=3(W_{\gamma}-1)/2$ is some $\gamma$-dependent constant. Taking this approximate formula as an indication of the evolution of $\Phi^{*}_{0}$ with decreasing $N$, we find that $\Phi^{*}_{0}\propto T^{1-\gamma}/(N-N^{*}_{cr})^{1/2}$. At $N>N^{*}_{cr}(\gamma)$, $\Phi_{0}$ vanishes at $T=0$ (we recall that we consider $\gamma<1$), but $N=N^{*}_{cr}(\gamma)$ the slope of $\Phi^{*}_{0}(T)/T^{1-\gamma}$ (and of $\Delta_{0}$) diverges. This divergence is consistent with the $T=0$ analysis, which indicates that at $N<N_{cr}$, given by Eq. (17), the system has superconducting order at $T=0$. This will change the system behavior at small temperature and frequencies compared to what we found above. We emphasize that the increase of $\Phi^{*}_{0}(T\to 0)$ is due to the contribution from fermions with $|\omega_{m}|\neq\pi T$, which give rise to the $S$ term in $1/N$ correction. This means that, as $N$ get reduced, fermions with Matsubara frequencies other than $\pm\pi T$ become progressively more involved in the pairing.

The $N^{*}_{cr}(\gamma)=3(W_{\gamma}-1)/2$ is an approximate form of critical $N$ and does not have to coincide with the actual $N_{cr}(\gamma)$, given by Eq. (17). We plot both functions in Fig.14. Interestingly, $N^{*}_{cr}(\gamma)$ and $N_{cr}(\gamma)$ show quite similar variation with $\gamma$.

We next consider the solutions for the pairing vertex and the self-energy in real frequencies. This will allow up to compute the spectral function $A(\omega)$ and the DOS $N(\omega)$.

III.3 Non-linear gap equation in real frequencies

The transformation of Elishberg equations for electron-phonon interaction from Matsubara to real frequencies using spectral decomposition method and analytical continuation has been discussed in several publications Combescot (1995); Marsiglio et al. (1988); Marsiglio and Carbotte (1991); Karakozov et al. (1991). We extend these result to our case with $\chi(\Omega_{m})=(g/|\Omega_{m}|)^{\gamma}$. The details of the conversion procedure are presented in the Appendix. The conversion procedure requires special care by two reasons. First, if one simply replaces $\omega_{m}$ by $-i\omega$, the bosonic propagator $\chi(\omega_{m^{\prime}}+i\omega)$ will have a set of branch cuts in the complex $\omega$ plane, along $\omega=i\omega_{m}+b$, where $b$ is real. One then need to add additional terms to the r.h.s. of the equations for retarded functions $\Phi(\omega)$ and $\Sigma(\omega)$ to cancel these singularities and restore analyticity. Second, we again need to eliminate singular contributions from the terms with zero bosonic Matsubara frequency. This is done in the same way as in the calculations along the Matsubara axis. Namely, we introduce new functions $\Phi^{*}(\omega)$ and ${\tilde{\Sigma}}^{*}(\omega)$ related to $\Phi(\omega)$ and ${\tilde{\Sigma}}(\omega)=\omega+\Sigma(\omega)$ as

[TABLE]

where $Q_{\omega}$ is singular (see Eq. (46) below), but $\Phi^{*}(\omega)$ and ${\tilde{\Sigma}}^{*}(\omega)$ are free from singularities. The gap function $\Delta(\omega)=\omega\Phi(\omega)/{\tilde{\Sigma}}(\omega)$ is equally expressed in terms of $\Phi^{*}(\omega)$ and ${\tilde{\Sigma}}^{*}(\omega)$:

[TABLE]

The equations on $\Phi^{*}(\omega)$ and ${\tilde{\Sigma}}^{*}(\omega)$ are the same as on $\Phi(\omega)$ and ${\tilde{\Sigma}}(\omega)$, but with additional terms which cancel out divergent contribution from $\chi(0)$. We have (see Appendix for details)

[TABLE]

where

[TABLE]

and $\chi^{{}^{\prime\prime}}(x)=\operatorname{sgn}(x)\frac{g^{\gamma}}{|x|^{\gamma}}\sin\frac{\pi\gamma}{2}$. In these equations, the solution of the Eliashberg set in Matsubara frequencies, i.e., $\Phi^{*}(\omega_{m})$ and ${\tilde{\Sigma}}^{*}(\omega_{m})$ are considered as inputs. The first term in each of the two equations is obtained by just replacing $\omega_{m}$ by $-i\omega$, and the second one cancels out non-analyticities. The last piece in the second term cancels out the divergent contribution from $\chi(0)$. Note that the subtraction of the divergence at $x=0$ has to be done before extending the model to large $N$. The function $Q(\omega)$, which determines the relations between $\Phi^{*}(\omega)$ and ${\tilde{\Sigma}}^{*}(\omega)$ and the original $\Phi(\omega)$ and ${\tilde{\Sigma}}(\omega)$, Eqs. (42), is

[TABLE]

where

[TABLE]

Equivalently we can express $\Phi(\omega)$ and ${\tilde{\Sigma}}(\omega)$ via $\Phi^{*}(\omega)$ and ${\tilde{\Sigma}}^{*}(\omega)$ as

[TABLE]

where

[TABLE]

In Eqs. (45-49) the branch cut of the square root is defined along positive real axis.

At $\omega=0$ we have

[TABLE]

The first term in the formula for ${\tilde{\Sigma}}^{*}(0)$ vanishes by symmetry, after summing up the contributions from positive and negative $\omega_{m}$.

We first consider large $N$. We assume and then verify that in this case ${\tilde{\Sigma}}^{*}$ is parametrically larger than $\Phi^{*}$ not only along the Matsubara axis but also along the real axis. To leading order in $1/N$ we then have for the self-energy

[TABLE]

where

[TABLE]

We plot $S_{\gamma}$ in Fig.11

For $\Phi^{*}(0)$ we find from Eq. (III.3)

[TABLE]

Using the fact that at large N the dominant contribution to the Matsubara sum comes from $m=0,-1$ and substituting the expressions for $\Phi^{*}(\pm\pi T)$ and $\Sigma^{*}(\pm\pi T)$, we obtain

[TABLE]

Then $D_{0}=\Phi^{*}(0)/{\tilde{\Sigma}}^{*}(0)$ is

[TABLE]

and the DOS at zero frequency is

[TABLE]

This agrees, up to a prefactor, with the result that we obtained along the Matsubara axis, by assuming that $D(\pi T)$ is comparable with $D(\omega=0)$.

We emphasize that $N(0)$ differs from the normal state value $N_{0}$ at all $T<T_{c}$, including $T=0$, where we expect superconductivity to disappear. We will show below that the limit $\omega\to 0$ and $T\to 0$ has to be taken carefully, and at any non-zero $\omega$ the DOS indeed transforms into $N_{0}$ at $T\to 0$. Still, strictly at $\omega=0$, $N(0)<N_{0}$. This is similar, but indeed not identical, to behavior of $N(\omega)$ in an ideal BCS superconductor, where $N(0)=0$ for all $T$ up to $T_{c}$, while $N(\omega\neq 0)$ approaches $N_{0}$ at $T\to T_{c}$.

We next move to finite $\omega$. For $\Phi(\omega)$, the second term in (44) still scales as $\Phi(\omega)/N$ and can be neglected. Evaluating the first term by summing up contributions from $m=0,-1$ at which $\Phi(\omega_{m})/|{\tilde{\Sigma}}(\omega_{m})|$ is the largest at large $N$, we obtain

[TABLE]

where

[TABLE]

Note that in this large $N$ approximation $\Phi^{*}(\omega)$ is real and even in $\omega$.

Because $\Phi^{*}(\omega)$ is small in $1/N^{3/2}$, the self-energy at finite $\omega$ remains the same as in the normal state, up to $1/N^{3}$ corrections:

[TABLE]

where

[TABLE]

The first term in $F_{\Sigma}(x)$ is real, the second is imaginary. At large $x$ (i.e., at $\omega\gg\pi T$), $F_{\Sigma}(x)\approx(x^{1-\gamma}/(1-\gamma))e^{i\pi\gamma/2}$. We plot the scaling functions $F_{\Phi}(x)$, $\operatorname{Re}[F_{\Sigma}(x)]$, and $\operatorname{Im}[F_{\Sigma}(x)]$ in Fig. 15.

We see that Im $\left[F_{\Sigma}\left(\frac{\omega}{\pi T}\right)\right]$ changes sign as a function of frequency (and then Im $[\Sigma^{*}(\omega)]$ also changes sign). This sign change is necessary because $\Sigma^{*}(\pi T)=0$ and $\operatorname{Im}[\Sigma^{*}(\omega)]$ are related by Kramers-Kronig(KK) formula,

[TABLE]

and the vanishing of the integral in (61) is only possible if $\operatorname{Im}[\Sigma^{*}(\omega)]$ has different sign at small and large frequencies. We verified numerically that the KK relation is indeed satisfied, see Fig.16. We remind in this regard that $\Sigma^{*}$ is the self-energy without the thermal contribution. For the full self-energy $\operatorname{Im}[\Sigma(\omega)]$ indeed remains positive for all frequencies.

Substituting the results for $\Phi^{*}(\omega)$ and ${\tilde{\Sigma}}^{*}(\omega)$ into $D(\omega)=\Phi^{*}(\omega)/{\tilde{\Sigma}}^{(}\omega)$, we obtain

[TABLE]

where at $\omega\leq g$, when the bare $\omega$ term is smaller than $\Sigma^{*}(\omega)$, i.e., ${\tilde{\Sigma}}^{*}(\omega)\approx\Sigma^{*}(\omega)$,

[TABLE]

The DOS is

[TABLE]

We see that the magnitude of $N(\omega)/N_{0}-1=\frac{1}{2}\operatorname{Re}D^{2}(\omega)$ is determined by the overall temperature-dependent factor in (62) and depends on $T/T_{c}$ ratio. However, the frequency dependence of $D(\omega)$ and of the DOS is determined by $F_{D}(\omega/(\pi T))$, which for any given $\gamma$ is a universal function of $\omega/T$ and does not depend on $T/T_{c}$. This implies that the characteristic frequency, at which $N(\omega)$ deviates from $N_{0}$, is determined by the temperature rather than by the magnitude of the superconducting gap, as was the case for a BCS superconductor.

Because $F_{\Phi}(x)$ is real,

[TABLE]

At small $x=\omega/\pi T$, $\operatorname{Re}F_{\Sigma}(x)\propto x^{2}$ and $\operatorname{Im}F_{\Sigma}(x)$ is finite. Then $\operatorname{Re}F^{2}_{D}(x)$ is negative. At $x$ where $\operatorname{Im}F_{\Sigma}(x)$ changes sign, $\operatorname{Re}F_{\Sigma}(x)$ is finite, hence for this $x$, $\operatorname{Re}F^{2}_{D}(x)$ is positive. In between $\operatorname{Re}F^{2}_{D}(x)$ then necessary changes sign. This in turn implies that $N(\omega)<N_{0}$ at small $x$ and exceeds $N_{0}$ at larger $x$. At even larger $x\gg 1$, $N(\omega)$ approaches $N_{0}$. Then, for any $\gamma$, $N(\omega)$ has a dip at $\omega=0$ and a hump at a characteristic frequency set by temperature, rather than by the gap itself. This frequency then increases with increasing $T$, in qualitative difference with a BCS superconductor, in which the maximum in the DOS is located at $\omega=\Delta(T)$, and shifts to a lower frequency with increasing $T$ because $\Delta(T)$ gets smaller. We plot $\operatorname{Re}F^{2}_{D}(x)$ in Fig.17 for two different $\gamma$. The hump at $\omega\sim T$ is clearly visible. The position of the hump shifts to a lower frequency with increasing $\gamma$ but remains at a finite $\omega$ even at $\gamma=1$.

On a more careful look, we find that there is still a small difference in the behavior of the DOS between $\gamma<1/2$ and $\gamma>1/2$. Namely, at $\omega\gg T$, $\operatorname{Re}\Sigma^{*}(x)=\cos{\pi\gamma/2}(x^{1-\gamma}/(1-\gamma))$ and $\operatorname{Im}\Sigma^{*}(x)=\sin{\pi\gamma/2}(x^{1-\gamma}/(1-\gamma))$. As a result, $\operatorname{Re}F^{2}_{D}(x)\propto\cos{\pi\gamma}$, i.e., it is positive at $\gamma<1/2$ and negative at $\gamma>1/2$. This implies that for $\gamma>1/2$ $N(\omega)$ crosses $N_{0}$ twice at $\omega=O(T)$ because $(\operatorname{Im}F_{\Sigma}(x))^{2}$ is larger than $(\operatorname{Re}F_{\Sigma}(x))^{2}$ at both large and small frequencies. The second crossing at $x\sim 7$ is seen in Fig. 15 for $\gamma=0.9$. Digging further into this issue, we find that for $\gamma>1/2$, $N(\omega)$ crosses $N_{0}$ one more time, now at $\omega\sim g\gg T$, when the bare $\omega$ term in ${\tilde{\Sigma}}^{*}(\omega)$ becomes relevant, and at highest $\omega$ approaches $N_{0}$ from above. To see this, we extend the analysis of the DOS to $\omega\sim g$. The calculation is straightforward and we only cite the result: the difference $N(\omega)/N_{0}-1$ at $\omega\sim g$ is proportional to $\cos{\pi\gamma}+(1-\gamma)^{2}(\omega/g)^{2\gamma}+2(1-\gamma)(\omega/g)^{\gamma}\cos{\pi\gamma/2}=0$. Solving this equation for $\gamma>1/2$, we find the sign change of $N(\omega)/N_{0}-1$ at $\omega=\omega_{1}\sim g$. We show this in Fig. 18.

In Fig.19 and Fig.20 we show the results of the full numerical calculation of the temperature evaluation of the gap $\Delta(\omega)$ and the DOS $N(\omega)$ for two values of $\gamma$: $\gamma=0.3$ and $\gamma=0.9$. For $\gamma=0.9$, $N=6$ is above $N_{cr}\sim 1.3$. For $\gamma=0.3$ we show the results for $N=6$, which is below $N_{cr}\approx 9.6$ (the numerical analysis for $N>N_{cr}$ for such small $\gamma$ is challenging). For $N<N_{cr}$, the behavior similar to the one at large $N$ exists above the crossover temperature $T_{cross}(N)$ (see Sec. IV) and we show the results only in this $T$ range. The value of $T_{cross}(N=6)$ for $\gamma=0.3$ is only $0.01T_{c}$, so the range of $T>T_{cross}$ is rather wide. The gap $\Delta(\omega)$ is complex even at very small $\omega$, in contrast to the conventional BCS-like behavior where $\Delta(\omega)$ is almost real up to some frequency which is approximately equal to this real value. For the DOS we clearly see that there is a dip in $N(\omega)$ at small frequencies and a characteristic frequency $\omega_{p}$ at which $N(\omega)$ approaches $N_{0}$ is set by the temperature.

A remark is in order here. The $\int d\omega N(\omega)$, with $N(\omega)$ as in Fig.20 does have some $T$ dependence. At a first glance, this contradicts the requirement that the total number of particles is a conserved quantity. In fact, there is no contradiction. The reasoning is that the momentum integration in Eliashberg equations is performed assuming particle-hole symmetry, i.e., neglecting contributions from energies of order $\mu$. There are additional contributions to the DOS from energies of order $\mu$, both in the normal and the superconducting state. They are not equal, because $\mu$ changes between normal and superconducting states Chubukov et al. (2016). This additional contribution must be included to ensure particle conservation.

We next consider the spectral function $A(\omega)=-(1/\pi)\operatorname{Im}[G(k_{F},\omega)]$. In terms of original $\Phi(\omega)$ and ${\tilde{\Sigma}}(\omega)$, we have

[TABLE]

Expressing $\tilde{\Sigma}(\omega)$ and $\Phi(\omega)$ via $\tilde{\Sigma}^{*}(\omega)$ and $\Phi^{*}(\omega)$, Eq. (49), we find

[TABLE]

where

[TABLE]

To leading order in $1/P$, $A(\omega)\propto\frac{1}{P}\operatorname{Re}\left[\frac{1}{\sqrt{1-\left(\Phi^{*}(\omega)/{\tilde{\Sigma}}^{*}(\omega)\right)^{2}}}\right]\propto N(\omega)/N_{0}$, i.e., the spectral function has the same dependence on $\omega$ as the DOS. Accordingly, at a finite $T$, $A(\omega)$ is non-zero for any frequency, and the position of the maximum in $A(\omega)$ scales with $T$ and remains at a finite frequency at $T_{c}$ (Fig.21). Like we said, this behavior has been termed as “gap filling”. If $P$ is finite, either because the system is at some distance from a QCP, or we probe $A(\omega)$ for fermions not connected by momenta at which static $\chi$ diverges (like near-nodal fermions in the cuprates, if a pairing boson is an antiferromagnetic spin fluctuation), the behavior of $A(\omega)$ depends on the interplay between $P$ and the other term in $L(\omega)$ in (68). If $P$ is smaller, $A(\omega)$ is given by (67), (68). Substituting the expressions for $\Phi^{*}(\omega)$ and ${\tilde{\Sigma}}^{*}(\omega)$ we find that in this situation $A(\omega)$ is peaked at zero frequency, as if the system was in the normal state (see Fig.21).

The analysis beyond the leading order in $1/N$ proceeds in the same way as for Matsubara frequencies. As $N$ gets smaller, the maximum in the DOS becomes more pronounced, and, at the same time, the DOS at zero frequency, $N(0)$ gets smaller. These modifications get larger as $N$ decreases towards $N_{cr}$ and eventually qualitatively change the system behavior at $N<N_{cr}$ and $T<T_{cross}$, as we show in the next Section.

IV The case $N<N_{cr}$

At smaller $N<N_{cr}$ analytical solution is difficult to obtain because there is no obvious small parameter, so our discussion will be based on numerical results.

IV.1 Non-linear gap equation in Matsubara frequencies

In Fig.22 we show the results for $\Phi^{*}(\omega_{m})$ and $\Delta(\omega_{m})$. We see that now $\Phi^{*}(\pi T)=\Phi^{*}_{0}$ and $\Delta(\pi T)=\Delta_{0}$ tend to finite values at $T\to 0$, i.e., show a “conventional” superconducting behavior. Because at $T=0$ Matsubara frequency is a continuous variable, and it does not select between $\pm\pi T$ and other frequencies, the development of a finite gap at $T=0$ implies that at $N<N_{cr}$ and a finite $T$, there should exist a $T$ range in which all Matsubara frequencies equally contribute to the pairing, i.e., fermions with $\omega_{m}=\pm\pi T$ are no longer crucial to the pairing. This is consistent with our earlier result that at $N<N_{cr}$, the transition temperature remains finite even we exclude fermions with $\omega_{m}=\pm\pi T$ from Eliashberg equations (${\tilde{T}}_{c}(N)$ in Fig.8).

In Fig.23 we show $\Delta(\pi T)=\Delta_{0}$ as a function of $T$. The temperature dependence of $\Delta_{0}$ is still non-monotonic, i.e., as $T$ is reduced below $T_{c}$, $\Delta_{0}$ first increases, and then drops below a certain $T$, before reaching a finite value at $T\to 0$. As $N\to 1$, the maximum in $\Delta_{0}(T)$ becomes shallow. The frequency dependence of $\Phi^{*}(\omega_{m})$ and of $\Delta(\omega_{m})$ at a given $T$ is monotonic, with the maximum at $|\omega_{m}|=\pi T$.

In Fig.23 we compare the behavior of $\Delta_{0}(T)$ at $N>N_{cr}$ and $N<N_{cr}$. Near $T_{c}$, the behavior in the two cases is the same, but at low $T$ $\Delta_{0}(T)$ at $N>N_{cr}$ continue decreasing, while $\Delta_{0}(T)$ at $N<N_{cr}$ saturates. The temperature at which the two curve separate marks the crossover between the conventional behavior at low $T$ and the behavior, undistinguishable from the one at $N>N_{cr}$, at higher $T$. In the higher $T$ region, the pairing can still be viewed as induced by fermions with Matsubara frequencies $\pm\pi T$. The crossover line $T_{cross}(N)$ ends at $T=0$ at $N=N_{cr}$, just like ${\tilde{T}}_{c}(N)$, but these two temperatures are not directly proportional to each other.

IV.2 Non-linear gap equation in real frequencies

We used the same computational procedure as at large $N$ and obtained $\Phi^{*}(\omega)$, ${\tilde{\Sigma}}^{*}(\omega)$, and $\Delta(\omega)$ along the real frequency axis. We present the results in Fig.24. We again see the crossover in the system behavior around $T_{cross}(N)$. At smaller $T$, the behavior of the gap function is conventional in the sense that $\Delta^{\prime}(\omega=0)$ is finite and $\Delta^{{}^{\prime\prime}}(\omega)$ emerges only above a finite frequency $\omega_{0}\approx\Delta^{\prime}(0)$. At higher $T$, at small frequencies $\Delta^{{}^{\prime\prime}}(\omega)\propto\omega$ and $\Delta^{{}^{\prime}}(\omega)\propto\omega^{2}$, i.e., the systems displays gapless superconductivity. The self-energy $\Sigma^{*}(\omega)$ is strongly reduced below $\omega_{0}$ at small $T$ compared to that in the normal state, but almost recovers the normal state value in the regime of gapless superconductivity (Fig.24).

In Fig.25 we show the behavior of the DOS $N(\omega)$. We see qualitative change of the behavior between $T>T_{cross}(N)$ and $T<T_{cross}(N)$. At smaller $T$, the DOS is similar to that in a BCS superconductor: it has a sharp peak at $\omega\approx\Delta(0)$ and nearly vanishes below the peak frequency. At $T$ increases but remains smaller than $T_{cross}(N)$, the position of the maximum in $N(\omega)$ moves to a smaller frequency because $\Delta(0)$ get reduced, i.e., the gap in the DOS “closes in”. However, at higher $T>T_{cross}(N)$, the DOS becomes non-zero at all frequencies, and the position of its maximum moves to a higher frequency and remains finite at $T=T_{c}-0$, i.e., the gap in the DOS “fills in”. We plot the variation of the position of the maximum in $N(\omega)$ with $T$ on the right side of each DOS in Fig.25.

The spectral function $A(\omega)$ shows the similar crossover (Fig.26). In the limit when the thermal contribution is large, it shows the same behavior as $N(\omega)$. In the opposite limit, $A(\omega)$ at $T<T_{cross}(N)$ has two sharp peaks at frequencies close to $\pm\Delta(0)$, and at $T>T_{cross}(N)$ it has a single peak at $\omega=0$

V Discussion

In this work we analyzed the interplay between the tendency towards fermionic incoherence and the tendency towards pairing near a quantum-critical point in a metal. We used the $\gamma$ model of dynamical fermion-fermion interaction mediated by a critical boson with susceptibility $\chi(\Omega_{m})\propto(g/|\Omega_{m}|)^{\gamma}$. We extended the model to $SU(N)$ global symmetry and used $N$ as a parameter. At large $N$, the interaction in the pairing channel is smaller by $1/N$ than the one in the particle-hole channel, which gives rise to a fermionic incoherence. Earlier work by some of us and others Wang et al. (2016) found markedly different behavior at $T=0$ and at a finite $T$. Namely, the calculations at $T=0$ showed that superconductivity develops if $N$ is smaller than some $\gamma$-dependent $N_{cr}$, while at larger $\gamma$ the system remains in a NFL normal state. On the other hand, computations of the onset temperature for the pairing $T_{c}(N)$ showed that $T_{c}(N)$ remains finite at any $N$ and the line $T_{c}(N)$ by-passes $N_{cr}$ (Fig.2). The authors of Wang et al. (2016) argued that this discrepancy is due to the fact that Eliashberg equations for spin-singlet pairing contain fermionic self-energy without thermal contribution (the self-action term), and this self-energy is large (and has NFL) form for all frequencies except for $\omega_{m}=\pm\pi T$, at which it vanishes. The existence of a finite $T_{c}$ for any $N$ then follows from the fact that the pairing interaction between fermions with $\pi T$ and $-\pi T$ is not countered by the self-energy and opens the gap $\Delta$ at these two frequencies at $T=T_{c}-0$. A non-zero $\Delta(\pm\pi T)$ then induces the pairing gap for fermions with other Matsubara frequencies.

In this communication we extended the analysis of the pairing problem to $T<T_{c}(N)$ and solved the non-linear gap equation. We analyzed the large $N$ limit analytically and solved the gap equation at smaller $N$ numerically. We first obtained $\Delta(\omega_{m})$ along Matsubara axis and used it to compute the Free energy and the specific heat. We found that the specific heat jumps at $T_{c}$, but at large $N$ the relative magnitude of the jump $\Delta C(T_{c}-0)/C_{n}(T_{c})$ is smaller by a factor $1/N^{2}$ than in a BCS superconductor. The behavior of the specific heat below $T_{c}$ is also rather unconventional, as specific heat approaches normal state form at $T\to 0$.

We then solved the gap equation along the real axis, using $\Delta(\omega_{m})$ as input. We obtained $\Delta(\omega)$ and used it to compute the DOS $N(\omega)$ and the spectral function $A(\omega)$. In a conventional BCS-type superconductor $A(\omega)$ and $N(\omega)$ are peaked at the gap value $\Delta(T)$, and the peak position shifts to a smaller $\omega$ as temperature increases towards $T_{c}$ (the gap “closes in”). We found that at $N>N_{cr}$, the behavior is very different – the position of the maximum in $N(\omega)$ increases linearly with $T$ and remains finite at $T_{c}$. The DOS remains finite at all frequencies, including $\omega=0$. At small $T$, $N(\omega)$ at small $\omega$ is reduced compared to $N(\omega)=N_{0}$ in the normal state, it displays a pseudogap behavior. As $T$ increases towards $T_{c}$, pseudogap just “fills in”.

The form of the spectral function $A(\omega)$ depends on the strength of the thermal contribution. In our model, thermal contribution diverged at a QCP. In this limit, $A(\omega)$ has the same frequency dependence as the DOS $N(\omega)$. Away from a QCP, when the bosonic susceptibility $\chi(\Omega)$ is not singular at $\Omega=0$, thermal contribution does not diverge. In the limit when the thermal contribution is weak, $A(\omega)$ at when $N>N_{cr}$ has a single peak at $\omega=0$. At $N<N_{cr}$, it has two sharp peaks at $\omega\approx\pm\Delta(0)$, when $T<T_{cross}(N)$, and a single peak at $\omega=0$ when $T>T_{cross}(N)$.

The issue we didn’t discuss in this work is whether gap fluctuations (transverse and longitudinal) destroy long-range superconducting order in some $T$ range below $T_{c}(N)$. Eliashberg theory, which we used, neglects gap fluctuations. It is very likely that in some range below Eliashberg $T_{c}(N)$ long-range superconducting order gets destroyed, and the actual $T_{c,act}<T_{c}$. We note in this regard that in our theory, the transformation from “gap closing” to “gap filling” in the DOS and the spectral function at $T\sim T_{cross}(N)$ is due to the fact that at $T>T_{cross}(N)$, the feedback from the pairing on the fermionic self-energy is weak. This last result does not actually rely on the existence of long-range superconducting order. If fluctuations destroy superconducting phase coherence, the feedback will be further reduced, but we emphasize that the feedback on fermions as small above $T_{cross}(N)$ already within the Eliashberg theory. The same argument holds for the transformation from two peaks at a finite frequency in $A(\omega)$ to a single peak at $\omega=0$.

The transformation from “gap closing” behavior at small $T$ to “gap filling” behavior at $T\sim T_{c}$ has been observed in high-$T_{c}$ cuprates, in the DOS Fischer et al. (2007) and ARPES measurements of the spectral function in the antinodal region Fischer et al. (2007); Reber et al. (2012); Kondo et al. (2013); *Kaminski2; Kanigel et al. (2006); *kanigel2; *kanigel3; nor ; Damascelli et al. (2003); *shen2; Johnson et al. (2001); Kordyuk and Borisenko (2006); *kordyuk2; He et al. (2014); Peng et al. (2013). Symmetrized data of MDC ARPES measurements a along particular direction of ${\bf k}$ in the near-nodal region showed the transformation from two peaks at a finite frequency to a single peak at $\omega=0$ (this is termed as the appearance of the Fermi arc). These results are consistent with our microscopic analysis for the DOS and also for the spectral function, if we assume that the thermal contribution is stronger in the antinodal region than in the near-nodal region. The strength of thermal contribution scales with the static bosonic susceptibility $\chi^{\prime}(0)$. Static $\chi^{\prime}(0)$ is larger for antinodal fermions in, e.g., spin-fluctuation models Abanov et al. (2003); Monthoux et al. (2007); Sca , where the interaction is peaked at momentum at or near $(\pi,\pi)$.

A final remark. In our analysis we didn’t include the dependence of the pairing vertex and the gap on the angle along the Fermi surface, e.g., $\cos{2\theta}$ dependence for the $d-$wave gap in the cuprate superconductors. The $d-$wave form of the pairing gap does not affect our results for the DOS and $A(\omega)$ in the antinodal region, as there the gap can be approximated by the constant. We note in this regard that our results for the development of the dip with increasing $T$ in the DOS and in the spectral function at large $P$, (e.g., representative results in Fig.4 and right panel in Fig.5) are quire consistent with DOS and ARPES data in the cuprates. We associate our result for smaller $P$ (i.e., smaller thermal contribution) with the system behavior closer to the nodes. This association is valid if the pairing interaction in the cuprates is the strongest at momentum transfers connecting antinodal points and weaker at momentum transfer along the diagonals, like in spin-fluctuation scenario. Our results then show (left panel in Fig.5) that in the near-nodal regime, the two peaks, originally separated by $2\Delta$ at a particular $k$-point on the Fermi surface transform into a single peak at $\omega=0$ as $T$ increases. This effect is well known as the development of the Fermi arc.

The modeling of the angular dependence of the gap $\Delta(\theta)$ is needed for the analysis how the spectral function evolves as a function of $\theta$ at a given $T$. To obtain this dependence in our data, we added $\cos{2\theta}$ factor to $\Phi^{*}(\omega)$ and solved the Eliashberg equations at a given $P$, $T$, and $\gamma$. At high $T>T_{cross}$, the evolution is similar to the one in Fig.21. Namely, near the node $A(\omega)$ has a single maximum at $\omega=0$, while in the antinodal region $A(\omega)$ has a dip at $\omega=0$ and a shallow maximum, whose frequency scales with $T$. At $T<T_{cross}$, $A(\omega)$ has two weakly separated peaks in the nodal region and strongly separated peaks in the antinodal region (Fig.27) This behavior and the one in Fig. 26 reproduce ARPES data in Refs. Reber et al. (2012); Kanigel et al. (2006); Kondo et al. (2013); *Kaminski2; Damascelli et al. (2003); *shen2; Johnson et al. (2001); Kordyuk and Borisenko (2006); *kordyuk2; He et al. (2014); Peng et al. (2013).

Acknowledgements.

We thank D. Dessau, A. Millis, N. Prokofiev, S. Raghu, G. Torroba, and A. Yazdani for useful discussions. This work by Y. Wu and AVC was supported by the NSF DMR-1523036.

Appendix: Analytic continuation from Matsubara axis to real frequency axis

In this Appendix we show the derivation of Eq. (44) for the pairing vertex $\Phi^{*}(\omega)$ and the self-energy $\Sigma^{*}(\omega)$ along real frequency axis. We follow Ref. Marsiglio et al. (1988) and use spectral decomposition approach. To avoid misunderstanding, here we explicitly keep the factors $i$ for Matsubara frequencies, i.e. define the interaction as $\chi(i\Omega)=\frac{g^{\gamma}}{|\Omega|^{\gamma}}$. For a general complex number $z$ the retarder $\chi(z)$ is we have

[TABLE]

Along real frequency axis ($z=\omega+i\delta$) we have

[TABLE]

By Cauchi theorem, the susceptibility at arbitrary $z$ can be expressed via $\operatorname{Im}\chi(x)$ as

[TABLE]

Along real frequency axis this reduces to KK relation $\operatorname{Re}\chi(\omega)=(2/\pi){\cal P}\int_{0}^{\infty}dx\operatorname{Im}\chi(x)x/(x^{2}-\omega^{2})$. In the calculations of $\Phi(\omega_{m})$ and ${\tilde{\Sigma}}(\omega_{m})$ we used the susceptibility $\chi(i\omega_{m}-i\omega_{m^{\prime}})$, weighted with ${\Phi(\omega_{m^{\prime}})}/{\sqrt{{\tilde{\Sigma}}^{2}(\omega_{m^{\prime}})+\Phi^{2}(\omega_{m^{\prime}})}}$ and ${{\tilde{\Sigma}}(\omega_{m})}/{\sqrt{{\tilde{\Sigma}}^{2}(\omega_{m^{\prime}})+\Phi^{2}(\omega_{m^{\prime}})}}$, respectively, and summed up over $\omega_{m^{\prime}}$ (Eq. II). These expressions cannot be converted to real frequency $\omega$ by just replacing $z=i\omega_{m}$ by $z=\omega$ in $\chi(z-i\omega_{m^{\prime}})=(-g^{2}/(z-i\omega_{m^{\prime}})^{2})^{\gamma/2}$ because this $\chi(z)$ has branch cuts in a complex plane of $z$ along $z=i\omega_{m^{\prime}}+z_{0}$, where $z_{0}$ is a real variable (Ref. Marsiglio et al. (1988)). Because of this complication, we have to implement the full spectral decomposition procedure. Namely, we depart from Eliashberg equations along Matsubara axis and use spectral representation to express $G(i\omega_{m},\bm{k})$ and $\chi(i\omega_{m}-i\omega_{m^{\prime}})$ via $\operatorname{Im}G(x,\bm{k})$ and $\operatorname{Im}\chi(x)$ along real axis as

[TABLE]

where ”R” stands for ”retarded”. We then explicitly sum over $\omega_{m^{\prime}}$ and integrate over ${\bf k}$ and obtain the expressions for ${\tilde{\Sigma}}(i\omega_{m})$ and $\Phi(i\omega_{m})$, in which the dependence on $\omega_{m}$ is only via $1/(i\omega_{m}-x-y)$. This form can be straightforwardly continued analytically to real frequency by just replacing $i\omega_{m}$ by $\omega+i0^{+}$.

For compactness, we do the calculations in Nambu formalism, in which one operates with the matrix Green’s function ${\hat{G}}(i\omega_{m},\bm{k})$ and treats both $\Sigma(i\omega_{m})$ and $\Phi(i\omega_{m})$) as elements of matrix self-energy ${\hat{\Sigma}}(i\omega_{m})$. The Eliashberg equation in Nambu formalism is

[TABLE]

where ${\hat{\tau}_{3}}$ is a Pauli matrix. $\hat{\Sigma}=\Sigma{\hat{\tau}_{0}}-\Phi\hat{\tau}_{1}$, and the matrix $\hat{G}(i\omega_{m},\bm{k})=-(i\omega_{m}-\hat{\Sigma}(i\omega_{m}))^{-1}$. The diagonal and off-diagonal elements of $\hat{G}(i\omega_{m},\bm{k})$ are conventional normal and anomalous Green’s functions.

Substituting the spectral representation (72) into (73) and performing the summation over $\omega_{m}$, we obtain

[TABLE]

Replacing $i\omega_{n}$ with $\omega+i0^{+}$ we obtain the self-energy at real frequencies

[TABLE]

We next express $\operatorname{Im}{\hat{G}}^{R}(x,\bm{k})/(\omega-x-y+i0^{+})$ via the full ${\hat{G}}^{R}(x,\bm{k})$ as

[TABLE]

and integrate over $x$ by closing the integration contour over the upper half-plane of complex $x$ Because ${\hat{G}}^{R}(x,\bm{k})$ is analytic in the upper half plane, the poles are at $x=\omega-y+i0^{+}$ and $x=i(2n+1)\pi T$ [these are the poles coming from $n_{F}(x)$]. Using the residue theorem, we find

[TABLE]

Using now $(1/\pi)\int dy\operatorname{Im}\chi(y)/(\omega-i\omega_{m}-y)=-\chi(\omega-i\omega_{m})$, we finally obtain

[TABLE]

Let’s now spit this matrix equation into the equations for the pairing vertex $\Phi(\omega)$ and conventional (non-anomalous) self-energy $\Sigma(\omega)$. Expressing $\hat{\Sigma}(\omega)$ as

[TABLE]

and substituting into the Dyson equation, we obtain

[TABLE]

where $\xi_{\bm{k}}$ is the fermionic dispersion. Expressing next $\int d^{2}k/(2\pi)^{2}=N_{0}\int d\xi_{\bm{k}}$, where $N_{0}$ is the DOS in the normal state, and integrating over $\xi_{\bm{k}}$, we obtain

[TABLE]

where ${\tilde{\Sigma}}(\omega)=\omega+\Sigma(\omega)$. Absorbing the density of states $N_{0}$ into $\chi$ and splitting $\hat{\Sigma}(\omega)$ into normal and anomalous components, we obtain

[TABLE]

where

[TABLE]

At a finite $T$ and small $y$, $n_{B}(y)\approx T/y$. At a QCP $\operatorname{Im}\chi(y)n_{B}(y)$ then scales as $T/|y|^{1+\gamma}$, and integrals over $dy$ in (82) diverge. The divergence, however, cancels out in the ratio $\Phi(\omega)/{\tilde{\Sigma}}(\omega)=\Delta(\omega)/\omega$ and can be formally eliminated by introducing new $\Phi^{*}(\omega)$ and $\Sigma^{*}(\omega)$ in which the divergent pieces are subtracted:

[TABLE]

Comparing (82) and (84) we see that

[TABLE]

where

[TABLE]

The ratio $\Phi^{*}(\omega)/{\tilde{\Sigma}}^{*}(\omega)$ is the same as $\Phi(\omega)/{\tilde{\Sigma}}(\omega)$, i.e., the gap function $\Delta(\omega)$ can be equally expressed via non-singular $\Phi^{*}(\omega)$ and ${\tilde{\Sigma}}^{*}(\omega)$. Furthermore, a little experimentation shows that $S_{\Phi}(\omega)$ and $S_{\Sigma}(\omega)$, given by (83), can be equally expressed via $\Phi^{*}(\omega)$ and ${\tilde{\Sigma}}^{*}(\omega)$, as

[TABLE]

By the same reason, $\Phi(\omega)$ and ${\tilde{\Sigma}}(\omega)$ can be expressed via $\Phi^{*}(\omega)$ and ${\tilde{\Sigma}}^{*}(\omega)$ in a manner similar to Eq. (85):

[TABLE]

where

[TABLE]

Equations (84) are free from divergencies and can be readily extended to $N\neq 1$, as we did in the main text.

Eqs. (84) have been solved numerically by iterations. For practical purposes, we found that in some cases the convergence is faster if we do calculations in two steps: first evaluate intermediate $\Phi^{**}$ and ${\tilde{\Sigma}}^{**}$, related to $\Phi$ and ${\tilde{\Sigma}}$ as in (85), but with $P=\int_{-\delta}^{\delta}dy\operatorname{Im}\chi(y)\frac{T}{y}$, where $\delta$ is some finite number, and then compute $\Phi^{*}$ and $\Sigma^{*}$ by adding the rest of the integral in $P$. The best convergence is achieved by adjusting the value of $\delta$.