Full-bandwidth Eliashberg theory of superconductivity beyond Migdal's approximation

We solve the anisotropic, full-bandwidth and non-adiabatic Eliashberg equations for electron-phonon mediated superconductivity by fully including the first vertex correction in the electronic self-energy. The non-adiabatic equations are solved numerically here without further approximations, for a one-band model system. We compare the results to those that we obtain by adiabatic full-bandwidth, as well as Fermi-surface restricted Eliashberg-theory calculations. We find that non-adiabatic contributions to the superconducting gap can be positive, negative or negligible, depending on the dimensionality of the considered system, the degree of non-adiabaticity, and the coupling strength. We further examine non-adiabatic effects on the transition temperature and the electron-phonon coupling constant. Our treatment opens a pathway to systematically study vertex correction effects in systems such as high-$T_c$, flat band and low-carrier density superconductors.


I. INTRODUCTION
The foundation for establishing the microscopic description of phonon mediated superconductors was laid by the pioneering work of Migdal on the electron-phonon interaction in metals, where his famous theorem was introduced [1]. The essence of Migdal's theorem lies in the effective electron-phonon coupling λ times Ω/ F being small, where Ω is the characteristic phonon frequency and F the Fermi energy. Under this assumption it is possible to treat the arising infinite Feynman series of vertex diagrams in a perturbative manner with an expansion parameter λΩ/ F . Migdal's theorem dictates that when Ω F (or in other words, in the adiabatic limit), vertex corrections become negligible and the series can be truncated up to first-order in the coupling. For typical metals, where Ω/ F ∼ 10 −2 , such an approximation to the electron self-energy (below referred to as Migdal approximation), is commonly believed to be valid even in materials characterized by strong-coupling λ O(1). Eliashberg generalized this formalism to the superconducting state [2] and thereby laid the foundation for the thus-far most successful description of a vast amount of superconductors that deviate from the weak-coupling limit of the Bardeen-Cooper-Schrieffer (BCS) theory [3][4][5][6].
The importance of vertex corrections has been studied theoretically only by a few groups [25][26][27][28][29][30], while numerical analyses were so far carried out mainly by Pietronero and coworkers [31][32][33][34][35]. Recently, quantum Monte Carlo simulations have been employed to address the issue [36]. Due to the huge numerical complexity of including vertex corrections in the electronic self-energy, studies so far have been subject to further approximations, such as confining all considerations near the Fermi surface. Here we shall refer to such type of calculations as Fermi surface restricted (FSR). This includes the original works by Migdal and Eliashberg, too, and it is the reason why existing theories are often formulated in terms of the coupling λ, instead of the actual scattering strength g 0 arising in the Feynman expansion (see Section II A for details).
Our current work is motivated by the observation that many superconductors, including high-T c 's, depart from the perfect adiabatic regime, hence there is no consensus about the smallness of vertex corrections to the bare electron-phonon scattering diagram. Here, we address the question about the validity of the adiabatic approximation, using a less approximative treatment compared to available literature. It is worthwhile to establish under which conditions an adiabatic full bandwidth (AFB) [37][38][39] approach gives sufficiently accurate results, when even FSR Migdal-Eliashberg calculations [4] can be used for a successful description, and when the usage of nonadiabatic Eliashberg theory is inevitable.
We assume conventional (s-wave) Cooper pairing and consider the first two Feynman diagrams for the electronic self-energy, from which it is possible to derive a set of self-consistent equations for the mass enhancement, the chemical potential renormalization, and the superconducting order parameter (given below in Section II A). We numerically solve this system of equations without further approximation, i.e. the full momentum and frequency dependence is kept. In Section II B we introduce the AFB, as well as the FSR equations which we solve to have reference points for comparison. We use an effective one-band tight-binding model up to next-nearest neighbor hopping for the electronic energies. Our final model framework spans a huge parameter space, which allows systematic variations in the dimensionality, the characteristic phonon frequency, the strength of the coupling, the electronic bandwidth and, most importantly the degree of non-adiabaticity. In Section III A we study the influence of vertex corrections on the superconducting gap in systems of different dimensionality. We subsequently confine ourselves to 2D systems for which we examine similar parameter space with focus on the superconducting transition temperature and the electron-phonon coupling constant in Sec. III B. Our final conclusions and a brief outlook are given in Sec. IV.

II. THEORY
We want to describe electron-phonon mediated superconductivity with a non-adiabatic and full bandwidth Eliashberg theory that goes beyond the commonly employed Migdal approximation. As a starting point, we consider an electronic dispersion ξ b k , rigidly shifted by a chemical potential µ, such that ξ k = ξ b k − µ. We use k as Brillouin zone (BZ) wave vector and assume for simplicity an effective one-band situation. In Fig. 1 we show an example of a 1D dispersion with (blue) and without (brown curve) chemical potential shift. We define the minimal energy as F = |min k ξ k |, which is also indicated in Fig. 1. From this we explicitly see that the 'shallowness', i.e., the depth of the band can be controlled by µ alone. Turning to the lattice vibrations, we assume an isotropic Einstein phonon spectrum with a characteristic frequency Ω, see Sec. II A for details. By defining α = Ω/ F we can associate the electronic and bosonic energy scales with a parameter which is representative for the degree of non-adiabaticity in our system.

A. Microscopic description
In the following we use b q and b † q as bosonic annihilation and creation operator. The fermionic analogues c k,σ and c † k,σ , with spin index σ, are hidden in Nambu spinors Ψ † k = c † k,↑ , c −k,↓ in the electron-phonon coupled Hamiltonian that reads In the above we include a purely electronic part, a bosonic term and a coupling between both. Our focus here is solely on the effect of the latter term on superconductivity, therefore for the sake of simplicity we will not include the impact of Coulomb repulsion on the pairing. Here we use the Pauli basisρ i , i = 0, 1, 2, 3, and describe the ion displacement by u q , while g q is the electron-phonon scattering strength. The electronic propagator as a function of imaginary time τ readŝ where T τ is the time ordering operator.
Let ω m = πT (2m + 1) be a fermionic Matsubara frequency with temperature T and m ∈ Z. Then the Green's function in Matsubara space obeys the Dyson equationĜ with the shorthand notation F (k, iω m ) = F k,m for any function F . The non-interacting Green's function is given by [Ĝ 0 k,m ] −1 = iω mρ0 − ξ kρ3 . For evaluating the electronic self-energyΣ k,m we consider in Fig associated with a factor |g q | 2 ≡ |g 0 f q | 2 , which we rewrite as product of scattering strength g 0 and form factor f q (max |f q | = 1). The latter carries the information of the momentum dependence of the interaction. The expansion parameter in the Feynman diagram series is therefore |g 0 | 2 .
Taking into account both diagrams of Fig. 2 and identifying the difference iq m−m = iω m − iω m as imaginary bosonic frequency, we find the anisotropic electronic selfenergy aŝ Due to momentum conservation we use q = k − k in Eq. (4). As indicated before, the scattering amplitudes |g k−k | 2 and phonon propagators D k−k ,m−m are assumed to be branch and band independent. The Einstein-like phonon spectrum that we consider here leads to with V k−k ,m−m the electron-phonon interaction kernel. In Eq. (5) the electron density of states at the Fermi level is given by N 0 , and From here we follow the standard recipe of Eliashberg theory and introduce the mass enhancement Z k,m , the chemical potential renormalization χ k,m , and the superconducting order parameter φ k,m by writinĝ By combining Eqs. (3)(4)(5)(6) it is possible to project on channelsρ 0 ,ρ 3 andρ 1 , and find the self-consistent equations for corresponding prefactors in Eq. (6). After some tedious algebra we obtain For brevity we define the vector γ T while the matrices in Eqs. (7-9) have the form The superconducting gap can be found from ∆ k,m = φ k,m /Z k,m . The diagram in Fig. 2(a) corresponds to the left, while the diagram in Fig. 2(b) corresponds to the right summand in the brackets of Eqs. (7)(8)(9). Taking a step further from existing theories we did not perform any further energy integrations or other approximations in deriving Eqs. (7)(8)(9), hence our treatment is formally exact up to second order in |g 0 | 2 . Keeping all dependencies as they are stated above we solve for Z k,m , χ k,m and φ k,m in an iterative self-consistent loop; the results are presented in Sec. III below and we give further numerical details in Appendix A. The algorithm for solving the non-adiabatic, anisotropic and full bandwidth Eliashberg equations has been included in the Uppsala superconductivity (UppSC) code [37,38,[40][41][42][43].

B. Simplifications to the model
Before going to our numerical results it is useful to briefly discuss some limiting cases, with which we want to compare. From Eqs. (7)(8)(9) it is straightforward to obtain the adiabatic limit, in which we can neglect all contributions from the vertex in Fig. 2(b). This results in the well-known adiabatic Eliashberg equations where we use the label (ad) to avoid confusion with the non-adiabatic functions. Consequently, we obtain the gap function ∆ k,m . In cases where only processes at the Fermi level are relevant to describe the interacting system it is possible to derive a further simplified, but still anisotropic theory. Using label (Fs), for Fermi surface, the mass renormalization and gap function for these systems can be given as In deriving Eqs. (17) and (18) one assumes an infinite electronic bandwidth of a system at half-filling [4]. In Sec. III, we perform a parameter space exploration following each of these three different approaches. We have the Fermi surface based calculation of Eqs. (17) and (18), the AFB Eqs. (14)(15)(16), and the self-consistent full bandwidth description, that includes both scattering diagrams in Fig. 2, in Eqs. (7)(8)(9).
For the electronic energies we use a tight-binding description up to next-nearest neighbor, reading For simplicity we assume equal hopping energies along all three spatial directions and write t The c i are conveniently used to control the dimensionality (from 1D to 3D) of our system, hence c x = 1 and c i =x ∈ [0, 1]. Inserting in Eq. (19) gives i=x,y,z j=x,y,z;j =i from which we find the electronic bandwidth as W = max k ξ k − min k ξ k . By fixing W we are able to obtain the nearest neighbor hopping from We further take t (2) = t (1) /2 as the next-nearest neighbor hopping. For an initial choice of characteristic phonon frequency Ω and degree of non-adiabaticity α one can deduce the shallowness of our dispersion ξ k by F = Ω/α ! = |min k ξ k |, which uniquely determines the global chemical potential µ.
Another degree of freedom in our calculations is the coupling strength, which is commonly measured as λ, and treated here as parameter that is varied from weak to strong coupling. As stated in Sec. II A, the momentum dependence of the interaction kernel is g q = g 0 f q , so we can write In the case of a FSR calculation, as in Eqs. (17) and (18), the coupling can be extracted by a double momentum average at the Fermi level: Combining Eqs. (22) and (23) allows us to solve for the scattering strength which can then be used to calculate a kernel corresponding to λ for the adiabatic or non-adiabatic full-bandwidth equations.
To summarize, within the here presented setup we have the freedom to choose the systems dimensionality via c i , the characteristic phonon frequency Ω, the degree of nonadiabaticity α, the coupling strength λ and the momentum structure of the electron-phonon interaction. All subsequent results follow from the considerations in this section.

III. RESULTS
The full parameter space introduced in Sec. II B is too big for a complete numerical analysis, especially for 3D systems. We confine ourselves to a characteristic phonon frequency of Ω = 50 meV, fix the electronic bandwidth at W = 1.5 eV and impose an isotropic coupling, f q = 1 q . Due to the fact that existing theories of non-adiabatic superconductivity are confined to the Fermi surface [31][32][33] we perform a variation in λ to make a comparison easier. By means of Eq. (24) we have a tool to translate this variation into the expansion parameter |g 0 | 2 , which is more meaningful in our full-bandwidth treatment. For the here-presented results we are interested in reaching a qualitative understanding rather than quantitative absolute numbers, since we stay on a model basis.

A. Dimensionality
We want to compare our self-consistent non-adiabatic results with AFB, and with FSR calculations. For this purpose we consider the maximum superconducting gap as function of coupling strength for different degrees of non-adiabaticity α, repeated for each possible dimensionality. Further we fix the temperature at T = 20 K in the current section.

3D systems
It is well established that Migdal's theorem is valid in three spatial dimensions, hence we expect vertex corrections to be small, provided that α 1. Choosing α = 0.05 and c x = c y = c z = 1, we test this aspect in Fig. 3(a) by comparing the adiabatic maximum gap (green circles) and the FSR results (blue crosses) with outcomes from our non-adiabatic Eqs. (7-9) (red stars). As directly apparent, the vertex corrections in such an adiabatic situation are indeed very small. Additionally we find good agreement also with ∆ (Fs) , which points towards purely Fermi surface based Cooper pairing and indicates that a FSR approximation is sufficient in this parameter regime. The observations hold true for all couplings we tested, hence the essential ingredients are the 3D character of the system and smallness of α. We checked that there is good agreement between all three approaches for α < 0.05.
Next we examine the effect of making the system less adiabatic, α = 0.1. In similar color code as before we show the outcomes in panel 3(b). For increasing coupling strength the difference between FSR and AFB calculations grows larger. Further we find for all λ values that ∆ (ad) in this setup is an underestimation of the vertexcorrected maximum gap size ∆. We observe that these trends persist when further increasing α in Fig. 3(c). Here ∆ (Fs) is far too small compared to non-adiabatic results. Note, that this property stems from the fact that the FSR theory of Eqs. (17) and (18) does not explicitly depend on ξ k , and hence it is independent of α. Further, we find ∆ > ∆ (ad) throughout the whole range of couplings, while their relative difference stays comparatively small. In both panels (b) and (c) we get a deviation of |(∆ − ∆ (ad) )/∆| ∼ 20% in the large coupling limit. This can be explained in terms of the expansion parameter |g 0 | 2 , which is smaller for panel (c) than for panel (b). Hence there is an increased importance of non-adiabatic effects which, however, are weighted less. It is worth mentioning that Fig. 3(c) may be relevant to the situation encountered in H 3 S where α ≈ 1 and λ ≈ 2 [20,21] and in SrTiO3 3 where α 1 and λ ≈ 0.1 − 0.4 [23].

2D systems
Next we treat the two dimensional case by setting c x = c y = 1 and c z = 0. For this situation there is no proof that vertex corrections are negligible, even for α 1. It is further exceptionally interesting to examine possible effects in this setup, because many high-T c superconductors are either quasi-2D [7,44] or pure 2D systems [9,10,45,46]. We choose similar parameters as before, α = 0.05, α = 0.1 and α = 1, the results of which we plot in panels (a), (b), and (c) of Fig. 4, respectively. The color code and choice of axes (coupling versus maximum superconducting gap) is the same as in Fig. 3. A common observation in all three panels of Fig. 4   FIG. 4. Comparison of maximum superconducting gaps in 2D systems. Our results are self-consistently obtained from the three algorithms discussed in this work, using similar color code and degrees of non-adiabaticity α as in Fig. 3. is a rather small deviation between results from FSR and AFB calculations. This hints towards small tendencies of Cooper pairing away from the Fermi level in the parameter space we examine here [37,39]. Turning to the maximum superconducting gap from our non-adiabatic theory we find ∆ < ∆ (Fs) < ∆ (ad) for all couplings considered in Fig. 4(a). By neglecting vertex corrections one therefore overestimates the gap size to an extent that depends on λ. As we discuss in Sec. III B below, this overestimation is not restricted to the superconducting gap, but translates also to the transition temperature.
Before discussing the intermediate case of Fig. 4(b), let us first turn to α = 1 in panel (c). Here we find an opposite trend to before, i.e. an underestimation of the gap for all λ, and for both adiabatic algorithms. Again referring to Sec. III B this trend is similarly true for T c . An example superconductor where Fig. 4(c) may be relevant is FeSe/SrTiO 3 where α ≈ 1.6 − 2 and λ ≈ 0.2 − 0.4 [10,37,47]. Despite the fact that we here assumed a plain electron-phonon coupling, instead of the small-q interaction that is at play in FeSe/SrTiO 3 , our results provide further support that the latter mechanism is capable of mediating the observed high-T c . Moreover, our findings coincide qualitatively with non-adiabatic, but FSR calculations carried out on small-q scattering in C 60 compounds [12,13]. Let us now turn to the intermediate situation α = 0.1 in Fig. 4(b). For sufficiently large couplings we retrieve the situation of α = 1, while results for small λ resemble more closely the situation in panel (a). It can hence be argued that the transition from the adiabatic (α 1) to the non-adiabatic regime (α ∼ 1) is smooth, where in an intermediate situation as in panel (b), the over-or underestimation of ∆ (ad) and ∆ (Fs) with respect to vertex corrected results depends on the coupling strength only. It is interesting to note, lastly, that cuprates [7] and iron pnictides [44] fall in this intermediate regime.

1D systems
For completeness, we briefly discuss here results when lowering the system's dimensionality to 1D via c x = 1 and c y = c z = 0. In Fig. 5 we show again the comparison of results from Eqs. (7-9) with outcomes of Eqs. (14)(15)(16) and Eqs. (17)(18), just as is done in Figs. 3 and 4. For α = 0.1 (the case α = 0.05 has similar trends), we observe very small effects due to vertex corrections [see Fig. 5(a)]. Although a slight deviation from ∆ can be found for large λ, all three curves lie almost on top of each other. From this we can conclude that AFB, and even FSR calculations in 1D are very accurate with respect to the maximum superconducting gap, provided that α O(0.1). Further we note that ∆ > ∆ (ad) > ∆ (Fs) for all couplings shown. An increased influence of non-adiabatic effects is found with enhanced α = 1 in panel 5(b). When comparing ∆ and ∆ (ad) the situation is similar as in Fig. 4(b), i.e. the adiabatic results are an over-or an underestimation depending on the coupling strength. In addition we find a result unique to this 1D simulation: Among our three approaches, FSR calculations lead to the largest gaps for all coupling strengths we consider. From Fig. 5(a) and Fig. 3(a) one sees that in 1D, Migdal's approximation stays valid for larger values of the non-adiabaticity parameter as compared to the 3D case. This result is unexpected since in 1D Migdal's theorem should not hold [4]. It is worth mentioning that in 1D the nesting properties of the Fermi surface (that now consists of two points) enhance the tendency of the system towards the formation of a Charge Density Wave (CDW) or Peierls instability [48]. In principle, such tendency should be taken into account by including both superconducting and CDW orders on equal footing in the Eliashberg equations. However, doing this in the presence of vertex corrections is out of the scope of the present work.

B. Possible implications for 2D systems
From results in Sec. III A we learn that, depending on system specifics such as dimensionality, an increase or decrease of the gap size due to non-adiabatic vertex corrections is possible. So far we did, however, focus only on max ∆ and similar conclusions about the transition temperature cannot be drawn without further investigation. Strictly speaking, an enhancement in ∆ might either lead to a larger T c , or simply to a change in 2∆ 0 /k B T c . With ∆ 0 = lim T →0 max ∆ , this ratio is a common measure for how strongly coupled a superconductor is. Henceforth we closer examine 2D systems (c x = c y = 1, c z = 0) with degrees of non-adiabaticity α = 0.05, α = 0.1 and α = 1, all at a coupling strength λ = 1.5. Note that an exploration with respect to T c in a comparable parameter space as in Sec. III A is hardly feasible due to tremendous computational costs. Figure 6(a) shows the electronic energies along high-symmetry lines of the BZ with colors as indicated in the corresponding legend. The associated Fermi surfaces are drawn in panel (b) of the same figure.
As evident from these graphs, by changing the chemical potential µ we accomplish shallow Fermi surface pockets, which result in an increased α since Ω = 50 meV is kept constant. In Fig. 6(c) we draw a comparison of the temperature dependent maximum gap for different α and for each of our three approaches as indicated in the legend.
Let us start with the results from AFB calculations, shown as solid lines with similar color code as in panels (a) and (b). Our results for T c and ∆ (ad) (T ) do not change significantly as function of α. With growing non-adiabaticity we detect a small decrease in the transition temperature and the zero-temperature superconducting gap ∆ (ad) 0 ∆ (ad) (T = 10 K). When considering ∆ (Fs) , shown as dotted blue line, the gap size and transition temperature decrease further, but are still in the same range as the full bandwidth calculations. The here-detected difference is due to Cooper pairing away from the Fermi level [37,39], which is neglected when solving Eqs. (17)(18), and is almost negligible in the current model system. From the discussion of our adiabatic calculations one can conclude that T Inclusion of vertex corrections introduces some rather drastic changes compared to the adiabatic picture; the results of including these are shown as dashed lines in Fig. 6(c). For α = 0.05 (plotted in orange), both ∆ 0 ∼ 15 meV and T c ∼ 75 K are heavily reduced compared to the adiabatic counterparts, while 2∆ 0 /k B T c is enhanced to ∼ 4.97, giving rise to a seemingly more strongly coupled behavior. From Fig. 4 we can readily see that ∆ 0 > ∆ (ad) 0 > ∆ (Fs) 0 when α = 0.1. In addition we see from data plotted by the dashed green line that T c is enhanced to approximately 120 K, which gives again a stronger coupling 2∆ 0 /k B T c ∼ 4.84 when combined with ∆ 0 ∼ 25 meV. For even shallower bands, α = 1, we find the dashed purple line in Fig. 6(c) to give ∆ 0 ∼ 31 meV and T c ∼ 135 K, and hence a very strong coupling situation of 2∆ 0 /k B T c ∼ 5.33. If we consider changes in the transition temperature with the degree of non-adiabaticity for the vertex-corrected theory, our results for the current model system suggest that T c increases with α. Interestingly the trend for the maximum gap size in the limit T → 0 is reversed when comparing with the AFB outcomes, i.e. we get ∆ We similarly find the largest (smallest) T c (T (ad) c ) for the highest (lowest) α.
It is interesting to see how the just-discussed trends in critical temperature compare to effects in the effective electron-phonon coupling constant, which is given by This quantity is a measure for the coupling strength, renormalized due to interactions, and it is in general not equivalent to the bare coupling λ = 1.5, except for purely Fermi-surface based calculations. Note that a large number of Matsubara frequencies is needed to numerically confirm λ (Fs) = λ. In case of the non-adiabatic treatment we know from Eq. (7) that contributions due to both scattering diagrams enter the mass renormalization in an additive way. It is therefore convenient to define where λ (1) m and λ (2) m arise from the first and second diagram, respectively. For comparison it is useful to choose a temperature larger than T  Table I.
From the first row we learn that non-Fermi surface processes give rise to coupling constants that are clearly different from the initial λ = 1.5. The trend of how λ is modified as we increase the non-adiabaticity is not apparent. We note, however, that α = 1 leads to the  Fig. 6(c). Let us now consider solutions to the non-adiabatic Eqs. (7)(8)(9). As already mentioned, λ (1) m corresponds to the first order scattering diagram of Fig. 2(a), and hence can be compared to λ (ad) m . By comparing the first and second row in Table I we observe that indeed λ (ad) m λ (1) m . However, the deviation between these two quantities grows as α increases. This is due to the increasing importance of the feedback of the non-adiabatic terms on the adiabatic ones within the self-consistent iterative loop. The vertex corrections introduce the correction λ (2) , which according to the third row in Table I increases significantly with α. For the most adiabatic situation shown here, α = 0.05, we find λ (2) m < 0, hence the overall coupling constant is drastically reduced. In case of α = 0.1 we similarly get a negative λ (2) m , but the magnitude is very small. Therefore the sum λ m = λ (1) m + λ (2) m is closer to the bare coupling constant λ. Increasing the non-adiabaticity further to α = 1, the contributions due to the second order diagram become non-negligibly positive and give rise, together with λ (1) m , to a stronger electron-phonon coupling constant.
Via a closer inspection of the results for λ m we can give a qualitative explanation of trends for the critical temperatures as they are observed in Fig. 6(c). In the AFB calculations all couplings λ (ad) m lie within a narrow window of ±0.1 around the bare coupling λ = 1.5. It is therefore intuitive that corresponding values of T (ad) c do also not differ drastically from each other. When solving the non-adiabatic Eliashberg equations (7-9) we get an additional coupling contribution λ (2) m , which heavily depends on α. For the most adiabatic situation the effective final coupling λ m is comparatively weak, hence the critical temperature is smallest for α = 0.05. By increasing the degree of non-adiabaticity, we get pairing contributions from λ (2) m which make the overall coupling stronger. This, in turn, leads to a rather pronounced enhancement in T c , as shown in Fig. 6(c).

IV. CONCLUSION AND OUTLOOK
We have investigated the influence of vertex corrections to the electronic self-energy in electron-phonon mediated, anisotropic and full-bandwidth Eliashberg theory. To our knowledge the present investigation is the first numerical study to approach the challenge of non-adiabaticity in superconductors without involving further approximation, such as Fermi surface averages [12,31,[49][50][51] or momentum space clustering [30,52]. Our calculations numerically confirm the validity of Migdal's theorem for 3D systems. For this dimensionality we have found, within the explored parameter space, that the AFB calculations always resemble the vertex corrected results to a good degree regardless of α, which is however by no means true for FSR calculations. Contrarily, the observed trends in 2D systems suggest that non-adiabatic contributions are only negligible for rather small coupling strengths. We found that adiabatic results for the superconducting gap and T c can both be an over-or underestimation of the vertex corrected ones. In effectively one dimensional systems we observe all corrections due to non-adiabaticity to be rather small. Adiabatic FSR, as well as AFB calculations start to become less accurate for α O(1) and large couplings.
The model introduced in Sec. II B gives rise to a huge parameter space, which we have partially explored here. Further work in this direction is in progress. There is a large amount of extensions and refinements, as well as applications, that go beyond the scope of the present work. One example is to study anisotropy effects, i.e. cases when the bare pairing interaction is momentum dependent [37,38]. This includes electronic pairing mechanisms. For example, it has been shown that no analogue of Migdal's theorem exists for spin fluctuations, which means that vertex corrections are of similar order as the bare vertex [53]. Lastly, the here-presented methods can in principle be made compatible with ab initio input from Density Functional Theory (DFT) based calculations [40,43,54]. The current caveat in this respect is of numerical nature, since the computational complexity (see Appendix A) prohibits the usage of multiple electronic bands and very dense momentum grids. Nevertheless, such a non-adiabatic treatment would be an important step towards an even more realistic description of superconducting materials.