Precursory Cooper Flow in Ultralow-Temperature Superconductors

Superconductivity at low temperature -- observed in lithium and bismuth, as well as in various low-density superconductors -- calls for developing reliable theoretical and experimental tools for predicting ultralow critical temperatures, $T_c$, of Cooper instability in a system demonstrating nothing but normal Fermi liquid behavior in a broad range of temperatures below the Fermi energy, $T_{\rm F}$. Equally important are controlled predictions of stability in a given Cooper channel. We identify such a protocol within the paradigm of precursory Cooper flow -- a universal ansatz describing logarithmically slow temperature evolution of the linear response of the normal state to the pair-creating perturbation. Applying this framework to the two-dimensional uniform electron gas, we reveal a series of exotic superconducting states, pushing controlled theoretical predictions of $T_c$ to the unprecedentedly low scale of $10^{-100} T_{\rm F}$.

Introduction.-Theconceptual elegance of the Kohn-Luttinger theorem establishing that the Fermi liquid state is unstable in high-angular momentum Cooper channels at low temperatures comes at the price of lacking accurate predictions as to which channels get unstable first and at what temperatures.Experimentally, each new discovery of the ultralow-temperature superconductor, be it lithium [1] or bismuth [2], exhibiting superconductivity at about 0.1 mK, emerges as a surprise.The potential for observing analogous phenomena in traditionally non-superconducting metals such as gold, copper, or sodium, as well as in low-density superconductors [3][4][5], adds to the scientific intrigue, with no a priori knowledge of what the critical temperature T c one should expect for a given system.
There is, however, a fundamental reason to expect that the desired answers can be controllably extracteddefinitely theoretically and, hopefully, experimentally as well-from the system's properties in the normal Fermi liquid regime at temperatures much lower than the Fermi energy, T F , but still many orders of magnitude higher than T c .Indeed, the (ultra-)low value of T c is due to the emergent Bardeen-Cooper-Schrieffer (BCS) regime of instability, when the strongly correlated at the ultraviolet level system gets renormalized into the Fermi liquid with a weak effective BCS interaction characterized by a dimensionless negative coupling constant |g| ≪ 1.The result is an exponentially small critical temperature T c ∼ T F e 1/g .The BCS nature of the transitions implies a rather characteristic temperature evolution of the pair susceptibility, χ 0 (T ), defined as the linear response to a static pair-creating perturbation.On the approach to the critical point, χ 0 (T ) should behave as [6,7] χ 0 (T ) ∝ 1/ ln (T /T c ) (T → T c + 0) .
While providing a proof-of-principle result for the idea of extracting T c from properties of the normal state at T ≫ T c , relation (1) turns out to be rather impractical, and sometimes even misleading, when it comes to a controlled quantitative analysis of (in)stability in a given pairing channel (for an illustration see the blue curve in Fig. 1).The reason is that ansatz (1) ignores a slowly evolving with temperature logarithmic prefactor (its physical origin is discussed below), making a naïve extrapolation of an apparent linear dependence of 1/χ 0 (T ) on ln T from high temperature to the temperature when it is supposed to hit zero very inaccurate, not to mention that it may predict finite T c for a Cooperstable channel (see a similar discussion in Ref. [14]).
Recently, a numeric method-the so-called implicit renormalization (IR)-allowing one to accurately predict T c from the field-theoretical properties of the system at T ≫ T c has been proposed in Ref. [14] and further developed in Refs.[15,16], with an application to the model of uniform electron gas.Despite unquestionable success, the IR approach encounters certain technical limitations and lacks a direct connection with what can be measured experimentally.Technical limitations of IR are most pronounced in the vicinity of the "quantum transition point" (QTP) at which a given channel undergoes a transition from Cooper-stable to Cooper-unstable regime.It is thus crucial to find a complementary to the IR approach and one compatible with experimental protocols.
In this Letter, we show that the desired solution is provided by nothing but a more accurate than (2) but still physically transparent expression for the pair susceptibility.Specifically, we find that the following (universal to all ultralow-temperature BCS-type superconductors) three-parametric ansatz-the so-called precursory Cooper flow, or PCF ansatz-perfectly captures the temperature evolution of χ 0 within a broad temperature range: The non-universal parameters c, g, and Λ are tied to the microscopic properties of the system, with Λ being the lowest relevant energy/frequency scale (we set ℏ = k B = 1), such as T F , Debye, or plasma frequency.Negative g implies the BCS transition at T c = Λe 1/g , while g > 0 implies its absence, with g = 0 corresponding to QTP.The logarithmic factor in the numeratordistinguishing (2) from (1)-has the same mathematical origin and, thus, the same expression as the "Tolmachev's logarithm" in the denominator.However, the two logarithms describe distinctively different physics.The one in the numerator is the pair susceptibility of an ideal Fermi liquid (a system with no coupling in the Cooper chan-nel), while the one in the denominator is responsible for Tolmachev's renormalization of the effective interaction [17,18].Sharp difference between the stable and unstable regimes develops only at |g| ln(Λ/T ) ≫ 1; otherwise, susceptibility increases in both regimes regardless of interactions.In the former case, χ 0 saturates to c/g at temperatures T ≪ T * ≃ Λe −1/|g| .These pair-correlations (diverging as g → 0) play a crucial role in the scenario of strange metal behavior discussed in Ref. [19].For a model with weak momentum-independent interaction, the expression (2) is readily obtained by Bethe-Salpeter summation of the Cooper-channel diagrammatic ladder.Far less trivial is our result shedding new light on the previous IR observations [ [14][15][16] that Eq. ( 2) also works in the case of dynamically screened Coulomb interaction with complex momentum and frequency dependence of the effective coupling in the Cooper channel.
In the context of ab initio calculations employing the PCF methodology, we introduce an optimized fieldtheoretical counterpart of the pair susceptibility, R 0 .As opposed to χ 0 , the flow of R 0 is free of the "confusing" ideal-Fermi-liquid logarithmic numerator and is characterized by only two parameters: [Consistency with (2) implies ln(Λ ′ /Λ) = 1/g ′ − 1/g.]This yields an exciting opportunity for precise theoretical and numerical determination of T c and QTP from normal state calculations using a minimal number of fitting parameters (see Figs. 1 and 2).With the precise method at hand, we performed a model study of the two-dimensional (2D) uniform electron gas (UEG) in the regime of weak-to-moderate interactions, which is interesting for its intrinsic (nophonons) superconductivity driven by the dynamically screened Coulomb interaction [20][21][22][23][24] as opposed to the original Kohn-Luttinger scenario [25].Our results (see Fig. 3) reveal a series of QTPs associated with ultralowtemperature superconducting instabilities that we can resolve down to 10 −100 T F .
Field-theoretical analysis.-Withoutloss of generality, we study the universal s-wave linear response scaling laws in a d-dimensional spherically symmetric homogeneous system (other channels and realistic superconductors are discussed in the Supplemental Material (SM) [26]).The linear response functions (2) and (3) originate from the two-electron Green's function with zero incoming momentum and frequency, G , where G k denotes the dressed one-electron Green's function.
We start by analyzing the analytic structure of R k as it follows from the self-consistent Bethe-Salpeter equation: where Γ is the particle-particle irreducible four-point vertex with zero incoming momentum and frequency.It encodes all effective pairing interactions, such as screened Coulomb potential.The second term in the r.h.s. of ( 4) is a sum of ladder diagrams generated by repeated products of Γ and GG, each carrying its own set of singularities.The finite-temperature cutoff of GG at the Fermi surface is responsible for the logarithmic flow that ultimately leads to BCS instability.Concurrently, the vertex function Γ has singular momentum dependence due to incomplete screening of long-range Coulomb interaction at any finite frequency.Possible interplay between the two singularities raises the question of whether the flow of the pair susceptibility still follows the same law as in the case of short-range interaction.
The key observation is that Coulomb singularity does not produce large terms when ΓGG is integrated over p, and the dominant contribution still comes from the BCS logarithm.That is, p Γ kp G p G −p = gk ln T + fk + O(T ), where gk and fk are temperature independent and regular in k functions, and the finite-T corrections vanish at least linearly with T .Further technical details are provided elsewhere [27].This observation allows one to parameterize ΓGG as where p δ p = 1), and the regular correction satisfies p ϕ kp = O(T ).By incorporating this form into Eq.( 4), we obtain the temperature dependence of the linear response (see SM [26]): where f k and g k are regular functions representing fk and gk renormalized by pair fluctuations.Remarkably, the logarithmic correction in the numerator vanishes in the low-energy limit k → 0, resulting in a simple relation for R 0 ≡ R k→0 given by Eq. ( 3) with g ′ = g 0 and Λ ′ = e −f0/g0 .The pair susceptibility, on the other hand, involves R k with finite k and, thus, retains the global logarithmic factor.Superconductivity in the 2D uniform electron gas.-In the absence of electron-phonon interaction, superconduc-tivity in this model is of the emergent BCS type, and the values of T c are supposed to be extremely low.The theory presented above and advanced numerical techniques not only allow us to study many paring channels but also accurately locate QTP when T c goes to zero.
Previous studies [15,22] revealed the existence of high orbital momentum ℓ superconducting states in the 3D UEG in the weak coupling limit when the random phase approximation (RPA) becomes controllably accurate.The pairing comes from dynamic nature screening in Coulomb systems, as discussed by Rieschel and Sham [20,28], and not from the static Kohn-Luttinger mechanism [25].In light of the pioneering work by Takada [21], we expect that the dynamic character of screening will play a crucial role in 2D as well.
In the weak coupling limit, the RPA vertex function Γ has the form where V q = 2πe 2 /q is the bare Coulomb repulsion and Π 0 (p − k, ω m − ω n ) is the RPA polarization.At this level of accuracy, we consider bare Green's functions 4), where m is the electron mass [Fermi temperature can be expressed as To efficiently solve Eq. ( 4) at ultralow temperatures, we employ the discrete Lehmann representation (DLR) [16,29,30] to radically reduce the computational cost of representing linear response and Green's functions from O(1/T ϵ) for uniform grids to O[ln(1/T ) ln(1/ϵ)], where ϵ is the error tolerance.The codes are available online [31].
We performed systematic calculations of R 0 for various orbital channels ℓ and values of r s .For each ℓ and r s , we determine the transition temperature by extrapolating 1/R 0 (T ) flows to zero using the least-squares criterion.Subsequently, for each ℓ, we extrapolate results for −1/ ln[T c (r s )/T F ] to zero to reveal critical values of r s .The phase diagram of competing superconducting states in the 2D UEG is shown in Fig. 3.
As expected, by increasing density (decreasing r s ), one suppresses T c in all orbital channels.While all channels are superconducting at r s ∼ 1, they successively go through QTPs so that only large ℓ channels stay superconducting at small r s .Accordingly, the orbital index of the dominant (highest T c ) channel also increases with density (we obtain it from crossing points between the T Data in the inset in Fig. 3 suggest that superconductivity in the 2D UEG survives in the high-density limit (r s → 0) for large enough ℓ.We emphasize that this outcome cannot be explained by the Kohn-Luttinger mechanism, because in 2D this mechanism simply does not exist at the RPA level [34].Similar to the 3D case, we are dealing with the consequences of dynamic screening in systems with long-range interactions.To explore how interaction range changes the picture, we took density corresponding to r s = 0.8 (when all channels are superconducting) and replaced Coulomb potential with the Yukawa one.We find that for screening length ∼ 1/k F , all channels mentioned in Fig. 3 are no longer superconducting.For more details, see SM [26].
Finally, we recalculated the phase diagram by accounting for renormalization of the Green's function within the G 0 W 0 approximation for the proper self-energy.Qualitative (and quantitative at small r s ) agreement between phase diagrams presented in the SM [26] and Fig. 3 demonstrates the robustness of our conclusions.
Conclusions and outlook.-Wehave shown that pair susceptibility (linear response to a static spatially uniform pair-creating perturbation) of the normal Fermi liquid features universal for all BCS superconductors temperature dependence, or "flow," Eq. ( 2), irrespective of the emergent pairing mechanism.The ansatz (2) applies to both stable and unstable pairing channels.In both cases, the higher-temperature part of the flow is the same, up to small corrections, and represents singular in the T → 0 limit response of an ideal Fermi liquid.A sharp difference between the stable and unstable cases develops only at exponentially low temperatures: the unstable channel hits finite-temperature singularity at T c while the stable channel develops non-trivial correlations suppressing the zero-temperature singularity.The T = 0 singularity survives only at the quantum transition points (QTPs) separating the stable and unstable regimes.
Using two-dimensional uniform electron gas as a paradigmatic model of intrinsic superconductivity mediated by dynamic screening of Coulomb interaction, we demonstrated that fitting the normal Fermi liquid data to the ansatz (2)-and its numeric counterpart (3)-allows one to accurately predict ultralow critical temperatures of unstable Cooper channels and locate QTPs.
We anticipate that our method for controlled quantitative predictions of ultralow-temperature Cooper instability (or its absence) from finite-temperature flows of the linear response to a spatially uniform pair-creating perturbation can be extended to cases where the perturbation is applied at the boundaries of the normal state.This would be the theoretical basis for experimental studies of precursory Cooper flows in the normal metallic state by using superconductor-normal metalsuperconductor tunnel junction setups [8][9][10][11][12][13] in the hightemperature range T c ≪ T ≪ T F .For metals such as Cu, Au, and Na, where superconductivity at ultra-low temperature remains uncertain, it would be equally informative to observe whether flows (2) correspond to negative or positive values of g in the s-channel.

ANGULAR MOMENTUM DECOMPOSITION IN TWO AND HIGHER DIMENSIONS
In an isotropic system, the particle-particle four-point vertex with zero incoming momentum and frequency Γ k,ωn p,ωm only depends on the angle between k and p (θ kp ≡ θ k −θ p), their moduli, and the Matsubara-frequencies difference ω n − ω m .This allows us to simplify the analysis by projecting Γ onto different orbital channels, ℓ.In this section, the frequency variables and Mataubara summation are omitted for brevity as they do not affect the decomposition.
Considering the Bethe-Salpeter equation for the linear response function R, where η k is the sourced term and F p = G p G −p R p , we can project it onto decoupled orbital channels ℓ via an expansion.After the projection to the angular momentum sector, R and Γ are solely dependent on the momen-tum amplitudes, allowing a unified treatment of generic dimension.
In two dimensions, the vertex function can be expressed as a Fourier expansion (2) Apart from the Γ expansion, Eq. ( 2), we also need to expand functions that depend only on one momentum variable: where O represents R, F , or η, O  1) can also be expanded in Fourier series: where p and R (ℓ) p = 2π 0 dθ p cos ℓθ pR p (the product of two Green's functions is isotropic).Therefore, the linear response function for each channel is given by arXiv:2303.03624v2 [cond-mat.supr-con]14 Sep 2023 Here by setting η (ℓ) ≡ 1 and projecting into channel ℓ, we break the U (1) symmetry and the rotation symmetry so that the result for different channels can be obtained separately.
In three or higher dimensions, the angular momentum decomposition can be achieved by a similar approach using expansion in the spherical harmonics.The vertex function can be expressed as an expansion in Legendre polynomials P ℓ (x) where x = cos θ kp and denotes the number of linearly independent Legendre polynomials of degree ℓ in d dimensions.Using the addition theorem for spherical harmonics [1] where ) is the solid angle in d dimensions and Y ℓm (θ) is the spherical harmonic function, vertex function can be expressed further on as Here we set m = 0 because the decomposed equation is independent of m for a system with rotation symmetry.The linear response function can also be decomposed as where Projecting the convolution term in Eq. ( 1) on the spherical harmonics, we obtain ) where the orthonormal property of the spherical harmonics as dθ kY ℓm (θ k)Y * ℓ ′ m ′ (θ k) = δ ℓℓ ′ δ mm ′ is utilized.Finally, we have the decoupled self-consistent linear response equation for each orbital channel in d(≥ 3) dimensions as where only the measure of momentum integral is different from two dimensions.We can unify both Eq. ( 5) and Eq. ( 11) with the following simplified expression where we introduce the Capital-letter momentumfrequency variable P = (∆p, ω m ) with ∆p = p − k f and G (2) P ≡ G P G −P .The integral is defined as dP ≡ T m d∆p.We have also defined Γ K,P = p 4π 2 Γ (ℓ) k,p in two dimensions and k,p in higher dimensions to absorb the measure of momentum integral.The ℓ label in R (ℓ) k is omitted.

SUPERCONDUCTING PHASES OF THE TWO-DIMENSIONAL ELECTRON GAS WITH G0W0 SELF-ENERGY RENORMALIZATION
We employed the precursory Cooper flow approach to investigate dynamic-screening-driven superconductivity in the two-dimensional uniform electron gas (2D UEG) model, parametrized by the Wigner-Seitz radius where k F is the Fermi momentum and a B is the Bohr radius.As shown in Fig. 1 of the main text, we obtained a rich phase diagram for several channels.
To verify the robustness of the pairing mechanism in the 2D UEG, here, we computed the same phase diagram using the RPA interaction for Γ and the G 0 W 0 renormalized propagator G r (referred to as RPA-G R ), as shown in Fig. 1.Comparison with the RPA-G 0 results presented in the main text reveals small shifts of the phase boundaries, but all key features of the RPA-G 0 phase diagrams are preserved for RPA-G R .Namely: (i) The critical temperature T c in all orbital channels is suppressed as density increases (r s decreases) and eventually becomes zero at the quantum transition point (QTP) r (ℓ) c , below which the 2D UEG no longer superconducts.(ii) On the logarithmic scale, −1/ ln[T c (r s )/T F ] still demonstrates a linear dependence on r s near the QTP, indicating that the BCS transition temperature obeys T (ℓ) c = T f e −1/λ ℓ with the coupling parameter λ ℓ = (r s − r (ℓ) c )/c ℓ , where c ℓ is a constant.(iii) The data for 1/ℓ versus r c in the inset of Fig. 1 suggests that superconductivity survives in the high-density limit for large enough ℓ.The qualitative agreement between RPA-G 0 and RPA-G R phase diagrams reflects the robustness of dynamic-screeningdriven superconductivity in the 2D UEG.

SUPERCONDUCTING PHASES OF THE TWO-DIMENSIONAL ELECTRON GAS WITH YUKAWA-TYPE INTERACTION
To investigate the impact of interaction range on the dynamic-screening-driven superconductivity in two dimensions, we consider the 2D electron gas with the Yukawa-type interaction V Y (q) = 2πe 2 / q 2 + q 2 0 , also known as a statically screened Coulomb interaction, where q 0 is the screening momentum.Specifically, we examine the behavior of the superconducting state in the 2D electron gas at r s = 0.8 when all channels exhibit superconductivity in the Coulomb system.We performed systematic and extensive calculations of the linear response function R for various channels ℓ and values of q 0 using RPA-G 0 and RPA-G R , respectively.
For each ℓ and q 0 , we determined the transition temperature T c by extrapolating the precursory Cooper flow to 1/R 0 (T c ) = 0. Subsequently, we used extrapolation to determine the critical values of q 0 (q 0c ), where −1/ ln[T c (q 0 )/T F ] becomes zero and superconductivity vanishes.Finally, we obtained the phase diagrams in the 2D electron gas with the Yukawa-type interaction under the RPA-G 0 [Fig.2(a)] and RPA-G R [Fig.2(b)] approximations, respectively.The qualitative agreement between the two phase diagrams indicates that the approximations are controlled.
As the range of interaction decreases (i.e., as q 0 increases), the superconducting critical temperature T c decreases for any given channel ℓ and eventually drops to zero at the QTP with the critical screening momentum q (ℓ) 0c .Since the rate of T c drop slows as ℓ increases, the dominant superconducting channel, determined from the crossing point of T (ℓ) c (q 0 ) and T (ℓ+1) c (q 0 ), increases as q 0 increases.For q 0 > 0.2k F , all orbital channels shown in Fig. 2 are no longer superconducting!These observations suggest that the long-range character of the Coulomb interaction is essential in dynamic-screening-driven super- TABLE II.Critical screening momentum q0c (in units of the Fermi momentum kf) of the 2D electron gas with a Yukawatype interaction under the RPA-G0 and RPA-G r approximations at rs = 0.4 and 0.8.Note that s-wave superconductivity is absent in the 2D electron gas with RPA-G r at rs = 0.8 and q0 = 0.
conductivity and that decreasing the interaction range suppresses superconductivity in the 2D electron gas.Furthermore, for q 0 → q − 0c , a linear relation exists between the inverse logarithmic scale −1/ ln[T c (q 0 )/T F ] and q 0 , indicating that the superconductivity near the QTP can be described by the emergent BCS theory with a coupling parameter λ ′ ℓ ∝ (q 0c − q 0 )/k f ≪ 1.

DATA TABLE OF QUANTUM TRANSITION POINTS
Results for quantum transition points r and left boundaries of the dominant superconducting channel r st (ℓ → ℓ + 1) for various orbital channels ℓ are given in Tab.I.The data for dimensionless critical screening momentum q (ℓ) 0c at r s = 0.8 are given in Tab.II.
2. Ultralow-temperature phase diagram of a 2D electron gas with rs = 0.8 and a Yukawa-type interaction VY(q) = 2πe 2 / q 2 + q 2 0 , calculated using the (a) RPA-G0 and (b) RPA-G R approximations.The superconducting temperature (Tc) as a function of the screening momentum (q0) is represented by lines going through data points for different channels, with the colored shadow indicating the dominant channel.For any given channel ℓ, T (ℓ) c decreases with increasing q0 until it disappears at the critical value q (ℓ) 0c .

DATA TABLE OF SUPERCONDUCTING TEMPERATURES OF THE TWO-DIMENSIONAL ELECTRON GAS
All the superconducting transition temperatures T  (in units of the Fermi temperature TF) for different rs from RPA-G0.The '-' symbol means that the Tc value is lower than the limit of double-precision floating-point number resolution.

8 FIG. 1 .
FIG.1.Finite-temperature flows of the largest s-wave eigenvalue, λmax, of the standard gap function equation (blue squares connected by line) and linear response functions χ0 (green triangles) and R0 (red circles) computed for the 2D UEG at rs = 1.156.The 1/χ0 and 1/R0 data are fitted using Eqs.(2)-(3), demonstrating excellent agreement.The role of the logarithmic numerator in Eq. (2) is clearly visible as the difference between the solid green line [fitted with ansatz (2)] and the dashed green line [fitted with ansatz (1)].The inset shows how linear in ln T scaling of 1/R0 is extrapolated to extract Tc in the f -channel at several densities.All 1/χ0 and 1/R0 flows undergo qualitative changes near T /TF ∼ 10 −1 , where TF is the Fermi temperature; this sets an energy scale Λ below which the attractive Cooper-channel interaction emerges.

T * /T f 10 − 5 T * /T f 10 FIG. 2 .
FIG.2.Temperature evolution of the standard pair susceptibility χ0 and modified pair susceptibility R0 of the 2D uniform electron gas in the s-channel for various values of rs.Red circles correspond to QTP rs = 0.6339, squares stand for stable regimes (g > 0), and triangles are used for the unstable regimes (g < 0).The lines are the fits with the ansatz (2) for χ0 and ansatz (3) for R0.(a) Function χ0(T ).For stable regimes, χ0(T ) saturates to a constant at T ≪ T * ≃ Λe −1/|g| ; for unstable regimes, χ0(T ) diverges at T = Tc; at the QTP, χ0(T ) diverges as T → 0. (b) Inverse χ0(T ) rescaled with the ideal-gas logarithmic factor.(c) Inverse R0.

FIG. 3 .
FIG.3.Superconducting phase diagram of the 2D UEG.For each channel, the line starting at QTP shows the (would-be) critical temperature.Critical values of rs for ℓ as large at 10 are presented in the inset.

FIG. 1 .
FIG. 1. Phase diagram of superconducting states in the 2D UEG within the RPA-G R framework showing s−, p−, d−, and f −wave superconducting states.Each orbital channel features a QTP with T (ℓ) c (rs → r (ℓ) c ) → 0. Critical values of rs for ℓ as large as 10 are presented in the inset.
in Tab.III with RPA-G 0 and Tab.IV with RPA-G r .

TABLE I .
The critical rs parameters for the 2D electron gas under the RPA-G0 and RPA-G r approximations determine the disappearance of superconductivity in channel ℓ at rc and the transition from the dominant channel ℓ to ℓ + 1 at rst.