Induced odd-frequency superconducting state in vertex-corrected Eliashberg theory

We show that vertex corrections to Migdal's theorem in general induce an odd-frequency spin-triplet superconducting order parameter, which coexists with its more commonly known even-frequency spin-singlet counterpart. Fully self-consistent vertex-corrected Eliashberg theory calculations for a two dimensional cuprate model, isotropically coupled to an Einstein phonon, confirm that both superconducting gaps are finite over a wide range of temperatures. The subordinate $d$-wave odd-frequency superconducting gap is found to be one order of magnitude smaller than the primary even-frequency $d$-wave gap. Our study provides a direct proof of concept for a previously unknown generation mechanism of odd-frequency superconductivity as well as for the generic coexistence of both superconducting states in bulk materials.


I. INTRODUCTION
Odd-frequency superconductivity describes a Cooper pair condensate in which the quantum mechanical wave function is odd under exchange of relative time between two electrons of a Cooper pair. As first proposed by Berezinskii [1], this ordered state is compatible with the Pauli principle as long as either the momentum space parity is even and the electrons are of spin-triplet type, or the parity is odd and one considers spin-singlet electron pairs, see e.g. Ref. [2] for a recent review. This type of superconductivity is clearly beyond the Bardeen-Cooper-Schrieffer (BCS) theory that neglects time-retardation effects between the electrons of the Cooper pair, but can be described naturally within Eliashberg's theory [3]. Historically, it has been argued that electron-phonon coupling can theoretically lead to odd-frequency superconductivity, provided that the interaction has a sufficiently strong odd momentum component as a result of coupling between electrons and acoustic phonons [4]. However, such a scenario was considered unlikely if no spin dependence enters into the interaction [5]. More recent theories have shown that signatures of an odd-frequency phase are indeed possible due to electron-phonon interactions only, depending on the interaction strength [6][7][8][9].
Odd-frequency superconductivity has been discussed recurrently in various forms [2,10]. In this discussion it is important to distinguish between existence of a superconducting condensate with an odd-frequency order parameter, i.e., with nonzero gap function, and oddfrequency pair correlations. The latter can occur when a frequency-even superconductor is placed in close proximity to a magnetic layer or when a magnetic impurity is placed on a conventional superconductor [10][11][12][13][14][15][16][17][18]. The spins of the conventional Cooper-pair electrons become rotated in the exchange field, leading to a spin-triplet component with time-odd parity [10,19]. Recent exper- * fabian.schrodi@physics.uu.se † alex.aperis@physics.uu.se ‡ peter.oppeneer@physics.uu.se iments have provided spectroscopic observations of oddfrequency pair correlations [17,[20][21][22][23]. The observation of an odd-frequency superconductivity phase in a bulk material remains however thus far elusive. The thermodynamic stability of an odd-frequency phase in the presence of an even-frequency order parameter has therefore been studied theoretically in several works [8,9,[24][25][26]. In many of the theoretical proposals for systems exhibiting this rather exotic type of Cooper pairing, there exists one particular degree of freedom that makes the occurrence of an odd-frequency superconducting pairing possible. As laid out in Ref. [2], odd-frequency superconductivity in bulk materials is often associated with spin degrees of freedom, or with the strong-coupling limit of electron-phonon interactions [6,26]. Other examples are pairing effects due to multiple band or electron orbitals [27,28], influence of magnetic fields [7,8], magnetic impurities [18], proximity effects [29,30], or due to edge-state induced ferromagnetic fluctuations [31].
Here we show that this conception is generally not complete. By deriving self-consistent and vertex-corrected Eliashberg equations, we prove that a finite oddfrequency spin-triplet order parameter can theoretically always coexist along with the 'standard' even-frequency spin-singlet superconducting gap. Vertex corrections to the electron-phonon problem generally become important if the ratio of phonon to electron energy scale is not small [32], and it recently has been shown that they are a potential candidate for producing unconventional symmetries in several classes of superconductors [33]. Our current results suggest that the phenomenology of vertex corrections is even richer, in that they can also be responsible for odd-frequency superconducting states. We benchmark our theory by solving the Eliashberg equations for a Holstein-like Hamiltonian with Einstein phonon and isotropic electron-phonon coupling, resembling the characteristics of a prototypical cuprate superconductor. The odd-frequency d-wave state is induced by its even-frequency counterpart, which is the primary order parameter, and is approximately one order of magnitude smaller in amplitude. Our results suggest that the coexistence of both condensates is possible due to a mixing of momentum space representations, triggered by the inclusion of vertex corrections into our Eliashberg framework.
In the following we assume for simplicity that the spatial parity is always even, hence we want to describe the most common Brillouin zone (BZ) symmetries of the superconducting gap function, such as s-wave or d-wave. We thereby leave the straightforward generalization to more exotic p-wave or f -wave states for future analysis. Consequently, we are left with spin and time (frequency) as the only two degrees of freedom, since we do not include multiple orbitals or bands in our formalism. Then the possible superconducting order parameters are either even in frequency and spin singlet, or odd in frequency and spin triplet.
The remainder of this paper is organized as follows: We start Section II by introducing the general mathematical framework of our theory, and present the vertexcorrected Eliashberg equations, as they were developed in Ref. [33], in Section II A. We show in Section II B that the results to those self-consistent equations always induce an odd-frequency order parameter, which describes Cooper pairs of spin-triplet electrons. The self-consistent set of nonadiabatic Eliashberg equations, including both even-and odd-frequency order parameters, is given in Section II C. We test our theory by using a one-band tight-binding description, which describes a prototypical cuprate superconductor. Section III contains an analysis of our numerical results with respect to momentum and temperature dependencies. Finally, we conclude with a short summary and discussion in Section IV.

II. VERTEX-CORRECTED ELIASHBERG THEORY
Let us assume that the atomic vibrations in the superconductor can approximately be described by a singlebranch optical phonon mode with frequency Ω. The purely electronic degree of freedom is characterized by a single-band dispersion ξ k , with k a BZ wave vector. We consider a momentum independent electron-phonon scattering strength g 0 , hence, the Hamiltonian of the system can be expressed as with phonon displacement u q = b † q + b −q , Pauli matriceŝ ρ i and Nambu spinor Ψ † k = c † k,↑ , c −k,↓ . We use b † q and c † k,σ as phonon and electron creation operators, with σ ∈ {↑, ↓} labeling the electron spin.
At temperature T we define fermion and boson Matsubara frequencies as ω m = πT (2m + 1) and q l = 2πT l, respectively, with m, l ∈ Z. For sake of brevity, we adopt from here on the four-momenta notation k = (k, iω m ) and q = (q, iq l ). For the Einstein phonon spectrum considered here the phonon propagator can be written as a simple Lorentzian in frequency space, D q = 2Ω/(Ω 2 +q 2 q ). The electron Green's function obeys the Dyson equation G k =Ĝ 0 k +Ĝ 0 kΣ kĜk , which equivalently can be written asĜ In Eq. (2),Σ k is the electron self-energy, and the bare electron Green's function is defined asĜ 0 k = [iω kρ0 − ξ kρ3 ] −1 . The electron self-energy including the two lowest order Feynman diagrams for electron-phonon scattering can be expressed aŝ where the electron-phonon interaction is defined via V q = g 2 0 D q [33][34][35][36]. The first-order vertex corrections are given by the second term in Eq. (3). In the following Section II A we derive self-consistent Eliashberg equations describing the interacting state of our model system, where we follow closely the notation of Refs. [33,37]. We show in Section II B that our theory allows for an induced subordinate odd-frequency spintriplet superconducting state. A complete mathematical framework for the coexistence of even-and odd-frequency Cooper pairs is given in Section II C.

A. Even-frequency, spin-singlet Cooper pairs
The vast majority of superconductors exhibits the formation of spin-singlet Cooper pairs, where all quantities transform even along the real and Matsubara frequency axis. Within our formalism, such a state can conveniently be described by the electron self-energy ansatẑ where the electron mass renormalization Z k , chemical potential renormalization χ k and superconducting order parameter φ k are introduced. Plugging Eq. (4) into Eq. (2) givesĜ with determinant As short-hand notation we introduce the symbols so that the electron Green's function can be written aŝ As has been shown in Ref. [33], Eqs. (3) and (8) can be used to derive a set of self-consistent vertex-corrected Eliashberg equations, reading For brevity, the above equations are defined in terms of the vector γ T k = γ The vertex-corrected Eliashberg Eqs. (9)-(11) exhibit richer phenomenology and more degrees of freedom than the standard theory, which only takes into account the lowest order Feynman diagram. For example, it has been shown by the current authors that superconductivity can be partially suppressed by vertex corrections [33], and can eventually lead to unconventional BZ symmetries of the pairing function [37,38].

B. Induced odd-frequency spin-triplet state
We now want to take a closer look at the electron selfenergy. The Pauli matricesρ i , i ∈ {0, 1, 2, 3}, form a basis of the 2 × 2 Nambu space employed here. Notice, however, that the ansatz forΣ k in Eq. (4) does not includeρ 2 . The neglect of this channel is justified due to the gauge freedom of Eliashberg theory, i.e. any function proportional toρ 2 describing spin-singlet, even-frequency Cooper pairs can be interpreted as a complex phase of the superconducting gap function, which will differ from φ k only by a proportionality factor. We are allowed to set the prefactor ofρ 2 identically zero, since no experimentally observable effects arise from this complex phase, except in Josephson tunneling [39] which is of no interest here. This is the well-established picture for standard Eliashberg theory [40,41], and remains valid even upon the inclusion of vertex corrections.
Considering the self-energy expression of Eq. (3), we can defineΣ k =Σ correspond to contributions from the first and second order Feynman diagrams, respectively. From our ansatz Eq. (4) it follows that the electron Green's function does not include any term proportional toρ 2 , therefore we obtain from Eq. (3), as expected. However, applying the same operation to the non-adiabatic correction yields which does generally not vanish. Due to the fact that ζ ind k ∝ρ 2 , even though our definition ofΣ k does not include theρ 2 channel, Eq. (16) describes an induced state. Similar induction mechanisms due to magnetic fields have been previously proposed within purely BCS pictures [42,43], however here retardation effects, that are captured by Eliashberg theory, are crucial for discussing induced odd-frequency superconductivity [9].
A closer inspection of the induced function ζ ind k reveals that each term in Eq. (16) contains γ (Z) , γ (χ) and γ (φ) exactly once, with different combinations of momentum and frequency dependencies. Consequently, ζ ind k = 0 in the normal state, because then φ k and therefore γ (φ) identically vanishes, compare Eq. (7). We can therefore safely conclude that ζ ind k describes physics of the superconducting state.
Next, we want to examine the frequency symmetry of ζ ind k , which we do here by neglecting the momentum dependence for the moment. We can safely assume that Z m is even along the Matsubara frequency axis, therefore Z m = Z −m−1 , and likewise for χ m and φ m . From the definitions in Eq. (7) it directly follows that γ For the sake of the argument, let us define which is representative for each term in Eq. (16). Note, that we write the definition of m 3 = m 2 − m 1 + m here explicitly. The function β m has fermionic Matsubara frequency symmetry, so we need to consider β −m−1 . Replacing m with −m − 1 in Eq. (17) m2−m1−m−1 . We therefore make the definitionm 1 = −m 1 − 1, or, equivalently, m 1 = −m 1 − 1. Since the summation over index m 1 includes all integers, we alternatively sum over indexm 1 , so as to get Next, we make a similar definition as before,m 2 = −m 2 − 1, leading to As a last step we invert the frequency index of each function on the second line of Eq. (19), using the aforementioned symmetries. After renaming the dummy indices m 1 andm 2 to m 1 and m 2 , respectively, we arrive at From here it is trivial to show that all terms in Eq. (16) obey this particular symmetry, hence ζ ind k is an oddfrequency pairing function. Note, that this does not imply a breaking of time reversal symmetry [2,44].
When focusing on the momentum dependence, it is more difficult to make concrete analytical statements. This is due to the fact that the BZ symmetries of φ and χ are a priori not known. We return to this aspect in Section III A when presenting our numerical results.
As mentioned before, we assume even momentum space parity and neglect orbital/band degrees of freedom in this work, hence we can restrict our symmetry classification of ζ ind k to spin and frequencies. We have shown analytically that the induced function is odd along the Matsubara frequency axis, it therefore describes spin triplet electron pairs [2] which coexist with the even-frequency, spin singlet Cooper pairs described by φ k . Needless to say, all this interpretation is needed only if ζ ind k is finite, which can not be guaranteed from the functional form of Eq. (16). Therefore, we derive in the following Section II C the self-consistent equations governing the coexistence of superconducting states φ k and ζ k . This is a necessary step for showing that the vertex-corrected electron-phonon interaction can support the spin-triplet odd-frequency state.

C. Equations for coexisting pairing channels
We repeat the derivation of Section II A in the most general way possible within the 2 × 2 Nambu space employed here. The electron self-energy given by Eq. (3) remains valid, while we modify the ansatz in Eq. (4) tô i.e., we now explicitly include theρ 2 channel. At this stage ζ k is just a mathematical definition, and we analyze its meaning later. With Eq. (21) at hand, the electron Green's function readŝ where we use the straight-forward definitions γ Due to the self-consistent inclusion of ζ k in our Eliashberg formalism, we get a total of four equations, reading

Instead of the three-component vectors used in Section II A, Eqs. (24)-(27) are defined in terms of
k are extended to a 4 × 4 pseudovector space, and read As a crosscheck, it is instructive to inspect Eqs. (24)-(27) with respect to their frequency symmetries. Starting with χ k and φ k , we know that the first summands of Eqs. (25) and (26), corresponding to the respective contributions due to the lowest order Feynman diagram, are even in ω k , since γ (χ) k and γ (φ) k are even in frequency. Further, in both equations the vertex corrections give rise only to products γ k1 · γ k2 · γ k3 that are even in ω k . This conclusion can be drawn by considering that γ are odd. The opposite holds when considering the functions ω k Z k and ζ k in Eqs. (24) and (27), respectively. Here, all first and second order contributions are odd along the Matsubara frequency axis, such that function Z k is even and ζ k is odd.
We note further that we can connect the current equations to the theory of Section II A by setting ζ k ≡ 0. As a consequence, also γ (ζ) k = 0 and we recover γ T k as 3-component vectors. Furthermore, for f = Z, χ, φ we need to consider only the upper 3 × 3 submatrix of Q In the following sections we apply the here-introduced theory to a cuprate model system. From here on, when referring to an induced state ζ ind k , we mean that the vertex-corrected Eliashberg equations of Section II A are solved self-consistently. Afterwards we calculate ζ ind k from Eq. (16) in a 'one-shot' calculation. On the other hand, when referring to function ζ k , we mean the numerical solution from the four coupled Eliashberg equations in Section II C. All calculations are performed with the Uppsala Superconductivity Code (UppSC) [8,[45][46][47][48][49], see also Appendix A for numerical details.

III. CUPRATE MODEL SYSTEM
We start the numerical analysis by considering a oneband tight-binding electron energy dispersion that is representative for the family of copper-based superconductors (cuprates). The electron energies are defined as with nearest and next-nearest neighbor hopping energies t (1) = 0.25 eV and t (2) = −0.1 eV. The chemical potential is fixed at µ = −0.09 eV. This parameter choice renders a well-nested Fermi surface (compare Fig. 2) characteristic for the cuprates, see e.g. Refs. [50,51] for a comprehensive overview. The Einstein phonon frequency and electron-phonon scattering strength are chosen here as Ω = 50 meV and g 0 = 150 meV, respectively.

A. Momentum dependence
Here we analyze the momentum structure of zerofrequency components as found from Eqs. (9)-(11) and Eq. (16) at T = 60 K. In Fig. 1(a)-(c) we show our results for Z = Z k,m=0 , χ = χ k,m=0 and φ = φ k,m=0 , respectively. As directly apparent from Fig. 1(b), the chemical potential renormalization has a nematic BZ structure, i.e., a mixture of s-wave and s ± -wave symmetry that breaks the C 4 to a C 2 rotational symmetry. The superconductivity order parameter φ has a d-wave shape with seemingly small nematicity [37]. After solving Eqs. (9)-(11) self-consistently we calculate the induced odd-frequency spin-triplet order parameter via Eq. (16). The result ζ ind = ζ ind k,m=0 is shown in Fig. 1(d). Similarly as φ [37], ζ ind obeys d-wave symmetry arising from the nesting properties of ξ k . Additionally, we learn that the induced function ζ ind is approximately one order of magnitude smaller than φ, and does not seem to have a nematic component. Note that the here-obtained (d-wave) odd-frequency order parameter is distinct from the oddfrequency s-wave, spin-triplet order parameter proposed originally by Berezinskii for 3 He [1].
For a closer investigation of the momentum symmetries we define seven different form factors, f (i) k ∈ {1, cos(k x ) + cos(k y ), cos(k x ) cos(k y ), cos 2 (k x ) + cos 2 (k y ), sin 2 (k x ) + sin 2 (k y ), cos(k x ) − cos(k y ), cos(k x ) cos(k y )[cos(k x )−cos(k y )]}, i = 1, . . . , 7. The form factors with indices i = 1, . . . , 5 belong to the A 1g representation, while i = 6 and i = 7 belong to the B 1g representation. Considering a function h, we can get the contribution in h due to BZ dependence f  In Fig. 1(e-h) we show the contributions due to each f (i) k for functions Z, χ, φ and ζ ind , respectively. Form factors belonging to the A 1g (B 1g ) representation are highlighted with purple (brown) background color. We have additionally tested other BZ symmetries, but the ones shown in Fig. 1 are clearly the most relevant.
From Fig. 1(e) we learn that the mass renormalization function has no finite projection in the B 1g channel, which is why Z is not nematic, compare also panel (a). Even though χ shows prevalent contributions from the A 1g representation, the projections for i = 6 and i = 7 are non-negligible, which is why the chemical potential clearly breaks C 4 rotational symmetry. A similar picture emerges for the nematic order parameter φ, see Fig. 1(g); here the dominant representation is B 1g , while small (but finite) contributions stem from A 1g . Turning to the induced term ζ ind , the projections in the A 1g channel are negligible to good approximation, which is why we observe a non-nematic d-wave state in Fig. 1(d).
So far we have shown that the induced order parameter is finite due to the vertex correction to Eliashberg theory, and results from a mixing of different momentum space symmetries. However, it remains to be shown that the actual (renormalized) interaction can support such an odd-frequency spin-triplet state, i.e., we need to solve the fully self-consistent Eqs. (24)- (27). We find that ζ k is in fact very similar to the 'one-shot' calculated ζ ind k . In Fig. 2(a) we show the induced function ζ ind as calculated before, projected on the renormalized Fermi surface, which is given by the condition On the other hand, using the fully self-consistent results for ζ = ζ k,m=0 does hardly lead to significant changes, as can be seen in Fig. 2(b). Both functions obey d-wave symmetry, and are of similar magnitude. Our calculations also reveal that the remaining functions Z, χ and φ are hardly affected, depending on whether the oddfrequency spin-triplet order parameter is included selfconsistently, or calculated in an isolated way.

B. Temperature evolution
In this section we discuss the temperature dependence of the superconducting state for our cuprate model system described by Eq. (32). For this purpose we define the experimentally observable superconducting gap ∆ k = φ k,m=0 /Z k,m=0 . This function describes an energy gap due to even-frequency spin-singlet Cooper pairs. Analogously, we define the odd-frequency spin-triplet gap function as η k = ζ k,m=0 /Z k,m=0 , or alternatively, η ind k = ζ ind k,m=0 /Z k,m=0 when we are concerned with the non-self-consistent induced order parameter ζ ind k,m=0 . We start by numerically solving Eqs. (9)-(11) as function of temperature, which leads to the blue circles in FIG. 3. Temperature dependence of the even-frequency spinsinglet (blue) and the odd-frequency spin-triplet gap function (times ten) (red). Open circles represent our computed results, while solid lines are obtained from fitting our data to Eqs. (35) and (36). (a) Results obtained from Eqs. (9)- (11) and (16). (b) Self-consistently computed temperature dependence from Eqs.  peratures at which ∆(T ) is finite we find the same BZ symmetries as in Section III A. The functional behavior of the even-frequency superconducting gap is accurately modeled by shown as solid blue curve. The variables a, b, and c in Eq. (35) are fitting parameters. The odd-frequency spintriplet order parameter η ind (T ) = max k |η ind k (T )|, as calculated from Eq. (16), is shown in red color in Fig. 3(a). Again, open circles represent our numerical results, while the solid line is found from a functional fit, here given as As direct comparison to η ind we also solved the fully self-consistent Eqs. (24)- (27) for η as function of T , leading to the curves shown in Fig. 3(b). It is directly apparent that results computed from the two different approaches only differ marginally. The zero-temperature limit of the even-frequency gap is in both cases fitted as ∆(T → 0) ≃ 10 meV, while the transition temperatures are given by T c ≃ 85 K. Further, we learn that the oddfrequency gap can equally accurately be determined in a single-shot or self-consistent calculation. This confirms the observations already made in Section III A, i.e., it is a good approximation to calculate ζ k (or equivalently η k ) via Eq. (16), requiring significantly less computational resources.
Our data reveals that the odd-frequency superconducting gap is approximately one order of magnitude smaller than its even-frequency counterpart. This finding is in agreement with previous adiabatic Eliashberg theory calculations of the odd-freuqency gap using ab initio input [8,9]. However, in contrast to these previous works where the odd-frequency order parameter gives rise to a paramagnetic Meissner effect, here the obtained odd-frequency superconductivity is expected to lead to a diamagnetic response. We suppose that this is the reason why it has been experimentally elusive so far, since both order parameters are expected to yield a diamagnetic Meissner effect. Therefore it is natural to associate the superconducting properties of the system with the dominant even-frequency order parameter only. Furthermore, it is apparent that η → 0 meV for small temperatures, which is consistent with earlier predictions [4,52], while the even-frequency gap approaches a finite value ∆(0) as T → 0 K. Note, that Fig. 3 contains few points for small temperatures, which stems from the necessary increase in the number of Matsubara frequencies to ensure convergence, becoming computationally unfeasible in the limit T → 0 K, see also Appendix A.

IV. SUMMARY AND DISCUSSION
We have shown that vertex corrections to the electronphonon problem in superconductors generally can lead to a coexistence of even-and odd-frequency superconducting states. Our analysis reveals that, due to a mixing of momentum space representations, the dominant spin-singlet even-frequency order parameter induces a spin-triplet odd-frequency superconducting gap, which, in comparison, falls short approximately one order of magnitude in amplitude. Both order parameters have the same critical temperature and are expected to yield a diamagnetic Meissner effect, hence, there is no straightforward way of verifying the induced odd-frequency spintriplet Cooper pairs experimentally.
Our numerical results confirm the existence of a finite odd-frequency d-wave order parameter in our cuprate model system, serving as a proof of concept. We cannot rule out that a similar picture can be found in other classes of superconductors, hence, it is possible that the induced state generically exists in a large variety of materials. Even though it has been argued in the past that odd-frequency superconductivity might only be realized for very anisotropic electron-phonon interaction mediated by acoustic phonons [4], or, alternatively, for spindependent couplings [5], our results are obtained for a completely isotropic electron-phonon coupling mediated by optical phonons. Our results are enabled by the renormalized vertex function, see Refs. [33,37,38], which is closely associated with Fermi surface nesting conditions. This allows for a self-consistent and strongly anisotropic interaction, supporting both even-and odd-frequency superconducting gaps to coexist. Our results furthermore imply that a stronger realization of an odd-frequency order parameter can be expected in nonadiabatic superconductors where vertex corrections beyond Migdal's approximation are important.
Lastly, it is currently an open question whether sys-tems with dominant odd-frequency order parameter exist. We speculate that, within the formalism used here, this is an unlikely scenario because the primary evenfrequency gap induces its odd-frequency counterpart. However, if we would assume odd-momentum space parity, the situation would be different. In this case, φ would still be dominant and describe spin-singlet odd-frequency Cooper pairs with e.g. p-wave BZ symmetry, while ζ would represent an induced spin-triplet even-frequency order parameter, a scenario that we leave for future investigations.