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

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 ($\sim 5$ to 10) ω_D , that is the upper characteristic frequency of the electron–phonon spectral function $\alpha^2F_{k,k^{^{\prime}}}(\omega)$. This enables to set a cutoff Matsubara point $\omega_{n_c}$, above which one can safely assume $Z(k,i\omega_n) = 1$ and $\phi^{ph}(k,i\omega_n) = 0$. Based on this observation, we rewrite equation (4) in the following form:

$$\begin{align} \phi^c(k)& = -\frac{1}{\beta}\sum_{k^{^{\prime}}}W_{k,k^{^{\prime}}}\left[\sum_{n^{^{\prime}} = -n_c}^{n_c}\frac{\phi(k^{^{\prime}},i\omega_{n^{^{\prime}}})}{\Theta(k^{^{\prime}},i\omega_{n^{^{\prime}}})}\right.\nonumber\\ & \quad \left.+\sum_{|n^{^{\prime}}|\gt n_c}\frac{\phi^c(k^{^{\prime}})}{\omega_{n^{^{\prime}}}^2+\varepsilon_{k^{^{\prime}}}^2+{\phi^c}^2(k^{^{\prime}})}\right], \end{align} \tag{ 9 }$$

where the summation for $\mid \omega_{n^{^{\prime}}}\mid \gt \omega_{n_c}$ in the second term has been simplified by replacing $Z(k,i\omega_n)$ with unity and $\phi(k,i\omega_n)$ with $\phi^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:

$$\begin{equation} \frac{1}{\beta}\sum_n \frac{1}{\omega_n^2+A} = \frac{1}{2}\frac{\tanh\left[\frac{\beta\sqrt{A}}{2}\right]}{\sqrt{A}}. \end{equation} \tag{ 10 }$$

Handling equation (9) by means of equation (10) leads to the expression:

$$\begin{align} \phi^c(k) = &-\sum_{k^{^{\prime}}}W_{k,k^{^{\prime}}} \Biggl\{\frac{1}{2}\frac{\tanh\left[\frac{1}{2}{\beta \sqrt{\varepsilon_{k^{^{\prime}}}^2+{\phi^c}^2(k^{^{\prime}})}}\right] }{\sqrt{\varepsilon_{k^{^{\prime}}}^2+{\phi^c}^2(k^{^{\prime}})}}\nonumber\\ &+\frac{1}{\beta}\!\sum_{n^{^{\prime}} = -n_c}^{n_c}\!\left[ \frac{\phi(k^{^{\prime}},i\omega_{n^{^{\prime}}})}{\Theta(k^{^{\prime}},i\omega_{n^{^{\prime}}})}-\frac{\phi^c(k^{^{\prime}})}{\omega_{n^{^{\prime}}}^2+\varepsilon_{k^{^{\prime}}}^2+{\phi^c}^2(k^{^{\prime}})}\right]\Biggr\}, \end{align} \tag{ 11 }$$

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 $\varepsilon_{k^{^{\prime}}}$. 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–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^{^{\prime}}}$ are expected to average out over a large energy scale.

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 $\lambda_{k,k^{^{\prime}}}$ and $W_{k,k^{^{\prime}}}$. For the electron–phonon coupling, the isotropic approximation is defined by averaging equation (7) over the Fermi surface:

$$\begin{equation} \alpha^2F(\omega) = \frac{1}{N(0)^2} \sum_{k,k^{^{\prime}}}\alpha^2F_{k,k^{^{\prime}}}(\omega) \delta(\varepsilon_{k})\delta(\varepsilon_{k^{^{\prime}}}), \end{equation} \tag{ 12 }$$

which, through equation (6), yields the isotropic kernel $\lambda(i\nu_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₂ [29–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–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 $\mu^*$ 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^{^{\prime}}}$ by its average over surfaces of constant energy (ε) in k-space, i.e.:

$$\begin{equation} W(\varepsilon,\varepsilon^{^{\prime}}) = \frac{1}{N(\varepsilon)N(\varepsilon^{^{\prime}})} \sum_{k,k^{^{\prime}}} W_{k,k^{^{\prime}}}\delta(\varepsilon_k-\varepsilon)\delta(\varepsilon_{k^{^{\prime}}}-\varepsilon^{^{\prime}}), \end{equation} \tag{ 13 }$$

where $N(\varepsilon) = \sum_k \delta(\varepsilon_k-\varepsilon)$. By averaging over k equations (2), (3) and (11), we obtain the following isotropic Eliashberg equations:

$$\begin{align} Z(i\omega_n)& = 1+\frac{\pi}{\beta}\sum_{n^{^{\prime}} = -n_c}^{n_c}\frac{\lambda(i\omega_n-i\omega_{n^{^{\prime}}})\omega_{n^{^{\prime}}}Z(i\omega_{n^{^{\prime}}})}{\omega_n\sqrt{\left[\omega_{n^{^{\prime}}} Z(i\omega_{n^{^{\prime}}})\right]^2+\phi^2(0,i\omega_{n^{^{\prime}}})}}, \end{align} \tag{ 14 }$$

$$\begin{align} \phi^{ph}(i\omega_n)& = \frac{\pi}{\beta}\sum_{n^{^{\prime}} = -n_c}^{n_c}\frac{\lambda(i\omega_n-i\omega_{n^{^{\prime}}})\phi(0,i\omega_{n^{^{\prime}}})}{{\sqrt{\left[\omega_{n^{^{\prime}}} Z(i\omega_{n^{^{\prime}}})\right]^2+\phi^2(0,i\omega_{n^{^{\prime}}})}}}, \end{align} \tag{ 15 }$$

$$\begin{align} \phi^c(\varepsilon) = & -\int d\varepsilon^{^{\prime}} W(\varepsilon,\varepsilon^{^{\prime}})N(\varepsilon^{^{\prime}}) \Biggl\{\frac{1}{2}\frac{\tanh\left[\frac{\beta}{2}\sqrt{\varepsilon^{^{\prime} 2}+{\phi^c}^2(\varepsilon^{^{\prime}})}\right]}{\sqrt{\varepsilon^{^{\prime} 2}+{\phi^c}^2(\varepsilon^{^{\prime}})}} \nonumber\\ &+ \frac{1}{\beta}\sum_{n^{^{\prime}} = -n_c}^{n_c}\left[ \frac{\phi(\varepsilon^{^{\prime}},i\omega_{n^{^{\prime}}})}{\Theta(\varepsilon^{^{\prime}},i\omega_{n^{^{\prime}}})}- \frac{\phi^c(\varepsilon^{^{\prime}})}{\omega_n^2+\varepsilon^{^{\prime} 2}+{\phi^c}^2(\varepsilon^{^{\prime}})} \right] \Biggr\}, \end{align} \tag{ 16 }$$

where

$$\begin{equation} \phi(\varepsilon,i\omega_n) = \phi^{ph}(i\omega_n)+\phi^c(\varepsilon), \end{equation} \tag{ 17 }$$

and

$$\begin{equation} \Theta(\varepsilon,i\omega_n) = \left[\omega_n Z(i\omega_n)\right]^2+\varepsilon^2+\phi^2(\varepsilon,i\omega_n). \end{equation} \tag{ 18 }$$

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 $\phi^c(\varepsilon)$ and $N(\varepsilon)$ 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₂, for which we use a 2-band generalization [37].

CaC₆ 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 $\mu^*$ approach to the Coulomb interaction is relatively accurate, although a $\mu^* = 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(\varepsilon,\varepsilon^{^{\prime}})$, 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).

Zoom In Zoom Out Reset image size

MgB₂ 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–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₂ 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].

The next system in our set is the quaternary borocarbide YNi₂B₂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₂B₂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₂B₂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 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(\varepsilon,\varepsilon^{^{\prime}})$ 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 $\mu^{*} = 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 $\mu^{*}$ is conventionally recommended [4].

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₂ and CaC₆, 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(\varepsilon,\varepsilon^{^{\prime}})$ 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 $\mu^*$ = 0.2 is required in the McMillan formula to reproduce the extrapolated experimental [47] T_c of approximately 11 K, whereas the standard $\mu^* = 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(\varepsilon,\varepsilon^{^{\prime}})$, 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].