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 $\Delta(\omega)$ must change sign at some (imaginary) frequency $\omega_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, $\Delta$ and $\omega_0$ must vanish in a correlated way such that $1/|\log[\Delta(0)]| \sim \omega_0^2$. We discuss the behavior of phase fluctuations near this transition and show that the correlation between $\Delta(0)$ and $\omega_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 towards a comprehensive scenario of superconductivity for a given material. In strongly correlated materials, like cuprates, iron-based, 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 orderparameter. 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 * dpimenov@umn.edu 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][8][9][10][11][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 sign-changing, 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 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/| log ∆| ∼ ω 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 (Abrikosov-Gorkov theory, [13][14][15][16]). In this situation, the destruction of superconductivity occurs via pair-breaking 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 = ω + iω , 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 pa iring 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Ṽ (ω − ω ) 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 odd-frequency ∆(ω). In this communication, we focus on the even-frequency 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 right hand side of Eq. (5) 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 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 ω, in f , we obtain 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 We assume and then verify that the transition is continuous, i.e., at f = f c − 0 + , ∆(ω) 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 signpreserving 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 nearperfect 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 where 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 artefact 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 .
We see that ω 0 vanishes as √ f c − f , 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.  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 mo- mentum 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.
Low-energy fluctuations around the mean-field solution correspond to slow variations of the phase of the order parameter θ(q): The expressions in the square brackets arise from small-θ expansion and Fourier transformation of the real-space phase factor 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): The expression for Π ∆ (q) in Eq. (24) can directly be obtained from the expansion in θ.
Multiplying both sides of Eq. (4) by V −1 (ω −ω)∆(ω) and integrating over ω,ω we find 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. We find where the derivatives are with respect to ω. A similar expression for n s was also obtained in Ref. [24].
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/ √ n s 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 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 |q|/∆(0). The form of the scaling function depends on the dimensionality. In 2D we have where we have approximated ∆(ω) ∆(0) and Λ/∆(0) = ∞, which holds to a good numerical accuracy. Evaluating the frequency integral we find: In 3D we have Figure 6. Generalized superfluid density. a n s (q) determined from Eqs. (35), (36) b κ(q), obtained in the same way as n s (q).
We plot n 2D s (q) and n 3D s (q) in Fig. 6a. We see that both functions decrease with increasing q, 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 Ω=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 superconductor-normal 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 signpreserving ∆(ω) 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. There are multiple possibilities to extend our analysis. One extension, which we leave for further research, is a potential co-existence of even-frequency and odd-frequency superconducting orders at f ≤ f c . Such a state spontaneously breaks time-reversal 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.

METHODS
Derivation of the expression for the superfluid density 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., [28]): Here, the first argument of ∆ contains the relative energy, and the second the total energymomentum 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 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 Tr ln −G −1 0 + Tr (G 0 · X) − 1 2 Tr (G 0 · X · G 0 · X) .
The first term is independent of θ, and the second term O(θ) cancels. To evaluate the third term O(θ 2 ), it is convenient to introduce center-of-mass coordinates as (Ω, q) = q = p − p and (ω, k) = k = p+p 2 . In these coordinates, where tr acts in the spinor space. After straightforward algebra, the combination of this term and (41) yields S θ from the main text, Eq. (24).
Derivation of the expression for the superfluid density from the particle-particle sus- ceptibility. An alternative way of deriving the expressions for the superfluid density n s and dynamical compressibility κ is by computing the full particle-particle susceptibility χ from Feynman diagrams. The basic building blocks for the diagrams are the normal and anomalous Green's functions, where ∆(ω) is chosen as real, α, β are spin indices, and σ y is a Pauli matrix.
We now compare the prefactors for O(1), O(|q| 2 ), and O(Ω 2 ) on both sides of this equation.
(54). Cancelling the identical terms on both sides, we solve for c 1 : In a similar fashion we obtain Combining (55) , (56) and (48), we obtain, to the leading order in Ω, |q|, where Code availability. The codes used to generate the numerical data are available upon request from the corresponding author.