Quantum phase transition in a clean superconductor with repulsive dynamical interaction

We consider a model of electrons at zero temperature, with a repulsive interaction which is a function of the energy transfer. Such an interaction can arise from the combination of electron–electron repulsion at high energies and the weaker electron–phonon attraction at low energies. As shown in previous works, superconductivity can develop despite the overall repulsion due to the energy dependence of the interaction, but the gap Δ(ω) must change sign at some (imaginary) frequency ω0 to counteract the repulsion. However, when the constant repulsive part of the interaction is increased, a quantum phase transition towards the normal state occurs. We show that, as the phase transition is approached, Δ and ω0 must vanish in a correlated way such that 1/∣log[Δ(0)]∣~ω02\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1/| \log [{{\Delta }}(0)]| \sim {\omega }_{0}^{2}$$\end{document}. We discuss the behavior of phase fluctuations near this transition and show that the correlation between Δ(0) and ω0 locks the phase stiffness to a non-zero value.


INTRODUCTION
Understanding the nature of the "pairing glue", which enables Cooper pair formation of fermions, is one of the key steps toward a comprehensive scenario of superconductivity for a given material. In strongly correlated materials, like cuprates, ironbased, heavy-fermion, and organic materials, the attractive pairing interaction is likely of electronic origin. Near a quantum phase transition, such an attraction often takes a more concrete form of an effective four-fermion interaction, mediated by soft collective fluctuations of the corresponding order parameter. Most often, the attraction emerges in a channel different from an ordinary s-wave, in which case the superconductivity is labeled as an unconventional one.
For more conventional metals the symmetry of the pairing gap is s-wave, and the attraction is believed to come from electron-phonon interaction. This is the backbone of the "conventional" BCS theory of superconductivity. Still, to fully understand the phononic mechanism of s-wave superconductivity, one must explain why it is not overshadowed by the Coulomb repulsion, which is seemingly much larger. The frequently cited explanation [1][2][3][4][5][6] is that the repulsive Coulomb repulsion is logarithmically renormalized down between the Fermi energy E F and the Debye energy Ω D (the Tyablikov-McMillan logarithm), and at energies below Ω D becomes smaller than the electron-phonon attraction, if the ratio E F /Ω D is large enough.
Upon closer examination, this explanation appears somewhat incomplete as Tyablikov-McMillan renormalization holds for the full interaction, i.e., for the sum of electron-electron and electron-phonon interactions, and under the renormalization this full interaction decreases, but does not change sign. It has been realized by several authors 4,7-12 that the underlying reason why electron-phonon superconductivity holds despite larger Coulomb interaction, is that the full interaction on the Matsubara axis (where it is real) is a dynamical one, V(Ω m ), and although a phonon-mediated attraction does not invert the sign of V(Ω m ), it nevertheless reduces it at frequencies below the Debye energy. It was argued that an "average" repulsive V(Ω m ) can be effectively eliminated from the equation for the pairing gap Δ(ω m ), by choosing a solution which changes sign as a function of ω m . This bears some similarity to how, for an electronic pairing, a static Coulomb repulsion is effectively eliminated by choosing a signchanging, non-s-wave spatial structure of the gap function.
A convenient way to model the dynamical V(Ω m ), suggested in refs. [8][9][10][11][12] , is to treat it as a sum of two parts: a constant repulsive part of strength f, representing the renormalized instantaneous Coulomb repulsion, and a frequency-dependent attractive part, due to electron-phonon interaction: where Ω 1 is of order of the Debye energy. A similar reasoning has been applied 10 to dynamically screened electron-electron interaction, where Ω 1 is of the order of plasma frequency. For f > 1, V(Ω m ) > 0 for all frequencies, yet for 1 < f < f c , superconductivity emerges below a finite T c , which contains f in the combination f =½1 þ const: f logðE F =Ω 1 Þ. For a given f and large enough logðE F =Ω 1 Þ, the Coulomb repulsion becomes logarithmically small, and one recovers the McMillan formula for T c . One the other hand, at a given E F /Ω 1 , at large enough f > f c , the repulsion becomes too strong and superconductivity vanishes. Obviously, T c and the magnitude of the gap Δ(ω m ) vanish at f = f c .
It is the goal of the present work to understand the nature of the T = 0 quantum phase transition between a superconducting state at f < f c and a normal state at f > f c . Specifically, we resolve the following puzzle: on the one hand, the gap Δ(ω m ) must change sign at some finite ω m = ω 0 , otherwise there would be no solution of the gap equation for f > 1. On the other hand, for any finite ω 0 , Δ(0) is non-zero, in which case the linearized gap equation does not have a solution as the pairing kernel contains an infrared-divergent Cooper logarithm, which is not regularized at T = 0 and therefore does not admit a solution. We argue analytically and check numerically that as f approaches f c from below, ω 0 → 0 and Δ → 0 in tune with each other, such that 1=j log Δj $ ω 2 0 . We also analyze the spectrum of gapless phase fluctuations near f = f c . We show that because of the relation between Δ and ω 0 , the superfluid stiffness remains finite as f approaches f c from below. This is in marked contrast with the behavior of the stiffness near the end point of superconductivity at T = 0 in a system with magnetic impurities . In this situation, the destruction of superconductivity occurs via pairbreaking due to the impurity-induced self-energy, and the superfluid stiffness gradually vanishes as the system approaches the T = 0 phase transition.
That superconductivity vanishes when ω 0 = 0 can also be interpreted from a topological viewpoint, because ω m = ω 0 is a center of a dynamical vortex: the anti-clockwise circulation of the phase of Δ(z), z ¼ ω 0 þ i ω 00 , around this point is 2π, refs. [17][18][19] . There is no way to eliminate this dynamical vortex as there are no anti-vortices in the upper half-plane of frequency (their presence would be incompatible with the analyticity of Δ(z)). Hence, as long as superconducting order is present, ω 0 must remain finite. The only possibility for a vortex to disappear without destroying superconductivity is when it moves to an infinite frequency. For the model of Eq. (1) this holds at f = 0, and for f < 0 the gap function Δ(ω m ) on the Matsubara axis is nodeless.
We note in passing that a vortex on the Matsubara axis gives rise to 2π winding of the phase of Δ(ω) on the real frequency axis, between ω = −∞ and ω = ∞. Such phase winding necessary leads to nodes in the real and imaginary parts of Δ(ω), which can be detected by spectroscopic experiments, e.g., ARPES 20 .

Model
We consider a spatially isotropic model of interacting spin-1/2 fermions at zero temperature in d dimensions, described by the effective low-energy action where Λ is a UV cutoff of order E F and ω are Matsubara frequencies (here and below we label Matsubara frequency as ω without subscript m). The interaction V(Ω) is taken to be a function of the energy transfer, but independent of momenta. We follow earlier works [9][10][11][12]21 and set V(Ω) to where ρ is the single-spin density of states at the Fermi surface, and Ω 1 is of the order of Debye energy for the electron-phonon case (the factor of 2 is introduced for notational convenience). In the following, we measure all energies in units of Ω 1 , and hence set Ω 1 = 1 in Eq. (3). Then, Λ ≫ 1. A discussion of the opposite low-density limit where E F , Λ ≪ 1 can be found in ref. 22 . For the known physical realizations of Eq. (3), f > 1, hence V(Ω) remains positive (repulsive) at all frequencies. For completeness, here we consider arbitrary f, but our key focus will still be on f > 1. For a generic f, V(Ω) is purely attractive for f ≤ 0, is attractive at small frequencies and repulsive at large frequencies for 0 < f < 1, and is purely repulsive for f ≥ 1, see Fig. 1a. The dimensionless λ parametrizes the overall strength of the interaction. We assume λ ≤ 1, this will allow us to neglect, at least qualitatively, the normal fermionic self-energy: One can show that the leading self-energy effect is a mere renormalization of the coupling constant λ → λ/(1 + 2λ).

Gap equation
To describe superconductivity, we perform a Hubbard-Stratonovich transformation in the spin-singlet, s-wave pairing channel and use a saddle point approximation. This procedure leads to the conventional Eliashberg equation for the gap function 23 , though without the additional contribution from the self-energy. On the Matsubara axis, we have The interactionṼðω À ω 0 Þ is real on the Matsubara axis, which allows us to set Δ(ω) to be real by properly choosing its phase. At the same time, because the interaction is a function of the frequency transfer, one can search for even-frequency and oddfrequency Δ(ω). In this communication, we focus on the evenfrequency solutions. For even-frequency Δ(ω) = Δ(−ω), the gap equation can be rewritten as It is obvious that for f > 1, whenṼ > 0, Δ(ω) must change sign at some frequency ω 0 as otherwise the left-hand side and the righthand side of Eq. (5) would have opposite signs. The value of ω 0 is chosen to minimize the effect of a repulsive f in Eq. (3). This has been discussed before [2][3][4][5][6]8,11,12 and we just state the results.
has a node as long as ω 0 < Λ. Second, optimizing ω 0 in the limit logðΛÞ ) 1, one obtains that a repulsive f effectively gets reduced to f =½1 þ λf logðΛÞ ! 1=λ logðΛÞ. This gives rise to the McMillan formula for T c / e À1=ðλÀμ Ã Þ , in which μ Ã % 1= log Λ is the contribution from the repulsion. Third, for any finite Λ, the repusive part of the interaction gets reduced, but cannot be completely eliminated. As a result, superconductivity exists at f smaller than some critical f c > 1.
In Fig. 1b we present the numerical solution of the non-linear gap equation (5) for some representative Λ and 1 < f < f c . We clearly see that Δ(ω) changes sign at some finite ω 0 . It reaches a finite value at ω = 0 and then saturates at some other finite value, of opposite sign at ω 0 ≪ ω < Λ. The numerical solution has been D. Pimenov and A.V. Chubukov obtained by a "damped iteration" method, in which only a certain portion of Δ(ω) is updated at each step of iterations. This method improves the convergence of the iteration procedure 12 .
Quantum phase transition towards a superconductor with nodeless Δ(ω) Before we proceed to the case f ≈ f c , we briefly discuss the transition towards the state with a nodeless Δ(ω). As stated above, this transition occurs at f = 0 when Λ = ∞, i.e., Eq. (5) holds at all frequencies. This transition can be classified as topological because it separates two states with and without a dynamical vortex. As f is reduced towards f = 0, ω 0 increases, i.e., the core of the dynamical vortex successively moves to larger ω. At f = 0 it reaches ω = ∞ and disappears.
We now argue that the dependence of ω 0 of f has a simple form This can be obtained as follows: Let Δ a (ω) be the solution of the gap equation at f = 0: SinceṼðωÞ is purely attractive for f = 0, Δ a has a fixed sign. At large Now let Δ b (ω) be the solution of the gap equation at small but finite f. To the leading order in f, we obtain, for large ω, At the node, Δ b (ω 0 ) = 0, hence ω 0 = f −1/2 . For large but finite Λ, the phase transition occurs when ω 0 = Λ. The corresponding critical f for a topological transition is then In Fig. 2 we checked these results by solving the gap equation numerically. The agreement between the numerical and analytical results is perfect.

Quantum phase transition towards normal state
We now consider the system behavior near the T = 0 transition towards the normal state at f = f c > 1, when the pairing interaction Vðω À ω 0 Þ is positive (repulsive) at all frequencies. We assume and then verify that the transition is continuous, i.e., at is infinitesimally small. Like we said in the Introduction, to understand this transition one has to resolve the following puzzle: if infinitesimally small Δ(ω) tends to a finite value at ω = 0, like, e.g., in Fig. 1b, the right-hand side of the linearized gap equation gives rise to a divergent Cooper logarithm. Because T = 0, the logarithmical divergence is not cut. The only way to avoid this divergence is to place the node (i.e., the vortex core) right at ω = 0. But then the gap becomes sign-preserving at all finite ω, and for such Δ(ω) there is no solution of the gap equation for a purely repulsive interaction.
As we now show, the resolution of this problem is to let both Δ and ω 0 vanish in a correlated way as f → f c from below. To simplify the analysis, we first note that for all f < f c , the gap function Δ(ω) is well approximated by a simple form A comparison with the numerical solution of the gap equation shows that this form is near-perfect for ω < 1 and matches the numerical results reasonably well for ω > 1, see Fig. 1b. Such agreement is sufficient to extract the leading behavior near the phase transition (see below). The coefficients Δ 0 , Δ 1 can be determined by inserting the ansatz (11) into (5) and expanding up to second order in ω. After a straightforward algebra we obtain In the limit Δ 0 , Δ 1 → 0, the equations simplify to The value of the critical f c can be determined by evaluating the determinant of the set (15) in the limit ℓ → ∞. We obtain The divergence of f c at a critical value of λL, which is evident from Eq. (17) is not an artifact of the approximation to Δ(ω), as we have checked numerically. Rather, it implies that by properly placing ω 0 , one can completely eliminate a constant repulsion f even when f is large. A detailed analysis of this effect will be presented elsewhere (in preparation).
The general trend that f c increases with increasing Λ is also in agreement with McMillan reasoning that the Coulomb repulsion is suppressed at large Λ. In the following, we focus on λL ≪ 1, in which case Evaluating the determinant again, but this time for a finite ℓ, we obtain to leading order in λL and f c − f: Using (16) we find that Δ 0 ¼ expðÀℓÞ vanishes exponentially fast as f ↗ f c . From the first equation in (15) we obtain The ratio is negative (hence ω 0 is finite) and progressively decreases when f approaches f c . Substituting Δ 0 /Δ 1 into (11), we obtain We see that ω 0 vanishes as ffiffiffiffiffiffiffiffiffiffiffiffi f c À f p , i.e., much more gradually than Δ 0 .
In Fig. 3, we verify the scaling forms of Δ 0 and ω 0 by extracting these two quantities from the numerical solution of the gap equation. The agreement between analytical and numerical results is quite good.
Phase fluctuations near critical f c For a more detailed characterization of the phase transition at f = f c , we now look at soft collective excitations in the system. These are phase fluctuations, which in the absence of long-range Coulomb interaction are Goldstone modes of the superconducting state. Our goal is to derive the superfluid density and the dynamical compressibility, which enter the propagator of phase fluctuations, as functions of q = (Ω, q), where q the total momentum of a Cooper pair and Ω is the total frequency. There are two ways to do this: either expand the action to second order in θ or analyze the pole structure of the full particle-particle susceptibility at small q. These two methods yield consistent results; we will focus on the first one in the remainder of this section as we discuss the other one in the "Methods" section.
To obtain the propagator of low-energy phase fluctuations, we introduce the total momentum q and the total frequency Ω of a Cooper pair. For convenience, we combine q and Ω into a (d+1)dimensional variable q = (Ω, q). In our mean-field solution the pairing involves fermions with frequencies ω and −ω and momenta k and −k, i.e., q is set to zero. In other words, the mean-field gap is a function of ω but not of q. This mean-field solution corresponds to a minimum of the Luttinger-Ward functional. States away from the minimum are described by a fluctuating pairing field (order parameter) that depends on both Ω and q. We illustrate this in Fig. 4.
The expressions in the square brackets arise from small-θ expansion and Fourier transformation of the real-space phase factor Z d dþ1 q ð2πÞ dþ1 exp½iθðxÞ exp½iðx Á qÞ ' where x is the center-of-mass coordinate. Inserting the expansion (22) into the Δ-dependent action, where the fermions have been integrated out, and expanding to the second order in θ, we obtain the following action for the θ field (see "Methods" for details):  Here, V À1 ðω À ω 0 Þ is the inverse of Vðω À ω 0 Þ defined by Z dω 1 2π Vðω À ω 1 ÞV À1 ðω 1 À ω 0 Þ ¼ 2πδðω À ω 0 Þ: The expression for Π Δ (q) in Eq. (24) can directly be obtained from the expansion in θ. Alternatively, it can be obtained diagrammatically as a particle-particle bubble with form-factors Δ 2 (ω) (see "Methods").
Multiplying both sides of Eq. (4) by V À1 ðω ÀωÞΔðωÞ and integrating over ω;ω we find that Z dωdω 0 ð2πÞ 2 ΔðωÞðÀV À1 Þðω À ω 0 ÞΔðω 0 Þ ¼ This expression coincides with the particle-particle bubble in the limit of vanishing q: As a result, the propagator of the field θ(q) is determined by Π Δ (q) − Π Δ (0). Expanding Π Δ (q) to second order in Ω, q, we obtain The coefficient n s in the second line of (29) is the superfluid density, normalized by the density of electrons in the normal state. It parametrizes the energy cost of spatial phase fluctuations of the order parameter and controls the supercurrent and magnetic response in the superconducting state (The prefactor 1/d comes from averaging over cos ½ðk; qÞ 2 in d dimensions). The factor κ is a dynamical compressibility which parametrizes the energy cost of temporal phase fluctuations. We find ΔðωÞ 2 ð3 À ΔðωÞΔðωÞ 00 Þ þ ω 2 ðÀΔðωÞΔ 00 ðωÞ þ 3Δ 0 ðωÞ 2 Þ h i : (31) where the derivatives are with respect to ω. A similar expression for n s was also obtained in ref. 24 .
For a weakly frequency-dependent Δ(ω), the expressions for n s and κ are the same as in BCS theory, n s = κ = 1. For a generic Δ(ω), the BCS results are correct by order of magnitude, but the actual values of n s and κ differ by BCS expressions by O(1), see Fig. 5.
At f ≤ f c , Δ(ω) is exponentially small, and the integrals are dominated by ω ≲ Δ(0). Because frequency variation of Δ(ω) occurs at a much larger scale ω 0 ≫ Δ(0), the gap in (30) and (31) can be approximated by a constant Δ(0). As a result, both n s and κ tend to BCS values n s = κ = 1. We verified this result in numerical calculations, see Fig. 5. The velocity of phase fluctuations also approaches the BCS value At f = f c + 0 + , n s and κ jump to zero. We emphasize that the behavior of n s near the T = 0 superconductor-normal state phase transition in a clean system at f = f c is different from the one at the T = 0 superconductornormal state phase transition due to magnetic impurities. There, n s gradually vanishes at a critical impurity concentration due to pair-breaking coming from the impurity-induced self-energy. As a result, the penetration depth $ 1= ffiffiffiffi ffi n s p diverges (refs. 13,14,16 ). At a technical level, this is because in the case of magnetic impurities the denominator of Eq. (30) contains an additional term proportional to the fermionic damping rate due to impurity scattering. In our case, such a constant term is absent. We expect, however, that it will appear if we extend the analysis of the superconductor-normal state phase transition to a finite magnetic field H. We therefore expect that at a finite field, n s will vanish at critical f c (H).
Still, the discontinuity of n s at f = f c in our case holds only for superfluid density evaluated at zero momentum q, or, more accurately, at |v F q| ≪ Δ 0 . To analyze the behavior at larger |q|, we define a generalized momentum-dependent superfluid density as This n s (q) is a scaling function ofq ¼ v F jqj=Δð0Þ. The form of the scaling function depends on the dimensionality. In 2D we have n 2D s ðqÞ ' (34) where we have approximated Δ(ω) ≃ Δ(0) and Λ/Δ(0) = ∞, which holds to a good numerical accuracy. Evaluating the frequency integral we find: n 2D s ðqÞ ¼ In 3D we have n 3D s ðqÞ ¼ 12 ffiffiffiffiffiffiffi ffi ffiffiffiffiffiffiffi ffi We plot n 2D s ðqÞ and n 3D s ðqÞ in Fig. 6a. We see that both functions decrease with increasingq, i.e., with decreasing Δ 0 for a given q.
At f = f c − 0 + , n 2D;3D s ðqÞ vanishes for any finite q. A suitably defined momentum-dependent compressibility κðqÞ $ ∂ 2 Π Δ ðΩ; qÞ=∂Ω 2 j Ω¼0 follows the same trend (Fig. 6b). This behavior is indeed fully expected as for |v F q| ≫ Δ(0), the system is effectively in the normal state, where the U(1) symmetry is preserved and gauge (phase) fluctuations do not cost any energy.

DISCUSSION
In this communication, we analyzed a T = 0 superconductornormal state transition for a model of fermions coupled by a frequency-dependent interaction V(Ω), which has a repulsive constant part f and an Ω-dependent attractive part. For f < 0, V(Ω) is fully attractive, and the system displays a conventional s-wave superconductivity with a sign-preserving Δ(ω) along the Matsubara axis. At 0 < f < 1, the interaction V(Ω) is attractive at small frequencies, but repulsive at large Ω. In this case, superconductivity is still present at T = 0, but the gap function has a node along the Matsubara axis, at some finite ω = ω 0 . Because a nodal point of Δ(ω) is a center of a dynamical vortex, the superconducting states at f < 0 and at f > 0 are topologically different. We analyzed the topological transition at f = 0 and argued that the vortex emerges at an infinite frequency at f = 0 and moves to a finite ω 0 at a finite f. This is an expected behavior, consistent with earlier analysis of a similar model 17 . We also analyzed how the critical f for such topological transition changes if we set a finite UV cutoff for the interaction.
The superconducting state with a sign-changing Δ(ω) persists also for f > 1, when V(Ω) becomes positive at all frequencies, and vanishes at a finite f c . The key part of our work is the analysis how the gap function Δ and the frequency ω 0 , where Δ changes sign, behave at f ≤ f c .
We found that the gap function at zero frequency, Δ(0), vanishes exponentially fast with (f c − f). The frequency ω 0 vanishes as well, but parametrically slower as ffiffiffiffiffiffiffiffiffiffiffiffi f c À f p . We argued that this parametrical difference between Δ(0) and ω 0 allows one to obtain a non-zero solution of the gap equation for all f < f c . We note that the transition at f c can be also interpreted from topological perspective, as the center of the dynamical vortex reaches ω = 0 at f = f c and would have nowhere to go if superconductivity persisted above f c .
We complimented the analysis of Δ(ω) near f c by the analysis of the propagator of phase fluctuations. We have shown that the superfluid density and the compressibility, which control momentum and frequency parts of the propagator of a phase field, deviate from the BCS values at a generic f < f c , but tend to the BCS values at f → f c and undergo a finite jump at f = f c + 0 + . We showed that this, however, holds only for the superfluid density (and the compressibility) defined at strictly zero momentum. We introduced a generic momentum-dependent n s (q) and showed that it gradually vanishes at f → f c for all |v F q| > Δ 0 . At f slightly below f c , this behavior holds or all |q| except the ones which are exponentially small in f c − f.
There are multiple possibilities to extend our analysis. One extension, which we leave for further research, is a potential coexistence of even-frequency and odd-frequency superconducting orders at f ≤ f c . Such a state spontaneously breaks timereversal symmetry.
Finally, our analysis is not constrained to electron-phonon interaction and is applicable to all cases when there is a near-constant repulsion and frequency-dependent, retarded attraction due to a boson exchange. Other interesting candidates for a boson are exiton-polaritons in a microcavity 25,26 or cavity photons 27 . Experimental cavity setups often come with a tuning knob which allows one to change the relative strength of repulsive and attractive components of the interaction (i.e., continuously change f in our model). This should allow one to observe phase transition at f c that we analyzed in this work.

Derivation of n s by expanding the action in θ
After the Hubbard-Stratonovich transformation is performed, the mixed boson-fermion action takes the form (for the following derivation, compare, e.g., ref. 28 ): (37) Here, the first argument of Δ contains the relative energy, and the second the total energy-momentum of the Cooper pair, compare Fig. 4. Integrating out the fermions, we obtain a purely bosonic action where the trace runs over energy-momenta and Nambu indices. From Eq. (38), the gap equation is simply derived by setting δS=δΔðωÞ ¼ 0. We look for mean-field solutions which have zero total energy-momentum, i.e., contain a delta-function δ (d+1) (q).
To find the action of the Goldstone-mode, we insert the expansion from Eq. (22) into (38). The Oðθ 0 Þ contribution from the S Δ term has the form Z dωΔðωÞV À1 ðω À ω 0 ÞΔðωÞ δ ðdþ1Þ ð0Þ : The term δ (d+1) (0) should be interpreted as volume factor. The OðθÞ-contribution cancels, and the Oðθ 2 Þ-contribution reads To expand the ðTr lnÞ-term, it is convenient to split Now, the trace of the logarithm can be expanded as Fig. 6 Generalized superfluid density. a n s (q) determined from Eqs.

DATA AVAILABILITY
The numerical data used in the analysis in this work are available upon request from the corresponding author.

CODE AVAILABILITY
The codes used to generate the numerical data are available upon request from the corresponding author.