Electron-phonon beyond Fr\"ohlich: dynamical quadrupoles in polar and covalent solids

We include quadrupolar fields beyond the Fr\"ohlich interaction in the first-principles electron-phonon vertex of semiconductors, and demonstrate their importance in calculations of carrier mobilities. Without such quadrupolar corrections, jump discontinuities for $\mathbf{q} \rightarrow \mathbf{0}$ remain in the short-range components. They lead to Gibbs oscillations in the interpolant, and affect the accuracy of the physical results. We apply our formalism to Si (non-polar), GaAs and GaP (polar). Electron mobilities converge much faster with the initial ab initio $\mathbf{q}$-mesh when dynamical quadrupoles are properly included.

The electron-phonon (e-ph) interaction plays a key role in the description of a variety of physical phenomena such as electronic transport, phonon-assisted light absorption, and phonon-mediated superconductivity [1]. In state-of-the-art ab initio methods, the e-ph coupling is described within density functional theory (DFT) by expanding the Kohn-Sham (KS) effective potential [2] in the nuclear displacements, and the vibrational properties are obtained with density-functional perturbation theory (DFPT) [3,4]. This DFPT-based computational scheme enables the calculation of screened e-ph matrix elements on a microscopic level with ab initio quality [1]. However, e-ph related properties require an accurate description of the coupling on very dense grids in the full Brillouin zone (BZ) thus rendering direct ab initio e-ph computations in real materials impracticable.
For this reason, different methods have been proposed to interpolate the e-ph matrix elements from initial coarse grids to dense ones. These approaches rely on localized basis sets such as maximally-localized Wannier functions [5][6][7] or atomic orbitals [8] to perform a Fourier interpolation of the e-ph vertex. Other works suggested to employ Fourier transforms to interpolate the local part of the scattering potentials and then use Bloch states to obtain the e-ph matrix elements on dense grids [9,10].
Despite formal differences in the treatment of the electron wavefunctions, all these approaches rely on the spatial localization of the scattering potentials for obtaining a reliable interpolation in the phonon wavevector q. In Fourier-based interpolation schemes, indeed, the signal in real space is assumed to be negligibly small beyond a certain distance so that a faithful interpolation in qspace can be obtained with a relatively coarse sampling of the BZ. In semiconductors and insulators, however, the incomplete screening of the external potential generated by the atomic displacements leads to long-range (LR) interactions. In polar materials, these interactions manifest in the long-wavelength limit (q → 0) by the LO-TO splitting of the optical frequencies [11] and the Fröhlich divergence of the e-ph matrix elements [12] that are not amenable to Fourier interpolation due to their non-analytic behavior. The treatment of the LR dipoledipole interaction in the Fourier interpolation of the dynamical matrix was developed in the early days of DFPT and the methodology is well documented [3,4]. Only recently, it was shown how to generalize the Fröhlich contribution to e-ph matrix elements to anisotropic materials and how to integrate it with Wannier-based interpolation schemes [6,13].
In this letter, we go beyond these contributions by demonstrating the importance of the next-to-leading order terms derived by Vogl [14]. We develop an approach to interpolate the short-range (SR) part of the scattering potentials on dense q-grids while the LR contributions are expressed in terms of the high-frequency dielectric tensor, Born effective charges (dipole potential), dynamical quadrupoles computed from first principles using the theory of spatial dispersion [15,16] (quadrupole potential) and the response to a homogeneous static electric field (local-field potential, also quadrupolar) already available in many DFPT implementations [3,4]. In nonpolar materials, the quadrupolar fields represent the leading contribution to the LR part, that should therefore be taken into account for accurate interpolations. In polar semiconductors, quadrupole terms are in principle small when compared to the divergent Fröhlich interaction but their non-analytical behavior for q → 0 implies a LR behavior in real space that will spoil any Fourier-based interpolation technique if not handled properly. Our formalism generalizes the previous approaches [6,13] by including contributions that are crucial to achieve reliable and accurate e-ph calculations in semiconductors. We demonstrate the importance of the quadrupole corrections by analyzing the convergence rate of electron mobilities in Si and GaAs. Additional results for GaP can be found in the Supplemental Material [17]. Without quadrupoles, carrier mobilities are plagued by an exceedingly slow convergence with respect to the initial ab initio q-mesh. Stable results can be achieved with reasonable q-meshes only when quadrupolar fields are separated from the SR part before the interpolation. Further details concerning our implementation including additional convergence studies and results for the electron mobility of GaP are given in our accompanying manuscript [18].
The key ingredients in e-ph computations are the coupling matrix elements g mnν (k, q) = ψ mk+q |∆ qν V |ψ nk with ψ nk the nk Bloch state and ∆ qν V the first-order variation of the KS potential V due to a phonon mode of wavevector q and branch index ν [1]. The scattering potential is defined as with ω qν the phonon frequency and e κα,ν (q) the α-th Cartesian component of the phonon eigenvector for the atom κ of mass M κ in the unit cell. V κα,q (r) is the latticeperiodic scattering potential computed with DFPT [19]. Following the same approach as Verdi et al. [13], the e-ph scattering potential is separated into SR (S) and LR (L) contributions: V L is supposed to include all the LR components so that V S is smooth in q-space and therefore tractable with Fourier interpolation. In the interpolation algorithm, V L is first subtracted from the DFPT potentials, then the Fourier interpolation is performed on the SR part only, finally V L evaluated at the interpolated q-point is added back [18]. In polar materials, the leading contribution to the LR part stems from the diverging Fröhlich-like potential [13]: with Ω the unit cell volume, G a reciprocal lattice vector, Z * κ the Born effective charge tensor, ∞ the highfrequency dielectric tensor and τ κ the position of the κth atom in the unit cell. The summation over repeated indices (β, η, δ and δ ) is implied in Eq. (3) and in the following, unless the sum is explicitly written.
Most investigations so far have been focusing on the treatment of Eq. (3) because the (integrable) divergence for q → 0 is expected to give an important contribution to the BZ integrals. However, as discussed by Vogl [14] and rederived in a DFPT context in Ref. [18], a careful analysis of the asymptotic behavior of the scattering potential in the long-wavelength limit reveals the presence of additional non-analytical terms besides Eq. (3). These are finite for q → 0 but their non-analytic behavior (discontinuities) will affect the spatial decay of the scattering potentials even when the dipole interaction given by Eq. (3) is properly accounted for. Both dipole and quadrupole terms can be included in the LR potential using the generalized expression [18]: where Q κα is the dynamical quadrupole tensor, Q is the electronic charge in atomic units, that is, −1, and v Hxc,Eγ is the change of the Hartree and exchange-correlation potential with respect to the electric field E in Cartesian coordinates (see Ref. [18] for the derivation). In Eq. (4), the term related to Q κα is non-zero even in non-polar semiconductors while the contribution associated to E is present only in systems with non-zero Born effective charges.
In the polar systems that we investigated so far, we observed that Q κα gives the most important contribution to Eq. (4) when compared to the E term [17]. For instance, ignoring the E term changes the electron mobility in GaAs by 0.1% and by 0.01% in GaP (see next paragraphs for more details about the mobility computations). Qualitatively, this behavior can be rationalized as follows. In the ionic limit, the potential change due to E can be understood as a local linear potential created around each atom. Its average must be zero, but there will be regions of constant positive shift and zones of constant negative shift. If the orbitals are centered on each atom (so that the charge density is symmetric around it), the symmetry of the charge and the anti-symmetry of the potential give a vanishing contribution. In other words, after the Fröhlich divergence, the spatial decay of the scattering potentials is governed by the dynamical quadrupoles Q κα at least as far as the materials considered here are concerned. Strictly speaking, it is not possible to generalize this observation on the predominance of Q κα to other non-polar systems although it is reasonable to expect that the same picture will be valid.
Note also that Q κα introduces additional dipolequadrupole and quadrupole-quadrupole terms at the level of the dynamical matrix [16] that are important to improve the accuracy of the interpolated phonon spectrum. Our calculations in Si, GaAs and GaP revealed that these extra terms in the Fourier interpolation of the dynamical matrix can accelerate the convergence of mobilities at the level of ∼1%, an effect that is smaller when compared to the error introduced by the Fourier interpolation of the e-ph potentials (discussed below) and yet larger than the correction due to the electric field. In the following, we will therefore focus on the contribution of Q κα while the effects associated to the electric field, the dipole-quadrupole and the quadrupole-quadrupole interactions in the dynamical matrix are discussed in [17].
We begin by analyzing the effect of Q κα on the interpolation of the scattering potentials and of the matrix elements. Then we discuss the convergence behavior of electron mobilities with respect to the ab initio q-grid.
The numerical values of Q κα are computed using the recent implementation of Royo et al. [15] that has been integrated with the e-ph part of Abinit [18,20,21].
In Figs. 1(a) and (b) we plot the average over the unit cell of the lattice-periodic part of the scattering potential, for selected atomic perturbations in Si and GaAs, along a high-symmetry q-path. The exact DFPT results (blue lines) are compared with those obtained with the models with (green) or without (red) quadrupole corrections of Eqs. (3) and (4) [22]. The other potentials, for all perturbations, are reported in the Supplemental Material [17]. In Si, the Born effective charges are zero and the imaginary part of the potential does not diverge for q → 0 (see dashed lines in Fig. 1(a)). In GaAs, the Fröhlich-like model in Eq. (3) correctly describes the divergence of the imaginary part of the potential close to Γ (see green dashed line in Fig. 1(b)). In both materials, however, the real part of the potential (solid lines) presents a jump discontinuity for q → 0. Note that the Fröhlich term alone does not capture this non-analyticity. On the contrary, if the quadrupole Q κα contributions are included through Eq. (4), the LR model reproduces these jumps as shown by the solid green lines in Figs. 1(a) and (b). Figure 1(c) shows the real part of the average of the scattering potential in Si interpolated from a 9 × 9 × 9 q-point grid onto the same q-path as in Fig. 1(a). If the Q κα terms are not removed from the input DFPT potentials, the Fourier interpolation introduces unphysical sharp oscillations for small q (see red line, FI). The correct behavior is properly reproduced only when Q κα is included (see green line, FI+Q). Small oscillations are still appreciable around Γ when a 9 × 9 × 9 q-mesh is used but these deviations have a limited effect on the final electron mobilities.
It is worth stressing that these considerations hold for any approach employing Fourier-based interpolations. The discontinuity of the matrix elements at Γ and the discrepancy between the interpolant and the exact DFPT results in the region around Γ have already been noticed for Si and diamond [8]. The resulting error was considered harmless under the assumption that it is always possible to improve the accuracy of the Fourier interpolation by densifying the initial ab initio q-mesh [8]. Unfortunately, this assumption does not hold in the presence of non-analytical behavior since, strictly speaking, an infinite number of Fourier components (i.e. an infinite number of real space lattice vectors R) would be needed to represent the discontinuous signal. To confirm our hypothesis, we modified the EPW code [7] to include the Q κα terms in the interpolation of the e-ph matrix elements in the Wannier representation in Si [23]. Figure 2 shows the e-ph matrix elements connecting the electron state at the valence band maximum at k = Γ to the other highest states of the valence band through the highest phonon mode at L as in Ref. [8]. The DFPT results (blue solid line) and the Wannier-interpolated results (red and green dashed lines) are compared in this figure. The interpolation of the matrix elements using the standard approach (i.e. without Q κα ) leads to systematic errors in the region around the Γ point associated to electronic transitions with small momentum transfer.
Including the dynamical quadrupoles in the model significantly improves the quality of the interpolant around the Γ point. Similar behaviour is also expected in e-ph calculations based on supercells and finite differences [10]. At this point, one question arises spontaneously: can we ignore such physical effects, and simply interpolate with an exceedingly dense grid around Γ? To answer this question, we need to quantify the error introduced by these spurious oscillations in the final physical results. Because of intra-valley transitions, the region around Γ is one of the most important for the description of eph scattering processes [17,18]. An accurate description in this region is crucial for reliable calculations of the phonon-induced electron linewidths, with Ω BZ the BZ volume, n qν and f mk+q the Bose-Einstein and Fermi-Dirac occupation functions and ε nk the energy of the electronic state nk. These linewidths represent the phonon-dependent ingredients needed to compute carrier mobilities within the self-energy relaxation time approximation [1,18,24,25]: where v nk,α is the α-th component of the velocity operator [18]. Figures 3(a) and (b) report the electron mobilities in Si and GaAs computed using plane waves with Abinit as described in Ref. [18], as a function of the initial q-grid. For each initial ab initio DFPT mesh, we interpolate the scattering potentials to obtain the lifetimes on k-and q-point meshes that are dense enough to reach convergence in the mobility within 5%. If the terms associated to the dynamical quadrupoles are included in the LR part of the potentials, a 9×9×9 (6×6×6) q-mesh is already sufficiently dense in Si (GaAs) for the correct interpolation of the potentials (green lines) whereas without Q κα the convergence is much slower and not even reached with a 21 × 21 × 21 q-mesh (red lines). Using q-grids typically reported in the literature (8 × 8 × 8) can therefore lead to significant errors at the level of the mobility. The error is around 10% in Si and goes up to 30% in GaAs where, due to the small effective mass, most of the scattering channels for electrons close to the band edge involve small momentum transfer. E-ph interpolations are usually performed starting from coarse q-meshes because both the memory requirements and the computational cost of the interpolation quickly increase with the number of q-points in the initial ab initio mesh hence it is not surprising that this behavior has been largely overlooked so far.
In conclusion, we have included quadrupolar fields beyond the Fröhlich interaction in the first-principles electron-phonon vertex of semiconductors. Their importance was demonstrated in calculations of mobilities. Accurate calculations of e-ph properties in semiconductors and insulators can be achieved with reasonably coarse ab initio q-meshes provided that a careful treatment of the long-range interaction associated to dynamical quadrupoles is properly taken into account. We presented a fully ab initio formalism that employs Q κα computed with DFPT to improve the description of the non-analytical behavior of the e-ph scattering potentials in the long-wavelength limit. Our approach improves state-of-the-art techniques that were mainly designed to cope with the Fröhlich-like interaction in polar materials. Since long-range macroscopic interactions play a key role in semiconductor physics, we believe that the quest for accurate ab initio descriptions of phonon-limited carrier mobilities should start from a proper treatment of these physical phenomena. used an energy cutoff of 25 Ry for the plane-wave basis set and a lattice parameter of 5.47Å. The matrix elements are interpolated starting from coarse Γ-centered 12 × 12 × 12 k-point and 6 × 6 × 6 q-point meshes. For the implementation of the quadrupolar corrections in EPW, we employed the quadrupole tensors computed by Abinit as the same computation is not available in QE. Because of the different pseudopotentials, lattice parameters, and computational settings, we had to adjust the value of Qκα to better describe the discontinuity, from 13.67 with LDA FHI [18]  As mentioned before, the term associated to E plays a very minor role if we focus on the average of the potential over the unit cell that corresponds to the G = 0 Fourier component.
To appreciate the effect of the electric-field term, we need to focus on the q-dependence of the Fourier components of the scattering potentials for small G = 0. In our tests we found that, for particular G-vectors, the quantity in Eq. S1 as a function of q presents (small) jump discontinuities for q → 0 and that the discontinuity is better described when the E-term is included in the LR model. The results are summarized in Figs. S9 and S10. We performed the same analysis for GaP (not shown) and found very similar behavior. At the level of the mobility, the inclusion of the electric-   the order of 1%, value that is anyway larger than the one observed from the E-term. We speculate that, in more "complicated" materials, the inclusion of these higher-order terms will play a more important role for the accurate description of the phonon dispersion in the long-wavelength limit. Figures S12 and S13 show the same comparison for GaP and GaAs.