Ab initio theory of plasmonic superconductivity within the Eliashberg and density-functional formalisms

We extend the two leading methods for the \emph{ab initio} computational descrip tion of phonon-mediated superconductors, namely Eliashberg theory and density fu nctional theory for superconductors (SCDFT), to include plasmonic effects. Furth ermore, we introduce a hybrid formalism in which the Eliashberg approximation fo r the electron-phonon coupling is combined with the SCDFT treatment of the dynam ically screened Coulomb interaction. The methods have been tested on a set of we ll-known conventional superconductors by studying how the plasmon contribution a ffects the phononic mechanism in determining the critical temperature (\tc). Our simulations show that plasmonic SCDFT leads to a good agreement between predict ed and measured \tc's, whereas Eliashberg theory considerably overestimates the plasmon-mediated pairing and, therefore, \tc. The hybrid approach, on the other hand, gives results close to SCDFT and overall in excellent agreement with exper iments.


I. INTRODUCTION
Superconductors that do not fit into the standard BCS class, have opened interesting routes for alternative pairing mechanisms, whose applicability is still under debate [1][2][3] . Developing a first-principles method for the accurate calculation of the critical temperature (T C ) would not only clarify the microscopic mechanisms of superconductivity, but also contribute to the search for higher temperature superconductors. It has long been suggested that the key to high temperature superconductivity might be a purely electronic mechanism, that directly exploits the Coulomb repulsion between the electrons to provide their pairing. Many investigations have addressed the role of paramagnetic spin fluctuations in iron-based 4-6 and copper-oxide high temperature superconductors 2,7 . Other proposals, instead, have focused on the effective attraction 8 appearing in the dynamically screened Coulomb interaction due to the exchange of excitons [9][10][11][12] or plasmons [13][14][15] . In particular, the plasmon mechanism has been extensively investigated arguing that it could induce or significantly enhance superconductivity in many and very different classes of systems. These include perovskitic oxides [16][17][18] , metalchloronitrides 19,20 , organic superconductors 21,22 and light-element systems such as lithium metal and high pressure hydride superconductors [23][24][25][26] .
For conventional superconductors, calculations of T C are commonly based on Eliashberg theory [27][28][29][30] . This is, in principle, a comprehensive theory of the superconducting state, including both electron-phonon and Coulomb effects. The usual application of Eliashberg theory to realistic systems is, however, oversimplified 31 in that the Coulomb interaction is assumed not to favor Cooperpair formation, and is reduced to a single parameter µ * 30-32 . The standard Eliashberg framework, thus, is not suitable for a quantitative description of superconductivity supported by electronic mechanisms. Unlike Eliashberg theory, the extension of density functional theory to superconductors 33 (SCDFT) does not involve any semi-empirical approximation for the Coulomb interaction, and enables calculations of T C entirely from first-principles. Nevertheless, SCDFT was formulated to address conventional superconductivity 34,35 , so that it employs a static screening of the Coulomb repulsion 36 . Recently, a generalization of SCDFT for applications to plasmonic superconductivity has been proposed 23,37 . However, in this theory plasmonic effects are included in the superconducting state, but neglected in the normal state.
In this work we extend Eliashberg theory (Secs. II) and SCDFT (Sec. III) to provide ab initio calculations of plasmonic effects on the superconducting properties of real materials. In both frameworks retardation effects in the phonon-mediated and screened Coulomb interactions are treated on the same footing by keeping their characteristic frequency dependence. By applying these methods in Sec. V, we study how the plasmon contribution affects phonon-induced superconductivity for a set of materials representing the main families of conventional superconductors. and anomalous (off-diagonal, F ) components describing, respectively, single-particle electronic excitations and Cooper pairs. The matrix Green's function is determined via the Dyson's equation: whereḠ 0 is the normal-state Green's function of the noninteracting electron system andΣ(k, iω n ) =Σ c (k, iω n ) + Σ ph (k, iω n ) is the electron self-energy associated with the screened Coulomb and phonon-mediated interactions. G 0 can be constructed from the Kohn-Sham (KS) states |k ≡ |kl and eigenvalues ε k of density functional theory (DFT) in the usual form where τ {0,...,3} are the Pauli matrices and the energy ε k is measured relative to the chemical potential. The pairing mechanism is dominated by the phononmediated interaction that, being retarded, overcomes the (almost instantaneous) Coulomb repulsion between the electrons. Since the phonon energy scale, set by the Debye frequency ω D , is much smaller than the electronic Fermi energy E F , the method relies on Migdal's theorem to treat the electron-phonon interaction accurately to order ω D /E F . The key approximation consists in retaining forΣ ph (k, iω n ) only the diagram for the dressed phonon exchange by the self-consistently dressed electron propagator (ḠW ). Due to the absence of an analogous theorem, the treatment of the Coulomb interaction is much harder. On the other hand, the possibility of a Coulomb enhancement of the T C is neglected. Coulomb effects are largely accounted for by normal state parameters, i.e., the electron and phonon quasiparticle energies, ε k and ω qν , and the screened electron-phonon coupling g kk ν . In addition, there remains a static screened Coulomb repulsion W (k, k ), which counteracts superconductivity. Within the Eliashberg approximation, the phonon and Coulomb contributions to the electron self-energy read, respectively, as and (4) Following the standard practice 30 , the anisotropic electron-phonon coupling λ k,k (iν n ) in Eq. (3) is defined by the spectral representation in terms of the Eliashberg function where N (0) is the electronic density of states at the Fermi level. Eq. (4) includes the subtraction of the exchangecorrelation potential of KS-DFT, v xc , so that the resulting Coulomb self-energy is purely off-diagonal. This prevents one from double counting Coulomb effects in the normal state, which are already included in the KS band structure ε k enteringḠ 0 . The total self-energy is more conveniently rewritten in terms of three scalar functions given by the coefficients of the Pauli matrix representation forΣ: These are the mass renormalization function Z(k, iω n ), the energy shift χ(k, iω n ), and the order parameter φ(k, iω n ). Through the Dyson's equation (1), the calculation ofḠ is reduced to solving three coupled equations for Z, χ and φ. In particular, the function ∆(k, iω n ) = φ(k, iω n )/Z(k, iω n ) plays the role of the superconducting energy gap, whereas the quantity χ(k, iω n ) leads to a shift of the chemical potential, which little affects the formation of the superconducting state. By neglecting χ, the equations of interest for the τ 0 and τ 1 components of the Eliashberg self-energy take the form where Note that, since retardation effects in the Coulomb repulsion are disregarded, Z(k, iω n ) is entirely determined by the phonon-mediated interaction 30 . Moreover, the Coulomb contribution to φ(k, iω n ), given by the second term of Eq. (9), is frequency independent. Several approximations are commonly employed in order to reduce the workload involved in solving the Eliashberg equations (8) and (9). Essentially, since the superconducting pairing occurs mainly within an energy window ∼ ω D around the Fermi surface, the equations are simplified by averaging over k and k in the electronic states on the Fermi surface as follows: Although being quite accurate for phonons, this approximation is not justified for the Coulomb interaction, which may remain large for energies up to E F . In practice, following the arguments of Morel and Anderson 40 , W (k, k )  Figure 1. Left: Eliashberg superconducting gap function φ for bulk Nb in the full dynamical (blue), weak-coupling dynamical (red) and static (black) Coulomb approach. Right: Eliashberg mass renormalization function Z and its decomposition (for the dynamical case) into Coulomb Z c and phononic Z ph components. All the quantities are computed in the low temperature limit.
in Eq. (9) can be replaced by a strongly-reduced pseudopotential µ * , with an energy cutoff ω c ∼ 10ω D , which effectively accounts for the Coulomb scattering of electrons far from the Fermi surface. The Morel-Anderson pseudo-potential is defined by the expression where µ ≡ N (0) W (k, k ) F S . In most applications, however, µ * is treated as a semi-empirical parameter fitted as to reproduce the experimental critical temperature. With the above mentioned approximations, the Eliashberg approach involves solving numerically the following isotropic equations: A. Plasmonic extension of Eliashberg Theory We go beyond Eq. (4) for the Coulomb self-energy by assuming theḠW approximation, i.e., we still neglect vertex corrections, but introduce the dynamical screening of the Coulomb interaction through the frequencydependent dielectric function. Hence, we consider the self-energy expression where the Coulomb potential W k,k (iν n ) is obtained from the symmetrized dielectric function GG (q, iν n ) as Here, Ω is the unit cell volume, G are reciprocal lattice vectors and q is the difference k − k reduced to the first Brillouin zone. The symmetric form of the dielectric function is defined by For numerical convenience, we rewrite W k,k (iν n ) in the form are the matrix elements of the bare Coulomb interaction and is the spectral function of the electronic polarization. Note that the second term on the right hand side of Eq. (18) is formally equivalent to the spectral representation of the anisotropic electron-phonon coupling (Eq. (5)). The screened Coulomb potential in its spectral representation can be separated into a static and a dynamical part, where the latter, given by incorporates plasma oscillations. We approximate Eq. (20) by its average taken over the corresponding surfaces of constant energy, ε, in k-space as It should be observed that Eq. (22) is a generalization of Eq. (11) for the conventional isotropic Eliashberg theory. By using Eqs. (22) and (23) for the Coulomb interaction in the expression for the self-energy, we obtain the following Coulomb contributions to the Eliashberg functions Z and φ: The influence of these terms on T C can be easily seen by considering that in the simple BCS limit one has , where 1+λ comes from the electron-phonon Z term in the Eliashberg equations.
Here, the dynamical contribution to the anomalous kernel φ c (given by the second term of Eq. (25)) enhances the Coulomb repulsion µ between the electrons in the energy scale of the plasmon frequency ω pl . On the other hand, since ω D ω pl , the effective Coulomb repulsion decreases from the original value µ * , favoring superconductivity (higher T C ). This effect, however, is counteracted by the Coulomb correction Z c to the effective mass, which adds up to the phononic term (1 + λ), contributing to the reduction of T C . Fig. 1 shows the Matsubara frequency dependence of the mass renormalization Z and gap function φ for bulk Nb. The inclusion of retardation effects in the Coulomb interaction leads in φ to large negative tails at high energy. Since the high-energy gap function is negative, the plasmonic coupling serves as an effective attraction, which, according to Eq. (25), increases the value of φ at the Fermi level. This effect is less pronounced when Z c is included. The right panel of Fig. 1 shows the decomposition of Z into phononic and Coulomb contributions. Z ph has a peak at low frequency with energy width of the order of the Debye energy and above ω D converges to 1. Z c , instead, which is non-zero only in the dynamical approach, decays at the plasmonic energy scale.
As evident from Fig. 1, the main difficulty in solving the Eliashberg Eqs. (24) and (25) is related to the fact that the integration (both in ω n and ε) has to be performed on a huge energy scale, and therefore cannot be tackled by brute force computation. Just to give an indication, at T = 2.8 K, the number of Matsubara points within the plotted energy range is of the order of 70 thousand, and reaching a tight convergence would require an even larger energy window of several hundred eV. To overcome this slow convergence problem in the numerical implementation of the equations we have adopted the following strategies: i) We have used a logarithmic ε integration mesh. This allows for a dense discretization at low energy, where variations in the func-tions have to be accounted for more accurately, but extends up to arbitrarily large energies with relatively few additional points. ii) Similarly, we have adopted a nonhomogeneous mesh of Matsubara points. Since Matsubara frequencies are fixed by the temperature, a non-linear mesh can be obtained by pruning points and redistributing their weight. The resulting Matsubara mesh at 2.8 K is indicated by orange ticks in the left panel of Fig. 1. iii) The dynamical Coulomb interaction itself depends on the (bosonic) Matsubara frequencies and on the energy. When computing the interaction from first principles, a huge computational cost is associated with the calculation of the matrix elements of the dielectric function at high frequency, with respect to the KS states at high energy. For this reason, we have introduced high-energy cutoffs in ω n and ε (typically of the order of 50-100 eV). Above this energy, the dielectric function of the material is replaced with that of the homogeneous electron gas in the plasmon-pole approximation. The parameters which enter the latter are fitted to the explicitly computed values of S(ε, ε , ω) at the cutoff, so to ensure a good overall match to the actual interaction. This approach not only reduces the numerical cost in computing the interaction, but also allows for the analytical integration of the Matsubara frequencies from the cutoff energy to infinity. We point out that these techniques do not introduce additional errors in the method. However, they involve convergence parameters that have to be carefully chosen in order to achieve the correct numerical result.

III. DENSITY FUNCTIONAL THEORY FOR SUPERCONDUCTORS
Density functional theory for superconductors (SCDFT) is an extension of conventional DFT for ab initio calculations of material-specific properties in the superconducting state 33 . The theory includes the superconducting order parameter χ sc (r, r ) as an additional density. The corresponding non-interacting KS system then reproduces, in principle exactly, both the normal density and the superconducting order parameter of the real system. In the so-called decoupling approximation (on which Eliashberg theory is also based), the KS system is fully determined by solving the BCS-like gap equation where E k = ε 2 k + |∆ s k | 2 are the KS excitation energies and β is the inverse temperature. The kernel of the equation consists of a diagonal part, Z k = Z ph k , and a nondiagonal part, K k,k . Z ph k plays the role of the renormalization function in the Eliashberg equations, whereas K k,k = K c k,k + K ph k,k , which includes both Coulomb and phonon-mediated effects, is responsible for the binding of the electrons in Cooper pairs. Compared to Eliashberg theory, SCDFT features two major advantages: (i) the treatment of the Coulomb repulsion does not resort on any empirical parameter µ * , (ii) all the Matsubara frequency summations are evaluated analytically in the construction of the exchange-correlation (xc) kernels. As in Eliashberg theory, phonon dynamics is properly included, but at the same time the gap equation retains the form of a static BCS equation. Hence, Eq. (26) allows one to account for the full anisotropy of materials at a low computational cost. However, the accuracy of the method is bound by the quality of the available functionals.
Making a connection to many-body perturbation theory, approximate xc kernels have been derived from approximations for the xc self-energy operator, via the Sham-Schlüter equation in Nambu space. The first SCDFT functional by Lüders, Marques and co-workers (LM) 34,35 employed theḠ s W approximation for the self-energy in the statically screened Coulomb repulsion and phonon-mediated interaction. By construction, this functional neglected higher order processes included in Eliashberg theory by the self-consistent dressing of the electron Green's function in theḠW self-energy. The LM approximation, thus, was not validated by Migdal's theorem, which made it of questionable accuracy for treating electron-phonon coupling effects. To solve this issue, Sanna, Pellegrini and Gross (SPG) 41 have recently introduced a parametrization of the functional based on the electron-phonon Eliashberg self-energy for a simplified (Einstein) phonon spectrum. The new SPG kernels, give superconducting transition temperatures and gaps in excellent agreement with experiments 42-47 , while still having a simple analytic form.
Further extensions and applications of the LM functional have addressed the description of superconductivity in the presence of magnetic fields 48,49 and in real space 50 , the inclusion of spin-fluctuations contributions to the pairing 51,52 and the treatment of the dynamical screening of the Coulomb interaction 20,[22][23][24]37 . In a first attempt to introduce plasmonic effects in SCDFT, Akashi and Arita 23,37 have proposed a dynamical correction to the pairing kernel K c k,k by retaining the frequency dependence of the Coulomb interaction at the RPA level in the exchange anomalous self-energy. The method, implemented in the multipole plasmon approximation, has given a systematic increase of the T C by 10−20% in compressed sulfur hydrates H 2 S and H 3 S, and by a factor of 2 in Al and Li under pressure.

A. Plasmonic SCDFT with mass term
The results of plasmonic Eliashberg theory for Nb ( Fig. 1) suggest that the Z c term, stemming from the diagonal part of the dynamical Coulomb self-energy, should play a major role in determining T C . Z c can be viewed as the Coulomb counterpart of the mass renormaliza-  Figure 2. Left: SCDFT gap function ∆s for bulk Nb in the dynamical (blue), weak-coupling dynamical (red) and static (black) Coulomb approach. Right: SCDFT Z kernel and its decomposition (for the dynamical case) into Coulomb Z c and phononic Z ph components. All the quantities are computed in the low temperature limit and using the SPG phononic functional.
tion enhancement, 1 + λ, that corrects the BCS predictions to their strong coupling values in Eliashberg theory 30,53 . Accordingly, it is expected to be relevant for strong electron-plasmon interactions. As discussed above, the recently developed SCDFT scheme for plasmonic superconductivity neglects this contribution by assuming a null diagonal kernel Z c and can, thus, be regarded as a weak-coupling plasmonic theory. In this section we propose a more general SCDFT approach, which also includes plasmonic corrections to the mass enhancement. By following a procedure analogous to that presented in Ref. 54 for the treatment of the electron-phonon coupling, we construct the SCDFT plasmonic kernels from theḠ s W self-energy in the screened Coulomb potential. In the isotropic approximation (Eqs. (22) and (23)), we obtain the following expressions: where the quantity I(ε, ε , ω) is defined in terms of the Fermi and Bose distribution functions f and b by In Fig. 2 we show the energy dependence of the calculated KS gap function and kernel Z for bulk Nb. As in Eliashberg theory, plasmonic contributions enhance the high-energy negative gap. The effect is much more pronounced in the weak-coupling approach, within which the value of the KS gap at the Fermi level almost doubles compared to the static and full dynamical cases. The kernels Z c and Z ph on the right panel of Fig. 2 have the shape of two over-imposed peaks. The sharper peak is strongly temperature dependent, and occurs very close to the Fermi level, at |ε| < 10 −4 eV . This energy scale is not directly related to that of the couplings, but arises from the constraint that the KS system should reproduce the interacting normal and anomalous densities 41 . On the other hand, the broader peak decays on the energy scale of the interactions, i.e., at the plasmon energy in Z c and at the Debye energy in Z ph .

IV. HYBRID ELIASHBERG
By virtue of Migdal's theorem 27,30 , the F W approximation for the anomalous self-energy in Eliashberg theory describes the phonon-mediated pairing very accurately. On the other hand, there is no a priori indication that the F W scheme improves over F s W for the treatment of plasmonic effects. Here, we consider a hybrid Eliashberg-SCDFT theory in which the Coulomb part of the pairing self-energy is in the F s W form, where F s is the KS Green's function which reproduces the superconducting order parameter in the Eliashberg approximation.
Since the KS system has the same anomalous density of the interacting system, the following equality holds 34,35,41 : where F s (ε, iω n ) = ∆ s (ε) ω 2 n + E 2 with E = ∆ 2 s + ε 2 . Using the Eliashberg Green's function F (ε, iω n ) = φ (ε, iω n ) /Θ (ε, iω n ) to compute χ sc and evaluating the Matsubara frequency summation on the right hand side of Eq. (31), yields: The obtained F s is then used as an input for the Eliashberg equation which determines the Coulomb gap func-tion, i.e., Eq. (25) is replaced with  Table I. Electron-phonon (λ) and electron-plasmon (Z c 0 ) coupling strengths for the test set of materials with associated experimental critical temperatures T exp C .

A. Material set
In order to assess the accuracy of the methods discussed above, we have investigated the effect of the electron-plasmon coupling on the transition temperatures of a set of conventional superconductors. Our set includes experimentally well-characterized systems chosen to cover a wide range of properties, i.e., elemental (Al, Sn, Ta, Nb and Pb) and binary phonon-mediated superconductors (TaC, ZrN, V 3 Si and CaC 6 ), ranging from weak to strong coupling. To keep the entire procedure ab initio, all the calculations have been performed at the theoretical lattice parameters obtained by means of the PBE functional 55 . The dielectric function entering the dynamical Coulomb kernel has been calculated within the random phase approximation (RPA), using the fullpotential LAPW code Elk 56 . Fig. 3 shows the electron  542  149  352  35  355  227  88  355  779   112  67  145  13  152  120  29  115  264   25  49  36  4.9  45  32  energy loss spectra in the low-q limit for the chosen set of materials. For both Al and Sn, which are free electronlike metals, one observes, similarly to the homogeneous electron gas, a single pronounced plasmon peak, centered respectively at 16 and 14 eV. All the other systems show a more complex spectrum, varying from the two-peak structure of V 3 Si, to the broad distributed structures of TaC. Apart from the low-energy plasmons of CaC 6 and ZrN, the main plasmonic structures are located at energies above 10 eV. The electron-phonon and electron-plasmon coupling strengths for the chosen materials are summarized in Tab. I, together with the experimental T C 's. The electron-phonon coupling is expressed in terms of the BCS-like coupling constant λ defined as the static limit of λ (iν n ) in Sec. II. The electron-plasmon coupling, being strongly energy dependent, cannot be reduced to a simple isotropic parameter, and is thus represented by the energy integrated quantity Z c 0 , defined as Z c (ε, iω n ) of Eq. (24) computed at the Fermi level (ε = 0) and ω n = 0.

B. Eliashberg
In Fig. 4 and Tab. II (columns A-C, L) the experimental values of T C for the chosen set of materials are compared to the values calculated within the Eliashberg approach from Eqs. (24) and (25) by employing the static and dynamical screening of the Coulomb interaction. For the dynamical case, the weak-coupling results obtained by neglecting the plasmon-induced mass renormalization term Z c are also shown.
It is evident that the static approximation gives better values for T C , whereas the plasmonic theory systemati-cally overestimates the experimental data. The inclusion of Coulomb retardation effects in theḠW approximation yields predicted temperatures that are on average two times bigger than the corresponding experimental values. Notably, the discrepancy between theory and experiments becomes huge when the term Z c is neglected. Since the inclusion of dynamical screening effects in the Coulomb interaction brings the theory a step closer to being exact, one would expect an improvement in the calculated values of the critical temperature. The apparent worsening of the results can be traced back to the neglect of vertex corrections in the Coulomb selfenergy diagram 57,58 and/or the breakdown of the RPA for W. Regarding this point, we should mention that going beyond RPA by using linear-response time-dependent DFT 59,60 within the adiabatic local density approximation in the calculation of the dielectric function does not improve significantly the quality of the results. These aspects will require further investigations. As a matter of fact, our ab initio treatment of the static Coulomb interaction in Eliashberg theory appears to be very accurate, confirming previous results along these lines 25,61,62 .

C. SCDFT
In this section we present the results obtained within the SCDFT framework. As for Eliashberg theory, we consider static, plasmonic weak-coupling (Z c = 0, Eq. (28) for K c ) and strong-coupling (Eqs. (27), (28) By using the phononic LM functional, the SCDFT results obtained within the static approximation for the screened Coulomb interaction underestimate the experimental values by an average error of 35%. The inclusion of plasmonic effects in both the normal and the superconducting state yields even lower T C 's, with an average error of 43%. On the other hand, if only the plasmonic contribution to the superconducting pairing is accounted for (i.e., Z c = 0), the theoretical results systematically overestimate the experimental data by a factor of 2. In spite of these deviations, plasmonic SCDFT with the phononic LM functional gives closer T C 's to experiments compared to Eliashberg theory.
As already mentioned, the SPG functional improves over the LM approximation and is comparable in accuracy to conventional Eliashberg theory in describing electron-phonon effects 41 . Employing this functional together with the static Coulomb kernel gives results in good agreement with the experiments. The agreement worsens considerably by including plasmonic effects in the weak-coupling approximation, as this leads to a sizable increase of the predicted T C 's. Nevertheless, adding the plasmonic renormalization mass factor Z c suppresses  the T C 's values and increases the overall accuracy. The average percentage error in this latter dynamical approach is less than 20%. For the chosen set of materials, this approximation turns out to be the most accurate, as reported in Tab II. However, it should be noticed that all the theoretical results have a non negligible intrinsic error due to the approximations made in calculating the phonon spectral function. For this reason it is not possible to precisely rank in accuracy the different methods. Nevertheless, we can say that plasmonic effects can be safely incorporated in the SCDFT scheme, as they introduce a relatively weak correction to the phonon-induced T C , which appears to be consistent with the experimental results.
In Sec. V B we have mentioned that the failure of plasmonic Eliashberg theory could be ascribed to the RPA screening or the neglect of Coulomb vertex corrections. The higher accuracy of plasmonic SCDFT, which employs the same Coulomb propagator, indicates that the RPA is not the main source of error. On the other hand, plasmonic SCDFT relies on theḠ s W approximation for the Coulomb self-energy, whereas Eliashberg theory amounts to the fully self-consistentḠW . This leads us to speculate that vertex corrections to the Coulomb selfenergy might be mostly cancelled by the self-consistent dressing of the KS electron Green's function inḠ s W .

D. Hybrid Eliashberg
From the results of plasmonic Eliashberg theory is evident that the F W approximation for the anomalous Coulomb self-energy significantly overestimates the T C . Since Eliashberg theory is a routinely used method for the prediction of superconducting properties, this appears as a major drawback. A viable alternative is the hybrid Eliashberg-SCDFT approach proposed in Sec. IV. This, in fact, employs the F s W approximation, which better describes the plasmonic contribution to the superconducting pairing.
The T C 's calculated for our test set of materials are collected in Tab. II (columns from J to K) and compared to the experimental data in Fig. 6. Consistently with all the previous weak-coupling calculations, the results without the plasmonic mass term tend to overestimate the T C . In this case the overestimation is, on average, by about 60%, considerably improving over Eliashberg theory. On the other hand, the fully dynamical approach leads to predicted temperatures that are very close to the experiments. | is the deviation from the experimental TC. For each method and approximation we indicate the average percentage error av. %|err|=(100 |err|/T exp C ), the average error (av. |err|), and the maximum error (max |err|) over the material set.

VI. CONCLUSIONS
We have presented an extension of Eliashberg theory and SCDFT to include the dynamical screening of the Coulomb interaction. Our analysis points at the importance of the plasmonic mass terms, which largely counterbalance the effect of the plasmon-mediated attraction in the Cooper pair. The computational cost associated with the inclusion of the frequency-dependent Coulomb interaction is made affordable by employing an energyresolved isotropic approximation and by setting nonlinear energy and frequency integration meshes. A hybrid Eliashberg-SCDFT scheme is also formulated, which combines the ME approximation for the electron-phonon coupling with the SCDFT treatment of the dynamically screened Coulomb interaction. The accuracy of the approximations employed in the different methods has been assessed by calculating the plasmon contribution to the critical temperature for a set of classic superconductors. Our simulations show that the SCDFT plasmonic kernels, combined with the phononic SPG functional, yield a good agreement between predicted and measured critical temperatures. Dynamical corrections turn out to be small but not negligible, being of the order of 10-15% of T C . Eliashberg theory, although accurate in the static limit of the screened Coulomb interaction, when plasmonic effects are included leads to a large overestimation of T C , by an average factor of 2. Dynamical Coulomb effects can be included in Eliashberg by adopting the hybrid approach, which gives results close to SCDFT and overall in excellent agreement with experiments.