Eliashberg theory with ab-initio Coulomb interactions: a minimal numerical scheme applied to layered superconductors

We present a minimal approach to include static Coulomb interactions in Eliashberg theory of superconductivity from first principles. The method can be easily implemented in any existing Eliashberg code (isotropic or anisotropic) to avoid the standard use of the semiempirical parameter μ∗ , which adds unnecessary uncertainty to T c predictions. We evaluate the prediction accuracy of the method by simulating the superconducting properties of a set of layered superconductors, which feature unconventional Coulomb effects: CaC6, MgB2, Li-doped β-ZrNCl and YNi2B2C. We find that the estimated critical temperatures are consistent with those from ab-initio density functional theory for superconductors, and in close agreement with the experimental values.


Introduction
Calculations of the critical temperature (T c ) of conventional superconductors are commonly based on Eliashberg theory [1][2][3][4]. 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 [5] in that the Coulomb interaction is reduced to a single semi-empirical parameter µ * [4][5][6]. In recent years the problem of developing a fully ab-initio Eliashberg theory accounting for Coulomb and phonon interactions on equal footing has been addressed in detail [7][8][9][10]. An extension of the Eliashberg approach to include dynamical Coulomb interactions is computationally very expensive as it involves the evaluation of slowly converging Matsubara frequency summations. Smart summation algorithms for the simulation of real materials have been developed, such as those presented in [9] and [10]. Nevertheless, the computational cost is still very high compared to density functional theory for superconductors (SCDFT) [10,11] and conventional Eliashberg codes in the µ * approximation [12].
On the other hand, for most phonon-driven superconductors plasmonic pairing contributions are negligible, and accurate predictions can be obtained by a proper treatment of the static Coulomb interaction [8,10]. We show that in this case the computationally demanding parts of the Matsubara sums can be evaluated analytically, similarly to how it is done in SCDFT [10,11,[13][14][15][16]. The approach that we present follows the lines of a scheme discussed by Scalapino, Schrieffer and Wilkins [5] and corresponds to the static limit of the hybrid Eliashberg-SCDFT approach by Davydov et al [10]. For classic elemental superconductors we do not expect this method to lead to significant differences with respect to the µ * = 0.1 rule, which was meant to describe these materials [4,17,18]. Deviations are foreseen for less common superconductors where Coulomb interactions have a non-trivial behavior, like Chevrel phases [19], multigap superconductors [20,21] and low dimensional systems. In this work we focus on the application of the method to a set of layered compounds, namely, CaC 6 , MgB 2 , Li 0.5 ZrNCl and YNi 2 B 2 C.
The paper is organized as follows: in section 2 we review the anisotropic Eliashberg theory and discuss our approach for an analytic Matsubara integration of high-energy Coulomb effects (section 2.1). Further simplifications are introduced by employing suitable isotropic approximations (section 2.2). The isotropic approach is used to estimate the superconducting properties of the chosen set of materials. Results are presented and discussed in section 3.

Eliashberg theory with static Coulomb interactions
Eliashberg theory is a many-body Green's-function approach for the description of conventional superconductors, where the pairing is driven by the phonon-induced attraction between the electrons [2-5, 12, 22]. For a recent review of conventional Eliashberg theory and its current extensions we refer the reader to [10]. Here, we use the same notation to summarize the key aspects of the method. Eliashberg theory relies on Migdal's theorem to treat the electron-phonon interaction accurately to order ω D /E F , where the phonon energy scale, set by the Debye frequency ω D , is assumed to be much smaller than the electronic Fermi energy E F . This amounts to treat the electron self-energyΣ(k, iω n ) up to second order in the electron-phonon coupling. On the contrary, there is no small parameter that allows for a satisfactory perturbative treatment of the Coulomb interaction. The possibility of a plasmon enhancement of T c is however neglected. The electronic band structure ε k is assumed to be weakly perturbed by the superconducting condensation and accurately described by normal-state density functional theory [23][24][25]. Additionally, one has to take into account a static Coulomb interaction W k,k ′ , which opposes the formation of the superconducting state. Under these assumptions, the Eliashberg self-energy is usually expressed in terms of the mass renormalization function Z(k, iω n ) and the superconducting order parameter ϕ(k, iω n ), through its decomposition into Pauli matrices (τ 0,1 ) [10]: where k = (k, n) is a combined momentum and band index and ω n are fermionic Matsubara frequencies.
The problem of calculatingΣ (and ultimately the one-electron Green's function) is then reduced to solving coupled equations for Z and ϕ. These equations are: where we have introduced the function: and split ϕ(k, iω n ) into phonon (ϕ ph ) plus Coulomb (ϕ c ) contributions. Note that, since retardation effects in the Coulomb repulsion are disregarded, ϕ c (k) is frequency independent and Z(k, iω n ) is entirely determined by the phonon-mediated interaction λ k,k ′ (iν n ) [4]. The electron-phonon coupling is usually defined by its spectral representation, in terms of the Eliashberg function: where N(0) is the electronic density of states at the Fermi level and g kk ′ ν are the electron-phonon matrix elements for the scattering between electronic states k and k ′ through a phonon with wave vector q = k − k ′ and branch index ν. W k,k ′ is given by the matrix elements of the screened Coulomb interaction with respect to the Kohn-Sham orbitals, where G are reciprocal lattice vectors and the static approximation for the dielectric matrix ϵ −1 GG ′ (q) is made.

Analytic summation of the Coulomb interaction
Solving the Eliashberg equations becomes computationally affordable if a small cutoff compared to E F can be introduced in the summations over the Matsubara frequencies.
The ω n sums in equations (2) and (3) converge within an energy scale of the order of a few times (∼ 5 to 10) ω D , that is the upper characteristic frequency of the electron-phonon spectral function α 2 F k,k ′ (ω). This enables to set a cutoff Matsubara point ω nc , above which one can safely assume Z(k, iω n ) = 1 and ϕ ph (k, iω n ) = 0. Based on this observation, we rewrite equation (4) in the following form: where the summation for | ω n ′ |> ω nc in the second term has been simplified by replacing Z(k, iω n ) with unity and ϕ(k, iω n ) with ϕ c (k). We note that the reduced Matsubara-frequency dependence of the argument allows for a semianalytic evaluation of this sum by using the relation: Handling equation (9) by means of equation (10) leads to the expression: where the remaining Matsubara frequency summation is restricted, for both phonon and Coulomb contributions, to a narrow energy window. Equation (11) represents a huge simplification, because it effectively separates the low-energy scale at which the electron-phonon coupling is active, from the high-energy scale associated with Coulomb renormalization effects. Concerning the summation over k ′ , this can be carried out differently depending on the value of ε k ′ . The integration over low-energy states would benefit from the use of smart interpolation algorithms, such as that developed by R. Margine and coworkers [26][27][28], which allows for a very fine yet efficient sampling of the electron-phonon scattering processes near the Fermi surface. Less accuracy is needed while integrating over high-energy states, as these essentially contribute via Coulomb effects, which are weakly k-dependent. Hence, the summation over k ′ can be limited to a sparse set of points, e.g. by means of the energy-dependent random sampling used in SCDFT [13]. Actually, one might even resort to a simplified isotropic treatment of the Coulomb interaction (see section 2.2), justified by the fact that anisotropy effects in W k,k ′ are expected to average out over a large energy scale.

Isotropic approximation
Despite the simplifications introduced so far, solving the Eliashberg equations still involves the significant cost of performing a k-integration able to resolve the properties of the electron-phonon coupling near the Fermi surface, and to account for the full extent of the Coulomb interaction. In order to overcome this numerical issue, most Eliashberg codes resort to a simplified approach by neglecting the k-dependence in both λ k,k ′ and W k,k ′ . For the electron-phonon coupling, the isotropic approximation is defined by averaging equation (7) over the Fermi surface: which, through equation (6), yields the isotropic kernel λ(iν n ). Results show that this approximation is quite accurate for the estimation of T c . In fact, the anisotropy in the electron-phonon coupling turns out to have no relevant effect on T c for the majority of bulk superconductors, possibly being truly essential only to describe MgB 2 [29][30][31]. On the other hand, it affects, often crucially, all those superconducting properties which depend on the excitation spectrum (such as tunneling behavior, thermodynamic properties and magnetic response) [32][33][34].
Handling the k-dependence of the Coulomb interaction is a more subtle problem. Formally, this dependence cannot be simply neglected because the integration in equation (4) would diverge logarithmically. For this reason the µ * approach, which involves the isotropic approximation, relies on the introduction of an arbitrary high-energy cutoff (usually chosen as E F ) [4,5,35,36]. The most sensible strategy [8,10,31] is to approximate W k,k ′ by its average over surfaces of constant energy (ε) in k-space, i.e.: where N(ε) = ∑ k δ(ε k − ε). By averaging over k equations (2), (3) and (11), we obtain the following isotropic Eliashberg equations: where and The resulting scheme is formally similar to the static limit of the hybrid SCDFT-Eliashberg method recently introduced in [10]. The main difference lies in equations (14) and (15), where the energy dependence of ϕ c (ε) and N(ε) has been neglected to allow for an analytic integration over ε, as in usual applications of conventional Eliashberg theory [4]. For the numerical tests discussed in section 3 we solve equations (14)- (18), with the only exception of MgB 2 , for which we use a 2-band generalization [37].

Superconducting properties of layered compounds within ab-initio Eliashberg theory
For most conventional superconductors the role of Coulomb interactions is almost trivial, as it is reduced to a scaling of the superconducting gap at the Fermi energy. This is due to the fact that by integrating W(ε, ε ′ ) over a large energy scale, specific material properties are often washed out and the resulting effect is nearly material independent. The striking evidence for this is that the McMillan formula for T c , which depends on the single Coulomb parameter µ * , works reasonably well for many materials, with the only caveat that µ * ≃ 0.1 should be used for sp-electron systems and µ * ≃ 0.14 for d-electron systems [4,36]. However, peculiarities in the electronic band structure may lead to significant exceptions. An interesting case recently discussed are Chevrel phases, where the band character changes abruptly close to the Fermi level, significantly weakening the Coulomb renormalization [19], and for which a naive application of the McMillan formula gives a T c more than 50% larger than the experimental value. Anomalous Coulomb effects are likely to occur when the Fermi level is close to Van Hove singularities causing W(ε, ε ′ ) to suddenly drop or increase. A typical example is the proximity of the Fermi level to a band gap, as in weakly doped insulators [48,49], or to a sharp peak in the density of state, like in high pressure sulphur hydride [7,50]. From purely theoretical arguments, an enhancement of Coulomb interactions is more likely to occur in low dimensional systems, because of the reduced metallic screening [51]. In this work we test our method on four layered systems, which have been formerly investigated within SCDFT [16,21,31,34,38,52,53]. The coupling parameters used in the calculations are collected in table 1: electron-phonon couplings are taken from the literature [31,38,39], while the Coulomb interaction W is computed within the RPA approximation using the Elk code [54].

CaC 6
CaC 6 is an intercalated graphite compound with an experimental critical temperature of 11.5 K [40]. Superconductivity arises from the strong electron-phonon coupling provided by both C and Ca phonon modes [55]. This coupling is strongly anisotropic with C-and Ca-related phonons acting selectively on the multiple Fermi surface sheets of the system [22,56]. This leads to a strong anisotropy of the superconducting gap [34], but has a rather small effect on T c (of the order of 1 K). The µ * approach to the Coulomb interaction is relatively accurate, although a µ * = 0.14 needs to be used in order to reproduce the experimental T c , whereas the conventional value of 0.1 yields an overestimation of 40%. By solving equations (14)-(18) using the ab-initio calculated couplings of figures 1(a)-(d), we obtain an excellent T c estimate of 10.7 K. This value is essentially identical to that provided by SCDFT and the full-scale Eliashberg equations of [10]. The unusually strong Coulomb renormalization effect is due to the block-matrix structure of the screened Coulomb interaction W(ε, ε ′ ), which effectively restricts the Coulomb renormalization to a 5 eV energy window [56], instead of employing the full valence bandwidth of 20 eV. This block structure is clearly visible in figure 1(b).

MgB 2
MgB 2 is probably the most remarkable phononic superconductor. It has a T c of 39 K and the unique feature of two, electronically well-distinguished superconducting gaps. The highest gap occurs in the hole-doped strongly covalent σ bands, while the lowest gap arises from the π bands of the B electrons. Electronic [57,58] and phononic [59][60][61] properties of this material have been discussed in detail in a number of excellent theoretical studies. Due to the different nature of σ and π states crossing the Fermi level, the Coulomb interaction acts differently on the two bands. Both states contribute significantly to metallic screening [20,21], however, electrons in the more spatially constrained σ states feel a 30% stronger Coulomb repulsion. The most simplified approach to MgB 2 needs to account for this two-band anisotropic structure in both the electron-phonon coupling and the screened Coulomb repulsion. Solving the Eliashberg equations with the coupling functions shown in figures 1(f)-(p), we obtain a T c of 33.5 K, which is about 15% smaller than the experimental value. A similar underestimation resulted from previous SCDFT calculations [31]. We note that such an error in the predicted T c by means of ab-initio methods is relatively common, and overall quite good, owing to the number of approximations involved. In the present case this underestimation is probably related to some inaccuracy in the calculated electron-phonon coupling. In fact, the average coupling that we compute (λ = 0.7) turns out to be smaller than that provided by a likely more accurate Wannier-function approach [62].

YNi 2 B 2 C
The next system in our set is the quaternary borocarbide YNi 2 B 2 C, whose crystal structure is shown in figure 1(o). The peculiarity of this material is the pinning of the Fermi energy to a Van Hove singularity, which leads to an extremely sharp peak (less than 100 meV wide) in the density of states ( figure 1(n)). This unique feature of YNi 2 B 2 C challenges the validity of the adiabatic approximation [1], the neglect of dynamical Coulomb effects [10,53], and the use of an energy independent electron-phonon coupling. Nevertheless, we have tested our approach on this system. The first problem that one encounters is the choice of the lattice parameters. Most ab-initio simulations are either carried out at the ab-initio relaxed lattice (typically in the local or semi-local density approximation [25]) or at the experimental lattice. Usually, the latter approach is slightly more accurate, but results do not depend critically on the choice. However, it turns out that for YNi 2 B 2 C this is not the case. By simulating the system at the experimental lattice we obtain an integrated electron-phonon coupling λ = 0.86, whereas using the LDA relaxed lattice gives λ = 0.56, that is a Table 1.    Li0.5ZrNCl ((p)-(t)). Panels (a),(f),(k),(p): α 2 F electron-phonon function [31,38,39]. Panels (b),(g),(l),(q): Coulomb function W(ε, ε ′ ); a diagonal cut (ε = ε ′ ) and a Fermi surface cut (ε ′ = 0) of W(ε, ε ′ ) are displayed in panels (c),(h),(m),(r). Panels large discrepancy. We do not delve here into this complex issue, and choose to use the electron-phonon coupling computed at the experimental lattice (shown in figure 1(o)), as it best describes the vibrational spectrum [39,63,64]. Consistently, we compute the static screened Coulomb interaction with the same structure. As shown in figures 1(l) and (m) the Coulomb repulsion is relatively smooth, with decreasing matrix elements at higher energies, not unlike the electron gas model [65]. However, since W(ε, ε ′ ) enters equation (16) jointly with the density of states, it plays a non trivial role in determining T c . We obtain a T c of 11.3 K, which underestimates the experimental value, reported to be within the range 14.6 K [44]-15.6 K [45]. A simple McMillan estimation with µ * = 0.1 predicts a T c of 15.2 K, in better agreement with experiments, even though, due to the d character of the density of states at the Fermi level, a larger µ * is conventionally recommended [4].

Li 0.5 ZrNCl
Lastly, we consider the lithium-doped ternary transition-metal nitride halide ZrNCl. Undoped ZrNCl is a layered semiconductor. The intercalation of lithium impurities into the Van der Waals gap (the interlayer space) leads, for x = 0.5, to the crystal structure shown in figure 1(t). This material has been studied thoroughly by Akashi and coworkers [38] within the SCDFT framework. Here, we use the electron-phonon coupling extracted from their reference and shown in figure 1(p).
The Coulomb interaction matrix that we computed is shown in figures 1(q) and (r). Compared to MgB 2 and CaC 6 , this function is significantly more structured, owing to the fact that, this system being a doped insulator, the metallic screening is poor. In fact at the band edges (that is close to the points where the density of states tends to zero in figure 1(r)), W(ε, ε ′ ) shows very high peaks caused by the 1/q divergence of the bare Coulomb interaction. Since the Fermi level is very close to one of the peaks, the Coulomb repulsion is strong. A µ * = 0.2 is required in the McMillan formula to reproduce the extrapolated experimental [47] T c of approximately 11 K, whereas the standard µ * = 0.1-0.14 gives a T c in the range 18-22 K, that is almost twice the experimental value. On the contrary, our Eliashberg approach yields an excellent T c estimate of 9.5 K. However, one should point out that due to the peaked-structure of W(ε, ε ′ ), the estimation of T c is very sensitive to parameters like the position of the Fermi level and properties of the dielectric function. Therefore, in a doped insulator of this kind, a 20% error in T c is usually expected [48]. As a test, we have used the same input functions to compute the critical temperature within SCDFT. This latter approach gives a T c of 13.1 K, that is about 15% higher than the experimental value. We observe that previous SCDFT-based calculations on this system [38] predicted a very low T c , highlighting a large discrepancy between theory and experiments. We can now assert that the problem had to be ascribed to inaccuracy of the functional available at the time [11].

Conclusions
We have presented a minimal yet efficient approach to include static Coulomb interactions in Eliashberg theory from first principles, which avoids the arbitrariness introduced by the use of the adjustable parameter µ * . The method is formally equivalent to the static limit of the ab-initio Eliashberg theory recently proposed in [10], but, unlike the more complex original scheme, allows for a straightforward implementation in conventional Eliashberg codes. In the derived equations, the effects of the Coulomb interaction outside the phonon frequency range are evaluated analytically, and the numerical integration is reduced, as for the standard approach, to the energy window 0 → ω nc ∼ 10 ω D . Both parts include the matrix elements of the screened Coulomb interaction, which can be computed using the random phase approximation for the dielectric matrix. In practice, this is the main additional computational step compared to conventional implementations. We have tested our method on the layered superconducting materials CaC 6 , MgB 2 , YNi 2 B 2 C and Li-doped β-ZrNCl, which present a non-monotonous behavior of the Coulomb repulsion, such to challenge the validity of the µ * approximation. The new approach turns out to be accurate as the estimated critical temperatures are consistent with SCDFT values. The agreement with experiments is excellent in terms of both T c 's and gaps for CaC 6 , MgB 2 and Li-doped β-ZrNCl, with a significant improvement in prediction quality with respect to the µ * formula. We observe, however, a large error (25%) in the predicted T c of YNi 2 B 2 C. We believe that this anomaly should be further investigated, as it stems from a difficulty in determining the electron-phonon coupling, and might point to strong non-adiabatic effects related to the presence of a sharp Van Hove singularity at the Fermi level.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).