Implicit renormalization approach to the problem of Cooper instability

TL;DR

This paper introduces an implicit renormalization method to accurately predict superconducting transition temperature and gap function from high-temperature data, bypassing traditional eigenvalue approaches and complex vertex calculations.

Contribution

It proposes an alternative eigenvalue problem within the implicit renormalization framework that improves the prediction of $T_c$ and the gap function without explicit vertex computations.

Findings

Accurate $T_c$ prediction from high-temperature data.

Eigenvalue problem formulation within implicit renormalization.

Solution via diagrammatic Monte Carlo without matrix inversion.

Abstract

In the vast majority of cases, superconducting transition takes place at exponentially low temperature $T_c$ out of the Fermi liquid regime. We discuss the problem of determining $T_c$ from known system properties at temperatures $T \gg T_c$, and stress that this cannot be done reliably by following the standard protocol of solving for the largest eigenvalue of the original gap-function equation. However, within the implicit renormalization approach, the gap-function equation can be used to formulate an alternative eigenvalue problem, solving which leads to an accurate prediction for both $T_c$ and the gap function immediately below $T_c$. With the diagrammatic Monte Carlo techniques, this eigenvalue problem can be solved without invoking the matrix inversion or even explicitly calculating the four-point vertex function.

Figures

Equations

$$V_{\rm RS}=\left\{\begin{array}[]{l}0,\quad\mbox{if}\;\>|\omega_{n}|>E_{c},\;\mbox{or}\;|\omega_{m}|>E_{c},\\ \\ g-gf\theta(\Omega\!-\!|\omega_{n}|)\theta(\Omega\!-\!|\omega_{m}|)\;\;\mbox{otherwise},\end{array}\right.$$

$$V_{\mathrm{eff}}(\omega )=g\left(1-\frac{{\Omega }_a^2}{{\omega }^2+{\Omega }^2}\right)=g\left(\frac{{\omega }^2+{\Omega }_1^2}{{\omega }^2+{\Omega }^2}\right),$$

$$\lambda (T){\Delta }_{\mathbf{k},{\omega }_n}=-T\sum _m\int \frac{d\mathbf{p}}{(2\pi {)}^d}{\Gamma }_{\mathbf{p},{\omega }_m}^{\mathbf{k},{\omega }_n}{G}_{\mathbf{p},{\omega }_m}^{(2)}{\Delta }_{\mathbf{p},{\omega }_m}.$$

$$\lambda (T){\Delta }_{{\omega }_n}=-\pi T\sum _{m}V({\omega }_m,{\omega }_n)\frac{{\Delta }_{{\omega }_m}}{|{\omega }_m|}.$$

$$\Delta ({\omega }_n)=a\theta (E_c-|{\omega }_n|)+(b-a)\theta (\Omega -|{\omega }_n|).$$

$$\ell =\pi T\sum _{E_c>|{\omega }_m|>\Omega }\frac{1}{|{\omega }_m|},L=\pi T\sum _{|{\omega }_m|<\Omega }\frac{1}{|{\omega }_m|}.$$

$$L=\ln [\Omega /\alpha T],\ell =\ln [E_c/\Omega ].$$

$$\displaystyle\left\{{\begin{array}[]{*{20}{l}}\lambda a=-g\ell a-gLb,\\ \lambda b=-g(1-f)Lb-gla.\end{array}}\right.$$

$${\lambda }_{\max }(T)=\frac{g}{2}\left[\sqrt{(L(1-f)+\ell {)}^2+4fL\ell }-L(1-f)-\ell \right].$$

$$f<1/(1+g\ell ),$$

$$L=\frac{1+g\ell }{g}\frac{1}{g\ell f-(1-f)}.$$

$$\lambda {\Delta }_n=-g\pi T\sum _m^{|{\omega }_m|<E_c}\frac{({\omega }_n-{\omega }_m{)}^2+{\Omega }_1^2}{({\omega }_n-{\omega }_m{)}^2+{\Omega }^2}\frac{{\Delta }_m}{|{\omega }_m|}.$$

$$\lambda \Delta (\omega )=-g{\int }_{\alpha T<|x|<E_c}\frac{dx\Delta (x)}{2|x|}\frac{(\omega -x{)}^2+{\Omega }_1^2}{(\omega -x{)}^2+{\Omega }^2}.$$

$$\Delta (\omega )=-g{\int }_{\alpha T<|x|<E_c}\frac{dx\Delta (x)}{2|x|}\frac{(\omega -x{)}^2}{(\omega -x{)}^2+{\Omega }^2}.$$

$$\Delta (0)=-g{\int }_{\alpha T<x<E_c}\frac{dx\Delta (x)}{x}\frac{x^2}{x^2+{\Omega }^2}.$$

$$\Delta (\omega )$$

$$\Delta (\omega )\approx \Delta (0)\left(1-gL\frac{{\omega }^2}{{\omega }^2+{\Omega }^2}\right).$$

$$1=-g\ell +g^2\ell L.$$

$$L=\ln \frac{\Omega }{\alpha T}=\frac{1+g\ell }{g^2\ell }.$$

$$L=\ln \frac{\Omega }{\alpha T}=\frac{1+g{\beta }_1}{g^2{\beta }_2-g^3{\beta }_3},$$

$${\beta }_1$$

$${\beta }_2$$

$${\beta }_3$$

$$\Phi (x)$$

$$\Delta (\omega )=\Delta (0)\left(1-\delta \frac{{\omega }^2}{{\omega }^2+{\Omega }^2}\right).$$

$$\delta =g{\Omega }^4{\int }_{\alpha T}^{E_c}\frac{dx}{x}\frac{{\Omega }^2-3x^2}{(x^2+\Omega {)}^3}\left(1-\delta \frac{x^2}{x^2+{\Omega }^2}\right).$$

$$1=-g{\int }_0^{E_0}dx\frac{x}{x^2+{\Omega }^2}\left(1-\delta \frac{x^2}{x^2+{\Omega }^2}\right).$$

$$\delta =\frac{1+g{\beta }_1}{g{\beta }_2}.$$

$$L=\ln \left[\frac{\Omega }{\alpha T_c}\right]=c_1+\frac{(1+g{\beta }_1)(1-c_{2}g)}{g^2{\beta }_2}$$

$$c_1={\gamma }_1+4{\gamma }_2+\ell -{\beta }_1,$$

$$c_2=3{\gamma }_2-4{\gamma }_3,$$

$${\gamma }_k={\int }_0^{E_c}\frac{dxx{\Omega }^{2k}}{(x^2+{\Omega }^2{)}^{1+k}}=\frac{1}{2k}\left[1-\frac{{\Omega }^{2k}}{(E_c^2+{\Omega }^2{)}^k}\right].$$

Full text

I Introduction

A conventional Bardeen-Cooper-Schrieffer (BCS) theory of $s$-wave superconductivity and its extensions to other pairing symmetries, and to strong coupling, assume that the pairing interaction is attractive in at least one pairing channel. In the BCS theory, the attraction comes from phonon exchange, in theories of non-ordinary $s$-wave superconductivity (e.g., $d$-wave superconductivity in the cuprates, $s^{+-}$ superconductivity in Fe-based metals, etc.), the attractive interaction between fermions is believed to originate from screened Coulomb interaction between electrons. In itinerant models of the electronic pairing, the attraction comes from an exchange of collective bosonic excitations in either spin channel (the spin-fluctuation exchange for superconductivity near the onset of magnetism Monthoux et al. (1991); Scalapino (2012); Abanov et al. (2003); Mazin et al. (2008); Kuroki et al. (2008); Chubukov et al. (2008)) or in charge channel (e.g., the exchange of nematic fluctuations near the onset of a nematic order Fradkin et al. (2010); Lederer et al. (2015, 2017); Berg et al. (2019); Klein and Chubukov (2018)).

Recent interest in $s$-wave superconductivity in SrTiO3  Schooley et al. (1964, 1965); Lin et al. (2014); Valentinis et al. (2017), Pb1-xTlxTe  Chernik and Lykov (1981), half-Heusler compounds Nakajima et al. (2015), and single-crystal Bi Prakash et al. (2017) re-ignited the discussion of another aspect of the pairing problem: the interplay between a weaker electron-phonon interaction and a stronger electron-electron repulsion Gurevich et al. (1962); Schooley et al. (1964, 1965); Chernik and Lykov (1981); Takada (1980); Ikeda et al. (1992); Grimaldi et al. (1995); Mahan (2000); Lin et al. (2014); Nakajima et al. (2015); Prakash et al. (2017); Edge et al. (2015); Ruhman and Lee (2016); Gor’kov (2016, 2017); Ruhman and Lee (2016, 2017); Lee (2015); Rademaker et al. (2016); Zhou and Millis (2016, 2017); Trevisan et al. (2018); Savary et al. (2017); Rowley et al. (2018); Coak et al. (2018); Wölfle and Balatsky (2018); Sadovskii (2018a, b); Aperis and Oppeneer (2018); Schrodi et al. (2018). BCS theory (and its finite-coupling version, known as Eliashberg theory Eliashberg (1960); Migdal (1958); Scalapino (1969); Carbotte (1990); Marsiglio and Carbotte (2008)) neglects electron-electron repulsion and considers only the phonon-mediated attraction. It has been argued long ago Tolmachov and Tiablikov (1958); Morel and Anderson (1962); Scalapino et al. (1966); McMillan (1968); Scalapino (1969); Carbotte (1990) that in a situation, when the Debye frequency $\Omega_{D}$ is much smaller than the Fermi energy $E_{F}$ (or, more precisely, the characteristic energy $E_{c}\leq E_{F}$, up to which one can expand the dispersion linearly near the Fermi level), the repulsive Coulomb interaction gets logarithmically suppressed by fermions with energies between $E_{c}$ and $\Omega_{D}$, and electron-phonon attraction, emerging at energies below $\Omega_{D}$, wins over the reduced Coulomb repulsion. The actual situation is more tricky, however, because the net interaction—the sum of the Coulomb repulsion and the electron-phonon interaction—is repulsive, and remains repulsive at all energies. The key effect of the electron-phonon interaction is that it makes the effective pairing interaction $V_{\rm eff}(\omega)$ dependent on the transferred frequency, interpolating between a smaller value at $\omega\to 0$ and a larger value at $\omega>\Omega_{D}$ Gurevich et al. (1962); Scalapino et al. (1966); McMillan (1968); Takada (1980); Coleman (2015); Ruhman and Lee (2016, 2017). A simplified model of such an interaction has been considered by Rietschel and Sham Rietschel and Sham (1983); Coleman (2015). They replaced the actual frequency-dependent $V_{\rm eff}(\omega=\omega_{m}-\omega_{n})$ along the Matsubara axis by a step-like $V_{\rm RS}(\omega_{m},\omega_{n})$ with separable dependence on $\omega_{m}$ and $\omega_{n}$:

[TABLE]

where $g$ is positive and $0<f<1$ [i.e., the interaction (1) is of repulsive character]. This $V_{\rm RS}$ has three values: $g(1-f)$ at small frequencies, $g$ at larger frequencies, and zero at very high frequencies. In the limit when $l=\ln[E_{c}/\Omega]$ is large, the analysis of the linearized gap equation shows Rietschel and Sham (1983); Coleman (2015) that $T_{c}$ is finite when $f>1/(1+gl)$ (see below). At large enough $l$ this holds for any $f>0$, i.e., for any interaction, which gets reduced below $\Omega$. The contribution of the “average” repulsion $g(1-f/2)$ to the gap equation is eliminated by sign change of the gap function $\Delta(\omega_{n})$ between $\omega_{n}<\Omega$ and $\Omega<\omega_{n}<E_{c}$, much like the on-site Hubbard repulsion gets eliminated from the gap equation for $s^{+-}$ superconductivity Maiti and Chubukov (2014).

This paper has two goals. First, we want to analyze superconductivity for more realistic repulsive interaction. Several recent studies of superconductivity in SrTiO3 and Bi argued that $V_{\rm eff}(q,\omega)$ can be viewed as the sum of a screened Coulomb interaction and an interaction with a gapped boson, dressed by the Coulomb potential, where a boson is a hybridized mode between a longitudinal phonon and a plasmon Mahan (2000); Ruhman and Lee (2016, 2017); Gastiasoro et al. (2019). We focus on the frequency dependence of $V_{\rm eff}(q,\omega)$ and neglect its momentum dependence. Specifically, we consider

[TABLE]

where $\omega$ is a running transferred Matsubara frequency, $\Omega$ is the frequency of a bosonic mode, and $\Omega_{a}\leq\Omega$, i.e., $\Omega^{2}_{1}=\Omega^{2}-\Omega^{2}_{a}<\Omega$. In the limiting case $\Omega_{a}=\Omega$, i.e., $\Omega_{1}=0$, the model reduces to the modified Bardeen-Pines model Bardeen and Pines (1955) (with a gapped boson at frequency $\Omega$ instead of an acoustic mode, as in the original Bardeen-Pines model).

The pairing interaction $V_{\rm eff}(\omega)$ is similar to $V_{\rm RS}(\omega)$ in the sense that it reduces to a larger repulsion $g$ at large frequencies and to a smaller repulsion $g(\Omega_{1}/\Omega)^{2}<g$ at small frequencies. However, in distinction to $V_{\rm RS}(\omega)$, the $V_{\rm eff}(\omega)$ form smoothly interpolates between the two limits, and is a function of a single variable $\omega=\omega_{m}-\omega_{n}$, rather than a separable function of $\omega_{m}$ and $\omega_{n}$.

We analyze superconductivity in this model analytically, in the weak coupling limit of small $g$. We show that the results are qualitatively similar to those of the Rietschel-Sham model in that at $\Omega_{1}=0$, $T_{c}$ is non-zero for arbitrary weak $g$, and for $0<\Omega_{1}<\Omega$, there is a threshold on $g$, below which superconductivity does not develop. However, the threshold value and the value of $T_{c}$ above the threshold are different from those in the Rietschel-Sham model. The absence of a threshold on $g$ at $\Omega_{1}=0$, i.e., at $\Omega_{a}=\Omega$ can be understood as the consequence of the fact that this is a boundary between a repulsive and an attractive interaction: at infinitesimally larger $\Omega_{a}$ the interaction becomes attractive at the smallest frequencies, in which case $T_{c}$ is finite at any $g$, like in BCS theory.

In most of realistic cases, an analytical solution is not possible and one has to rely on a numerical procedure of determining $T_{c}$ and the gap function immediately below $T_{c}$. This brings us to the second goal of our work—to set up the computational protocol to obtain critical parameters in the generic case of frequency and momentum dependent interaction. In general, one has to solve for the largest eigenvalue $\lambda_{\rm max}$ of the gap-function equation and obtain $T_{c}$ from the condition $\lambda_{\rm max}=1$  Scalapino (1969); Gladstone et al. (1969). In practice, $T_{c}$ at weak coupling is small, and to reach it, one needs to consider a very large number of Matsubara points (and momenta points near the Fermi surface). This creates a serious computational challenge. The conventional recipe in this situation would be to assume that at a temperature, at which a numeric simulation is done, the gap function is already saturated to its value in the superconducting state immediately below $T_{c}$. Below we will call this a critical gap function. Under this assumption, the functional form, and, correspondingly, the flow of $\lambda_{\rm max}$ with $T$ is logarithmical, allowing one to approximate $\lambda_{\rm max}$ as $\lambda_{\rm max}=g\rho_{F}\ln[\Omega/T]+{\rm const}$, where $\rho_{F}$ is the density of states per spin component. Using this relation, one can extrapolate $\lambda_{\rm max}(T)$ from $T\gg T_{c}$, where $\lambda_{\rm max}$ can be obtained with a manageable number of frequency and momentum points, to a much smaller $T$, and obtain $T_{c}$ from the condition $\lambda_{\rm max}(T_{c})=1$.

We argue that, while this approach works well for the case of an attractive potential, when the gap function does not change sign as a function of frequency, it fails for the case of a frequency dependent repulsive interaction, like $V_{\rm RS}(\omega_{m},\omega_{n})$ or $V_{\rm eff}(\omega)$. The key reason is that for a repulsive interaction, the ratio of the gap function in the frequency range where it is positive, and the one where it is negative, by itself depends on the logarithm of temperature, and this additional logarithmical dependence cannot be neglected in the extrapolation procedure. This leads to a non-linear dependence of $\lambda_{\rm max}$ on $\ln[\Omega/T]$, and makes the extrapolation from $T\gg T_{c}$ unreliable for determining $T_{c}$. We show this explicitly for both the Rietschel-Sham and $V_{\rm eff}(\omega)$ models.

We introduce the new protocol, which overcomes this complication. We call it an implicit renormalization approach. Specifically, we re-formulate the eigenvalue problem in such a way that the new largest eigenvalue, ${\bar{\lambda}}_{\rm max}(T)$, which still obeys ${\bar{\lambda}}_{\rm max}(T_{c})=1$, remains linear in $\ln[\Omega/T]$ in the whole $T$ range between $\Omega$ and much smaller $T_{c}$. We show that this allows one to determine $T_{c}$ with high precision by the extrapolation from higher $T$. Furthermore, the eigenvector of the implicit renormalization protocol is straightforwardly related to the critical gap function. After re-weighting its low-frequency part with the factor ${\bar{\lambda}}_{\rm max}(T)$, the former becomes equal to the latter. We apply the implicit renormalization method to both the Rietschel-Sham model and the model with $V_{\rm eff}(\omega)$ and in both cases find that it works with a remarkable accuracy.

The structure of the paper is the following. In the next section we present the analytical solutions for $T_{c}$, first for the exactly solvable Rietschel-Sham model and then for the model with more realistic $V_{\rm eff}(\omega)$. In Sec. III we discuss the extrapolation problem and introduce the implicit renormalization approach to determine $T_{c}$ in both models. In Sec. IV, we explain the method in more general terms and how to apply it to metals, when the irreducible interaction in the Cooper channel is given by a four-point vertex $\Gamma$, which depends on both frequency and momentum deviations from the Fermi surface. We argue that even a static repulsive interaction, which is weaker near the Fermi surface than away from it, may give rise to a finite $T_{c}$. In this Section we also discuss the Diagrammatic Monte Carlo (DiagMC) method and explain that to successfully apply the implicit renormalization scheme one does not need to know $\Gamma$ explicitly—all computed objects are no more complex than the single-particle Green’s function while all integrals and summations over the diagrammatic space are performed stochastically. In Sec. V, we provide further discussion and present conclusions.

II Analytical solution for $T_{c}$ in models with frequency-dependent repulsive interaction

The eigenvalue problem for the gap function $\Delta_{{\mathbf{k}},n}$ in $d$ dimensions reads [see Fig. 1(a)-(b); we use units such that $\hbar=1$ and $k_{B}=1$]:

[TABLE]

Here $\Gamma$ is a four-point vertex with zero incoming momentum and frequency, irreducible in the particle-particle (Cooper) channel, $\omega_{n}=\pi T(2n+1)$ and $\omega_{m}=\pi T(2m+1)$ are fermionic Matsubara frequencies, $G^{(2)}$ is the product of two single particle Green’s functions, $G^{(2)}_{{\mathbf{p}},m}\equiv G_{{\mathbf{p}},m}G_{-{\mathbf{p}},-m}$. To simplify the discussion, we consider $s$-wave pairing and assume that $\Gamma^{{\mathbf{k}},n}_{{\mathbf{p}},m}$ depends on frequencies $\omega_{n}$ and $\omega_{m}$, but does not depend on momenta; for the Green’s function we take the form $G_{{\mathbf{p}},m}=1/(i\omega_{m}-\xi_{p})$ and assume that the quasiparticle dispersion $\xi_{p}$ can be linearized around the Fermi surface for all $|\xi_{p}|\leq E_{c}$. Under these assumptions, the momentum integration can be carried out exactly, and the eigenvalue problem reduces to

[TABLE]

We first consider the exactly solvable Rietschel-Sham model with $V(\omega_{m},\omega_{n})=V_{\rm RS}(\omega_{m},\omega_{n})$ and then discuss analytical approach to the model with frequency-dependent repulsive interaction $V_{\rm eff}(\omega_{n}-\omega_{m})$.

II.1 Rietschel-Sham model

The pairing interaction in the Rietschel-Sham model, $V_{\rm RS}(\omega_{m},\omega_{n})$, is given by Eq. (1). We recall that this interaction has a step-like form and is separable between $\omega_{m}$ and $\omega_{n}$. The step-like $V_{\rm RS}$ equals $g(1-f)$ at small frequencies, $|\omega_{m}|,|\omega_{n}|<\Omega$ and equals $g$ at large frequencies, $\Omega<|\omega_{m}|<E_{c},~{}\Omega<|\omega_{n}|<E_{c}$.

For such pairing potential, the gap function $\Delta(\omega_{n})$ also displays a step-like behavior and can be parameterized as

[TABLE]

The eigenvalue problem is then reduced to solving the quadratic equation for the ratio $a/b$. It is convenient to define

[TABLE]

To logarithmic accuracy,

[TABLE]

Here $L$ is a conventional Cooper logarithm, for which $\alpha=0.882$ is the factor that one needs to add to convert the discrete sum $2\pi T\sum_{\pi T}^{\Omega}\omega_{m}^{-1}$ into the $\int_{\alpha T}^{\Omega}d\omega/\omega$, and $\ell$ (sometimes called Tyablikov-McMillan logarithm) accounts for the reduction of the repulsive interaction at energies between $E_{c}$ and $\Omega$.

Using these notations and assuming that both $L$ and $\ell$ are large, one can cast the eigenvalue problem into a compact form

[TABLE]

The transition temperature is determined by the condition $\lambda_{\rm max}(T_{c})=1$, where $\lambda_{\rm max}$ is the largest eigenvalue of (10). Solving for $\lambda_{\rm max}$ we obtain

[TABLE]

We see that $\lambda_{\rm max}$ is positive and monotonically increases with decreasing temperature, as long as $f$ is non-zero, i.e., as long as the interaction varies between small and large frequencies. However, we also see that $\lambda_{\rm max}(T)$ is not simply proportional to $L$, i.e., the increase of $\lambda_{\rm max}(T)$ with decreasing $T$ is not simply logarithmical. Moreover, at $T=0$, when $L\to\infty$, $\lambda_{\rm max}$ saturates at $\lambda_{\rm max}(T=0)=g\ell f/(1-f)$. When this limiting value drops below unity, i.e., when

[TABLE]

the system remains in the normal state down to $T=0$. This holds when $f$ is non-zero and $l$ is finite. If we fix $f$ and keep increasing $l$, we find that for large enough $l$, i.e., for large enough separation between $E_{c}$ and $\Omega$, $T_{c}$ is still finite. Setting $\lambda_{\rm max}(T_{c})=1$ in (11), we obtain

[TABLE]

The existence of a finite $T_{c}$ at large enough $\ell$ agrees with the general reasoning that at $E_{c}/\Omega\gg 1$, a repulsive Coulomb interaction gets strongly reduced at energies comparable to $\Omega$ and does not overshadows an attraction due to phonon exchange.

II.2 The model with $V_{\rm eff}(\omega)$

We next discuss how the results get modified if we replace separable step-like $V_{\rm RS}(\omega_{n},\omega_{m})$ by $V_{\rm eff}(\omega)$ from Eq. (2), which is a continuous function, and depends on the frequency transfer $\omega=\omega_{m}-\omega_{n})$ rather that separately on $\omega_{m},\omega_{n}$.

The equation for the gap function now reads:

[TABLE]

The continuous-frequency version of (14) is

[TABLE]

It is instructive to consider separately the case $\Omega_{1}=0$, like in Bardeen-Pines model Bardeen and Pines (1955), and the case when $\Omega_{1}$ is finite.

II.2.1 The case $\Omega_{1}=0$

We first solve for $T_{c}$, for which $\lambda_{\rm max}(T_{c})=1$, and then obtain $\lambda_{\rm max}(T)$. To find $T_{c}$, we compare the gap equation at $\omega=0$ and at a finite $\omega$. Setting $\lambda=1$ in (15) yields

[TABLE]

At zero frequency we have

[TABLE]

Because the interaction vanishes at $x=0$, the lower limit of integration over $x$ can be safely set to zero.

At a finite $\omega$, we single out logarithmically divergent term from $\int dx/x$ and write

[TABLE]

One can easily make sure that the last term is free from infra-red singularity, hence the lower limit of the integration over $x$ can be set to zero. This last term scales as $\omega^{2}$ at small $\omega\ll\Omega$ and reduces to $\Delta(0)g\ell$ at $\omega\gg\Omega$. When $\Omega\ll E_{c}$ and $\ell=\ln[E_{c}/\Omega]$ is large, it can be approximated, to logarithmic accuracy, by $g\Delta(0)\ell\omega^{2}/(\omega^{2}+\Omega^{2})$. It then cancels out the equivalent term in the first line in (18), and Eq. (18) reduces to

[TABLE]

Substituting this $\Delta(\omega)$ into the r.h.s. of (17), we obtain the equation on $T_{c}$:

[TABLE]

Hence,

[TABLE]

One can verify that this is exactly the same expression as in Eq. (13) at $f=1$, which in Rietschel-Sham model corresponds to the vanishing of the repulsion at small frequencies.

Returning to (19), we see that the gap function $\Delta(\omega)$ changes sign at $\omega_{c}=\Omega/\sqrt{gL-1}=\Omega\sqrt{g\ell}$ and saturates at high frequencies to $-\Delta(0)/(g\ell)$.

At $\ell=O(1)$, the analysis of $T_{c}$ becomes more involved as the prefactors for $g$ and $g^{2}$ terms in (21) are determined by internal frequencies $x\sim\Omega$ rather than $E_{F}\gg x\gg\Omega$. One way to proceed is to solve the gap equation (16) in direct perturbative expansion in $g$ (see Refs. Karakozov et al. (1975); Dolgov et al. (2005); Wang and Chubukov (2013); Marsiglio (2018)). The computation is lengthy but straightforward. We skip the details and present the result:

[TABLE]

where

[TABLE]

At large $E_{c}/\Omega$, we have $\beta_{1}\approx\beta_{2}\approx 2\beta_{3}\approx\ell$. Keeping terms of order $g\ell$ but neglecting terms of order $g$, we reproduce Eq. (21). At $E_{c}\geq\Omega$, all three constants $\beta_{i}=O(1)$, and Eqs. (22) and (23) yield the result for $T_{c}$ with corrections of order $g$. Note that the leading term in the r.h.s. of (22) is $1/(g^{2}\beta_{2})$, hence the leading exponential dependence of $T_{c}$ is $T_{c}\propto E_{c}e^{-1/(g^{2}\beta_{2})}$. We will see that the presence of $g^{2}$ in the exponent is specific to the case $\Omega_{1}=0$. For finite $\Omega_{1}$, the exponential factor contains $1/g$ rather than $1/g^{2}$.

Eq. (22) gives correct $T_{c}$ in the asymptotic limit $g\ll 1$, but at realistic $g\leq 1$ higher order terms in $g$ may become substantial. To estimate $T_{c}$ at $\Omega\leq E_{c}$ and $g\leq 1$ we take as an input the numerical solution of the actual gap equation (14). It shows (see Fig. 2) that the functional form of $\Delta(\omega)$ is consistent with (19) to reasonably good accuracy, but the prefactor for $\omega^{2}/(\omega^{2}+\Omega^{2})$ is different from $gL$. We then search for the approximate solution using the functional form

[TABLE]

The parameter $\delta$ is determined self-consistently, by substituting (24) into the r.h.s. of (16), expanding it in $\omega$, and matching the $\omega^{2}$ terms. This yields

[TABLE]

Substituting further (24) into (17) we obtain

[TABLE]

Equations (25) and (26) determine $T_{c}$ and $\delta$ self-consistently. The result is

[TABLE]

and

[TABLE]

where

[TABLE]

and

[TABLE]

are numerical coefficients of order unity (approaching, respectively, $3/2$ and $1/12$ in the $E_{c}/\Omega\gg 1$ limit), and

[TABLE]

For small $g$ and large $\ell=\ln[E_{c}/\Omega]$, Eqs. (28) and (22) agree to leading order, but, predictably, differ in corrections of order $g$.

To compare the two formulas, we estimate $T_{c}$ for the same parameters that were used in Fig. 2. Eq. (28) yields $T_{c}=1.8\times 10^{-7}$, and Eq. (22) yields $T_{c}=3.3\times 10^{-7}$. The two values are quite close and also close to the numerical result $T_{c}=2.41\times 10^{-7}$. This situation holds for smaller values of $E_{c}/\Omega$. For the same $g=0.228$ and $E_{c}=\pi$ but $E_{c}/\Omega=e$, Eq. (28) yields $T_{c}=2.8\times 10^{-18}$, and Eq. (22) yields $T_{c}=4.4\times 10^{-18}$, providing accurate estimates of the transition temperature $T_{c}\approx 2\times 10^{-18}$ obtained numerically.

The perturbative and the self-consistent computational schemes can be easily extended to obtain temperature dependence of the largest eigenvalue $\lambda_{\rm max}(T)$. All one needs to do is to substitute $g$ with $g/\lambda(T)$ in Eq. (16) and re-express the condition for $T_{c}$ at $\lambda_{\rm max}=1$ as the equation on $\lambda_{\rm max}(T)$. Within the self-consistent scheme we obtain (for $E_{c}/\Omega\gg 1$)

[TABLE]

The dependence on $L$ is manifestly non-linear.

II.2.2 Finite $\Omega_{1}$

The computations at a non-zero $\Omega_{1}$ in $V_{\rm eff}$ in Eq. (2) proceed in a similar way. We skip the details and present the results. At large $\ell$, keeping powers of $g\ell$ but neglecting powers of $g$, we obtain

[TABLE]

One can easily check that this is equivalent to Eq. (13) for Rietschel-Sham model, if we identify $f$ with $1-\Omega^{2}_{1}/\Omega^{2}$. We see that there is a threshold on finite $T_{c}$ at $\Omega_{1}/\Omega=[g\ell/(1+g\ell)]^{1/2}$. For larger $\Omega_{1}$, frequency variation of the interaction is not sufficient, and $T_{c}=0$. For large enough $g\ell$, $T_{c}$ is finite.

For $\Omega\leq E_{c}$, we again obtain $T_{c}$ in direct perturbative expansion in $g$. The perturbative analysis is only valid at small $\Omega_{1}/\Omega$, otherwise small $g$ are below the threshold. Expanding in $g$ and in $\Omega_{1}/\Omega$, we obtain

[TABLE]

where $\beta_{i}$ are the same as in (23) – corrections to $\beta_{i}$ due to non-zero $\Omega_{1}$ account for the terms, which are smaller than the ones that we kept in (34).

III The numerical analysis of the eigenvalue problem and the breakdown of the logarithmic flow

We now discuss in more detail the numerical solution for $T_{c}$. Like we said in the Introduction, the hallmark of the weak-coupling BCS theory is logarithmic in temperature flow of the largest eigenvalue, $\lambda(T)$, of the linearized gap-function equation (3), see, e.g., Scalapino (1969); Gladstone et al. (1969). We present it graphically in Fig. 1(b) using the notion of the four-point vertex function $\Gamma$ irreducible in the Copper channel. If $\Gamma$ is replaced with a negative constant $-g/\rho_{F}$, where $\rho_{F}$ is the density of states on the Fermi-surface per spin component (i.e., if the interaction is attractive), then the eigenvalue increases with decreasing $T$ as

[TABLE]

Even if $T_{c}$ is extremely small, one can predict its value from higher temperature data by using linear in $\ln[\Omega/\alpha T]$ extrapolation. Such an extrapolation appears to be an indispensable part of any fully ab initio approach in view of the technical challenge of explicitly dealing with the energies/frequencies ranging from the bandwidth down to $T_{c}$.

However, in the case of frequency dependent repulsive interaction, the behavior of the largest eigenvalue $\lambda_{\rm max}(T)$ is not purely logarithmical with $T$, see Eq. (11) for the Rietschel-Sham model and Eq. (32) for the model with $V_{\rm eff}(\omega)$. In both cases, the logarithmic scaling of $\lambda(T)$ fails because it was based on the implicit assumption that the eigenvector does not depend on $T$, while in our case, the $a/b$ ratio in the Rietschel-Sham model and the position of sign change of $\Delta(\omega)$ in the model with $V_{\rm eff}(\omega)$ varies with temperature.

So far, all our considerations were done for momentum-independent $\Gamma$. However, the singular part of the $G^{(2)}$-function is fully symmetric with respect to its frequency dependence and the dependence on the distance to the Fermi surface. An immediate conclusion then is that systems with substantial dependence of the $s-$wave component of $\Gamma$ on the magnitude of the momentum will feature similar properties, if a repulsive interaction is weakened near the Fermi surface. If this is the case, then even a static repulsive interaction may give rise to a finite $T_{c}$ by exactly the same mechanism as the one we discussed above.

III.1 Implicit renormalization approach

At this point one might get an impression that the idea of numerically extracting $T_{c}$ from an eigenvalue/eigenvector problem in terms of genuine $\Gamma$ and $G^{(2)}$ at $T\gg T_{c}$ is hopeless, and the only meaningful way to proceed is the pseudopotential (explicit renormalization) approach (see Refs. Scalapino (1969); Gladstone et al. (1969); Morel and Anderson (1962)). Nevertheless, it turns out that the eigenvalue/eigenvector problem in terms of genuine $\Gamma$ can be reformulated in such a way that the resulting eigenvalue $\bar{\lambda}(T)$ does feature the desired simple logarithmic flow. Moreover, $\bar{\lambda}(T)$ is essentially the eigenvalue of the renormalized problem despite being obtained directly from the genuine $\Gamma$ without explicitly constructing the pseudopotential counterpart for the latter.

III.1.1 The Rietschel-Sham model

The Rietschel-Sham model (1) is perfectly suited for introducing our method, especially for tracing the intrinsic connection to—and the fundamental structural difference from—the pseudopotential approach. First, let us change the parametrization of the gap function:

[TABLE]

This way we separate $\Delta$ into two distinctively different parts: the high-frequency part, $|\omega_{n}|>\Omega$, (fully described by the parameter $a$) and the the low-frequency part, $|\omega_{n}|<\Omega$, (fully described by the parameter $s$). With the new parameterization—but precisely the same eigenvector and eigenvalue, the problem (10) is reformulated as

[TABLE]

Make the formal replacement $\lambda\to 1$ in the second equation:

[TABLE]

Now make a straightforward observation that the requirement $\bar{\lambda}(T_{c})=1$ leads to precisely the same $T_{c}$—with precisely the same eigenvector $[a(T_{c}),s(T_{c})]$—as the original system (39). The utility of replacing (39) with (42) becomes immediately clear by eliminating the high-frequency part $a$ between the two equations:

[TABLE]

The flow of $\bar{\lambda}$ is obviously linear in $L$. While the value of $T_{c}$ is, of course, independent of the method, see Fig. 3, the crucial difference is in the simplicity of extrapolating data towards exponentially low temperature—in more complex models one may not be able to solve for eigenvalues below $T\gg T_{c}$. For example, if we scale the coupling constant $g\to\epsilon g$ to smaller values to ensure that the condition (12) is satisfied, the flow of $\lambda_{\epsilon}=\epsilon\lambda$ remains qualitatively the same at $T\gg T_{c}$, see Fig. 4. Reliable prediction of $T_{c}$ under these conditions would be nearly impossible; moreover, one would be left wondering whether the model features an effective attraction at low energies and ultimately goes SC. In contrast, under the same conditions, negative $\bar{\lambda}_{\epsilon}$ would immediately signal that the model is not SC in the corresponding channel.

By the very fact of eliminating the high-frequency part we understand that Eq. (43) corresponds to the peseudopotential theory with the temperature flow of $\bar{\lambda}(T)$ controlled by the effective coupling constant $g_{\rm eff}$:

[TABLE]

We indeed see that the structure of $g_{\rm eff}$ reproduces the Tolmachov-McMillan logarithm Tolmachov and Tiablikov (1958); Morel and Anderson (1962); Rietschel and Sham (1983); i.e. it demonstrates that a repulsive interaction is renormalized to a smaller value $g/(1+g\ell)$ at low frequencies $(\omega_{n},\omega_{m})\lesssim\Omega$, and SC instability is possible if this renormalized value is smaller than the bare low-frequency attractive term $fg$. Diagrammatically, this result follows from the ladder summation for the effective low-frequency $\Gamma$-function shown in Fig. 1(c): $\Gamma_{\rm eff}=(g-gf)-g^{2}\ell+g^{3}\ell^{2}-\dots$

While for the utterly simple Rietschel–Sham model (1) the distinction between the formulation (42) and the explicit pseudopotential formulation (43) is merely nominal, the difference becomes profound for any realistic case. Here the equivalents of the numbers $a$ and $s$ are the high- and low-frequency parts of (momentum- and frequency-dependent) $\Delta$. Correspondingly, the system (42) becomes the system of coupled integral equations (with the kernels given by genuine $\Gamma$). At not too low temperature, the resulting eigenvector/eigenvalue problem remains solvable by techniques of DiagMC—even without explicitly evaluating $\Gamma$ (see Sec. IV). In contrast, the elimination of the high-frequency part required for going from (42) to (43) would face the challenge of numerically constructing the vertex function $\Gamma_{\rm eff}$ out of the multi-variable and multi-scale vertex function $\Gamma$ (putting aside the non-trivial problem of obtaining $\Gamma$ from first principles in a correlated system).

So far, we discussed the protocol of extracting $T_{c}$, but not the critical gap function. It turns out, however, that the desired solution is immediately related to the eigenvector/eigenvalue problem of the implicit-renormalization scheme, by simply multiplying the low-frequency part of the eigenvector, obtained at a given temperature $T$, by the factor $\bar{\lambda}(T)$. This relationship is readily traced with the model (1) when the low-frequency part of the solution is described by a single parameter $s$. We introduce the re-weighted quantity

[TABLE]

and re-express (42) in terms of $s_{*}$, thus arriving at the eigenvalue/eigenvector problem for $(s_{*},a)$:

[TABLE]

where, by construction, $g_{\rm eff}=\bar{\lambda}(T)/L$. The problem (48) is free of any temperature dependence, meaning that the vector $(s_{*},a)$ is temperature-independent. We now recall that at $T=T_{c}$, ${\bar{\lambda}}=1$, hence $s_{*}=s$, see Eq. (45). This implies that $(s_{*},a)$ coincides with the critical gap function.

III.1.2 The model with $V_{\rm eff}(\omega)$

The same computational scheme can be use to obtain $T_{c}$ in the model with $V_{\rm eff}(\omega)$, Eq. (2). The formulation of the implicit renormalization approach is as described above, but its practical implementation is different for the final step. Here we describe it using general vector-matrix notations (closely following Ref. Rietschel and Sham (1983)), when the original $T_{c}$ problem can be written in a compact form as

[TABLE]

First, the gap function is decomposed into two complementary parts (low- and high-frequency projections), $\vec{\Delta}\equiv\vec{\Delta}^{(1)}+\vec{\Delta}^{(2)}$, such that $\Delta^{(1)}_{n}=0$ for $|\omega_{n}|>\Omega_{c}$ and $\Delta^{(2)}_{n}=0$ for $|\omega_{n}|<\Omega_{c}$. Correspondingly, the $\hat{A}$-matrix is decomposed into four complementary matrixes, $\hat{A}=\hat{A}^{(11)}+\hat{A}^{(22)}+\hat{A}^{(21)}+\hat{A}^{(12)}$, such that $\hat{A}^{(11)}$ and $\hat{A}^{(22)}$ have zero matrix elements between low- and high-frequency subspaces, while the only non-zero matrix elements of $\hat{A}^{(21)}$ and $\hat{A}^{(12)}$ are those connecting low-to-high and high-to-low frequency subspaces, respectively. In analogy with (42), we then consider the eigenvector-eigenvalue problem

[TABLE]

Along the same lines, the relationship between this problem and the explicit renormalization approach is readily established by formally substituting $\vec{\Delta}^{(2)}=-[\hat{I}+\hat{A}^{(22)}]^{-1}\hat{A}^{(21)}\vec{\Delta}^{(1)}$ (the equality implied by the second equation) into the first equation. This yields, in analogy with Eq. (44),

[TABLE]

where $\hat{B}=\hat{A}^{(11)}-\hat{A}^{(12)}[\hat{I}+\hat{A}^{(22)}]^{-1}\hat{A}^{(21)}$ is the renormalized kernel in the Cooper channel. Its diagrammatic expansion in terms of the bare $\Gamma$ has the same structure as that of the series shown in Fig. 1(c).

When the problem (52) is solved—the procedure is described in the next section—for model (15) with $\Omega_{c}\lesssim\Omega$, the result is a nearly perfect linear dependence of $\bar{\lambda}$ on $\ln\Omega_{c}/T$ at low temperature, see Fig. 5.

In Fig. 5, we also compare the numerical solutions for ${\bar{\lambda}}$ for the case $\Omega_{1}=0$ to the analytical solution, Eq. (32), and the numerical solution of the the original eigenvalue problem, Eq. (15), using the same model parameters as in Fig. 2. We see that ${\bar{\lambda}}(T)$ is linear in $\ln T$ all the way down to $T_{c}$, and this allows one to determine $T_{c}$ in a controllable way by extrapolating ${\bar{\lambda}}(T)$ from higher temperatures. We illustrate this in Fig. 6, where we show the results of the extrapolation using various high-frequency cutoffs $\Omega_{c}$. We see that the extrapolation of ${\bar{\lambda}}(T)$ yields the correct $T_{c}$ if we set $\Omega_{c}<\Omega$. On the other hand, we clearly see from the Figure that the non-linear dependence of the original $\lambda_{\rm max}(T)$ on $\ln T$ results in a large overestimate of $T_{c}$, if the extrapolation is done by using the original $\lambda_{\rm max}(T)$, or by setting the cutoff at $\Omega_{c}\gg\Omega$. It is also clear that if we were to choose $\Omega_{c}$ much smaller than $\Omega$, it would be close to the extrapolation interval, and this would result in increased systematic error.

Finally, we want to make sure that the vector $[\vec{\Delta}_{*}^{(1)},\,\vec{\Delta}^{(2)}]$, wehre [cf. (45)]

[TABLE]

yields the critical gap function. Rewriting the problem (52) in terms of $[\vec{\Delta}_{*}^{(1)},\,\vec{\Delta}^{(2)}]$, we get

[TABLE]

The structure of this problem is similar to that of (48), with an important distinction. While in Eq. (48) the temperature dependence is absent, Eq. (57) becomes effectively temperature independent only within the leading logarithmic accuracy, when the logarithmic factor brought by the integration with the kernels $\hat{A}^{(11)}$ and $\hat{A}^{(21)}$ is compensated by matching behavior of $\bar{\lambda}(T)$. This means that $[\vec{\Delta}_{*}^{(1)},\,\vec{\Delta}^{(2)}]$ reproduces the critical gap function only within the leading logarithmic approximation. To recover an accurate result for the critical gap function one has to extrapolate $[\vec{\Delta}_{*}^{(1)}(T),\,\vec{\Delta}^{(2)}(T)]$ to $T=T_{c}$.

IV Implicit renormalization approach: techniques of implementation

IV.1 Iteration scheme

Fully ab initio calculation of $T_{c}$ in metals using Eq. (52) faces two technical problems. To begin with, the vertex function has two four-dimensional (momentum-frequency) indexes and its full tabulation, including high-energy scales $\gg E_{F}$, is challenging; without simplifying assumptions the effort is about the square of what is required for tabulating the Green’s function $G_{{\mathbf{p}},m}$. [Here we work with the extended momentum space; otherwise ${\mathbf{p}}$ has to be understood as a composite label based on the momentum in the first Brillouin zone and band index.] Suppose that this problem is taken care of.

Next, we need to solve the equation for $\vec{\Delta}^{(2)}$ [the second equation of the system (52)] on a fine four-dimensional grid covering all relevant energy scales from $\gg E_{F}$ to $\ll\omega_{D}$. At this point, inversion of the $\hat{I}+\hat{A}^{(22)}$ matrix will pose a serious, if not unsolvable, problem because of the huge matrix size (which can be reduced by employing coarse-graining description only at the expense of accuracy and additional technical complexity). An attempt to solve the equation by standard iterations,

[TABLE]

will certainly fail because for Coulomb systems, the largest positive eigenvalues of $\hat{A}^{(22)}$ are very large as noted by Rietschel and Sham Rietschel and Sham (1983). The negative eigenvalues of $\hat{A}^{(22)}$ are all smaller than unity by the very statement of the problem—otherwise the system will go SC at temperature above $\Omega_{c}$.

It turns out that the desired solution can be always found by a simple modification of the iteration scheme. The idea follows from known convergence properties of “damped” iterations, namely, $f=-af-b$ with $a>-1$ can be solved by substituting into the r.h.s. the average of all previous iterations:

[TABLE]

(see, for example, Ref. Prokof’ev and Svistunov (2007)). One can easily verify an extremely fast convergence of this scheme even for large values of $a$. It is thus guaranteed that

[TABLE]

will quickly converge to the desired solution, which in its turn can be subsequently used to find the largest eigenvalue of (52) by the standard power method:

[TABLE]

Here $\hat{B}\vec{\Delta}^{(1)}$ is a short-hand notation for the r.h.s. of the first equation of the set (52) once we use for $\vec{\Delta}^{(2)}$ the solution of the second equation in (52).

In practice, we have two (internal and external) cycles and work with the original matrix $\hat{A}$. In the internal (damped iterations) cycle, $\hat{A}\vec{\Delta}$ multiplication involves the full vector $\vec{\Delta}$ but the result is used to compute and average only the high-frequency components; the low-frequency part remains intact. After convergence, the same $\hat{A}\vec{\Delta}$ multiplication is used in the external cycle to obtain the new vector, estimate the eigenvalue $\bar{\lambda}$, and normalize the low-frequency components. The norm $||\vec{\Delta}^{(1)}||$ can be defined in a number of ways, say, through the inner product.

To improve efficiency, all equations should be projected on different symmetry channels, and their largest eigenvalues be determined independently. This is important for correctly predicting the outcome of close competition between two or more channels, because the order of the largest eigenvalues can potentially change with temperature.

IV.2 Diagrammatic Monte Carlo method

Within the DiagMC approach, where statistics is accumulated and averaged as a matter of principle, damped iterations are realized naturally Prokof’ev and Svistunov (2007). The internal cycle is automatically realized by running a self-consistent scheme when statistics for high-frequency components of $\Delta$ is accumulated by sampling the diagrammatic space for the r.h.s., $\Gamma G^{(2)}\Delta$, where $\Delta$-function is known from averaging previously collected statistics.

The DiagMC approach also circumvents the difficulty of tabulating the vertex function explicitly. Formally, the $\Gamma$-function is represented by the series of Feynman diagrams and thus its evaluation involves summation over diagram orders and topologies, as well as multi-dimensional integrals and sums over internal momenta and frequencies. In DiagMC, all diagrammatic space parameters, including external variables, are sampled stochastically without systematic bias. In this sense, the r.h.s. of Eqs. (52) are subject to the same DiagMC simulation as electron self-energy or polarization function. At no point one has to worry about handling an object more complex than the single-particle Green’s function. In fact, the gap function is more simple than $G$ because it lacks singular momentum and frequency dependance, see Fig. 2.

V Discussion and Conclusions

In this communication we considered superconductivity in systems where the interaction is repulsive, but depends on frequency, and is weaker at smaller frequencies than at larger ones. This is a typical systems in a metal, where the interaction in the particle-particle channel is a combination of a stronger repulsion due to Coulomb interaction, and a weaker attraction due to electron-phonon interaction.

There are at least three characteristic energy scales in metals: the Fermi energy $E_{F}$, the plasmon frequency $\omega_{p}$, and the Debye frequency $\Omega_{D}$. Screening of the Coulomb interaction develops at frequencies below $\omega_{p}$, but for small momenta $q\ll k_{F}$ it is not complete down to $\omega\leq v_{F}q$, where $v_{F}$ is the Fermi velocity. The electron-phonon interaction becomes important at $\omega\lesssim\Omega_{D}$. In the conventional pseudopotential approach, all effects of repulsive Coulomb interaction, including its momentum and frequency dependence, are absorbed into a single semi-phenomenological dimensionless coupling $\mu_{*}$ (see Ref.  Giustino (2017) for a recent review on this issue). Superconductivity develops if a dimensionless $\lambda_{ph}$ due to electron-phonon interaction exceeds $\mu^{*}$.

In our analysis of superconductivity somewhat different approach, inspired by recent studies of $s$-wave superconductivity in SrTiO3 and Bi. Namely, we assumed, following earlier studies Mahan (2000); Ruhman and Lee (2016, 2017); Gastiasoro et al. (2019), that the effective interaction in the particle-particle channel can be viewed as the sum of a screened Coulomb interaction and an interaction with a gapped boson, which is a hybridized mode between a longitudinal phonon and a plasmon. We considered s-wave superconductivity and focused on the frequency dependence of the effective pairing interaction, and neglected its momentum dependence.

We considered the two models: the Rietschel-Sham model with a step-like pairing interaction, and the model with a continuous pairing interaction $V_{\rm eff}(\omega)$ (Eq. (2)). In both models the pairing interaction reduces to a larger constant $g$ at high energies and to a smaller constant $g(1-f)$ at the lowest energies. We found that $T_{c}$ is finite already at arbitrary small $g$, if $f=1$, i.e., if the interaction vanishes at the smallest frequencies. However, if $f$ is smaller than one, $T_{c}$ is finite only if $g$ exceeds a certain threshold. The threshold values and the results for $T_{c}$ for the two models are similar, but not equivalent. For the model with a continuous interaction $V_{\rm eff}(\omega)$ we computed $T_{c}$ using various computational schemes, and analyzed the interplay between the threshold value and the ratio $\Omega/E_{c}$, where $\Omega$ is the boson frequency and $E_{c}\leq E_{F}$ is the frequency, up to which one can expand the fermionic dispersion to linear order in $k-k_{F}$. We assumed that $\Omega<E_{c}$ and showed that the threshold value is reduced when $\Omega/E_{c}$ gets smaller.

We also discussed in all detail the protocol for numerical computation of $T_{c}$ for systems with frequency dependent repulsive interaction, particularly for the cases when $T_{c}$ is small, and one needs to extrapolate the results for the largest eigenvalue of the gap equation, $\lambda_{\rm max}(T)$, from temperatures $T\gg T_{c}$ to $T=T_{c}$. We demonstrated that within the standard setup, used for the systems with frequency independent attractive interaction, such an extrapolation is not possible due to non-linear dependence of $\lambda_{\rm max}(T)$ on the logarithm of temperature. This non-linear dependence emerges because frequency-dependent repulsive interaction undergoes significant changes at intermediate energy/momentum scales.

To set the new protocol, we noticed that the numerical calculations necessary involve the energy cutoff at $\Omega_{c}$. In most calculations, this cutoff is set at $\Omega\ll\Omega_{c}\ll E_{F}$. This is made for purely technical reasons—to have better momentum resolution near the Fermi surface. However, with this choice of $\Omega_{c}$ one faces the problem of extrapolating non-linear in $\ln(\omega_{D}/T)$ data towards low temperature, if $T_{c}$ happens to be much smaller than the lowest possible temperature in the calculation. Another, more general drawback, is that for $\Omega_{c}\ll E_{F}$, the effects of Coulomb interactions are no longer included at the fully ab initio level and the subsequent calculation contains an unknown systematic error, not to mention that momenta satisfying $q<\Omega_{c}/v_{F}$ are treated inadequately. Moreover, in multi-band or strongly anisotropic systems the effects of Coulomb interaction cannot be described by a single parameter.

We argued that accurate evaluation of $T_{c}$ from Fermi-liquid properties at $T\gg T_{c}$ in correlated systems, where the BCS regime is an emergent phenomenon involving multiple energy scales, cannot be achieved unless the renormalization scale $\Omega_{c}$ is made smaller (much smaller) than $\Omega$. This choice, however, brings about two technical problems when it comes to the practical implementation of the method following the protocol described in Ref. Rietschel and Sham (1983). One is the necessity to know the full vertex function in a broad frequency and momentum range; the amount of information is about the square of that required for knowing the single-particle Green’s function. Even if $\Gamma^{{\mathbf{k}},n}_{{\mathbf{p}},m}$ can be tabulated without approximations and any loss of accuracy, solving for high-frequency components of the gap function by matrix inversion would be impossible because of the huge matrix size. Finding the solution by standard iterations will not work either because negative eigenvalues for the full problem are largest in modulus and for Coulomb systems will exceed unity already at $\omega\ll T\ll E_{F}$ Rietschel and Sham (1983).

A protocol for extrapolating numerical data towards $T_{c}$ from higher temperatures—applicable to first-principle description of real metals—has to adequately capture the physics of the emergent weakly-interacting effective theory. We have formulated the so-called implicit renormalization approach and demonstrated that it provides a simple, efficient, and unbiased protocol for solving the extrapolation problem. The scheme has a built-in tool of controlling the systematic error of extrapolation (see Fig. 6)—the only systematics of the otherwise numerically exact method. The implicit renormalization approach is perfectly compatible with the diagrammatic Monte Carlo techniques, allowing one to solve the corresponding eigenvalue problem without invoking the matrix inversion or even explicitly calculating the four-point vertex function $\Gamma$. The implicit renormalization protocol also allows one to obtain the correct gap function immediately below $T_{c}$.

Throughout the paper, the separation of the gap function into the low-energy part $\vec{\Delta}^{(1)}$ and the higher-energy part $\vec{\Delta}^{(2)}$ was performed in the frequency domain. Our approach, however, can be readily extended to include this separation for both the frequency and momentum variables. For example, the condition on the “low-energy” regime, where $\vec{\Delta}({\bf p},\omega)=\vec{\Delta}^{(1)}({\bf p},\omega)$ can be as simple as $\xi_{p}^{2}+\omega^{2}\leq\Omega_{c}^{2}$. Outside this range, $\vec{\Delta}({\bf p},\omega)=\vec{\Delta}^{(2)}({\bf p},\omega)$.

Acknowledgements.

We thank M. Gastiasoro, R. Fernandes, D. Maslov, and A. Millis for fruitful discussions. The work by AVC was supported by the Office of Basic Energy Sciences U. S. Department of Energy under award DE-SC0014402.