First-principles predictions of Hall and drift mobilities in semiconductors

Significant progress on parameter-free calculations of carrier mobilities in real materials has been made during the past decade; however, the role of various approximations remains unclear and a unified methodology is lacking. Here, we present and analyse a comprehensive and efficient approach to compute the intrinsic, phonon-limited drift and Hall carrier mobilities of semiconductors, within the framework of the first-principles Boltzmann transport equation. The methodology exploits a novel approach for estimating quadrupole tensors and including them in the electron-phonon interactions, and capitalises on a rigorous and efficient procedure for numerical convergence. The accuracy reached in this work allows to assess common approximations, including the role of exchange and correlation functionals, spin-orbit coupling, pseudopotentials, Wannier interpolation, Brillouin-zone sampling, dipole and quadrupole corrections, and the relaxation-time approximation. A detailed analysis is showcased on ten prototypical semiconductors, namely diamond, silicon, GaAs, 3C-SiC, AlP, GaP, c-BN, AlAs, AlSb, and SrO. By comparing this extensive dataset with available experimental data, we explore the intrinsic nature of phonon-limited carrier transport and magnetotransport phenomena in these compounds. We find that the most accurate calculations predict Hall mobilities up to a factor of two larger than experimental data; this could point to promising experimental improvements in the samples quality, or to the limitations of density-function theory in predicting the carrier effective masses and overscreening the electron-phonon matrix elements. By setting tight standards for reliability and reproducibility, the present work aims to facilitate validation and verification of data and software towards predictive calculations of transport phenomena in semiconductors.


I. INTRODUCTION
The ability of metals and semiconductors to transport electrical charges is a fundamental property in manifold applications, ranging from solar cells, light-emitting devices, thermoelectric, transparent conductors, photodetectors, photo-catalysts, and transistors [1][2][3][4][5] . Predicting transport properties from first-principles 6 thus offers a powerful tool for the design of new materials and devices.
The aim of this work is to to assess the reliability and predictive power of novel first-principles methods, using III-V and group IV semiconductors as a representative benchmark. Already in 1966, in a seminal work Cohen and Bergstresser 7 computed the band structure of 14 semiconductors with the diamond and zincblende structure. In 2013, Malone and Cohen 8 repeated this task computing the quasiparticle bandstructure in the manybody GW approximation, including spin-orbit coupling (SOC) effects. Last year, Miglio et al. 9 studied the zeropoint renormalization of the electronic band gap of 30 semiconductors, mostly of the zincblende, wurtzite and rocksalt crystal structure. Here, we deliver an accurate and efficient calculation method for carrier mobilities, including spin-orbit coupling, dynamical quadrupoles and magnetic Hall effects, focusing on prototypical semiconductors. Our benchmark includes the following ten semi-conductors: diamond, Si, GaAs, 3C-SiC, AlP, GaP, cubic BN, AlAs, AlSb and SrO. For five of these no firstprinciples calculations of carrier mobility have been reported to date.
Due to the complexity of computing carrier mobilities fully from first principles, the first calculation appeared only in 2009 10 46 .
In this work, we perform an in-depth analysis of the predictive accuracy of the first-principles Boltzmann transport equation (BTE) for computing the carrier drift mobility and the Hall mobility, and we compare directly to experimental data. In brief, we find that spin-orbit coupling plays an important role for hole mobilities, that a local approximation to the velocity matrix elements (also used by us in earlier work) is not accurate enough, and that optimizing the construction of the Wigner-Seitz cell used in Wannier interpolations can accelerate the convergence of the Brillouin zone integrals. We confirm the recent findings 16,22,25 that dynamical quadrupoles, the next order of correction to dynamical dipoles, can significantly affect carrier mobilities, and we show that this effect arises predominantly from the change of the vibrational eigenmodes induced by the quadrupole term of the dynamical matrix.
In this work we also find that Wannier-function interpolations (especially those associated with the conduction bands) converge more slowly than previously estimated with respect to the sampling of the coarse grid employed in the Wannier representation. As a result, in order to obtain carrier mobilities with an accuracy of 1% for fixed fine grids, it is often necessary to use coarse grids including up to 18 3 points, still providing a speed-up of two to three orders of magnitude with respect to brute force sampling. We also find that a good indicator of convergence is provided by the mobility effective mass, as defined below. The sampling of Wannier-interpolated quantities require between 80 3 and 250 3 k and q grids, with the exception of the electron mobility of GaAs which requires a 500 3 grid due to its very small electron effective mass. We also show that the mobility and Hall factor converge linearly with grid spacing, allowing for a simple extrapolation 47 which achieves an accuracy of 1%.
This work answers crucial methodological and physical questions. On the technical side, it establishes welldefined criteria for high-accuracy calculations of transport coefficients, clarifying the role and the level of (i) the approximations in the underlying first-principles calculation to extract materials parameters (exchange-andcorrelation functionals, pseudopotentials, lattice parameters and SOC), (ii) the interpolation procedures for quasiparticle dispersions and interactions, and (iii) the numerical techniques that ensure the efficient convergence of transport coefficients. On the physical side, the ability to describe with high accuracy transport phenomena provides detailed access to the carrier dynamics and to the inherent limits of the electrical transport properties of key prototypical semiconductors. Crucially, since usually it is the Hall mobility that is measured rather than the drift mobility, this work allows a detailed and unambiguous comparison with experimental data.
The manuscript is organized as follows. In Section II we briefly present the theory for computing phononlimited drift mobilities. We then extend the theory to include a small finite external magnetic field, to access the Hall mobility. We also show how to efficiently compute mobility using interpolations based on maximallylocalized Wannier functions (MLWF) 48 . In particular, we discuss how to compute exact velocities and estimate quadrupole tensors, and include dipole-dipole, dipolequadrupole and quadrupole-quadrupole corrections during the Fourier interpolation of the interatomic force constant, as well as how to include dynamical quadrupoles in the interpolation of the electron-phonon matrix element. Finally, we discuss the use of adaptive broadening for improved convergence.
In Section III, we present the computational methodology. We first report the methodology and parameters used, then outline the Wannier interpolation procedure for electrons, phonons, and electron-phonon matrix ele-ments. Next, we report our interpolation of the electronphonon matrix elements and deformation potential including dynamical quadrupole, and show that the interpolated values reproduce density functional perturbation theory calculations. We then discuss typical coarse and fine grid convergence rates. We conclude the section by showing that the mobility and the Hall factor converge linearly with the spacing of the fine grid, allowing for a linear extrapolation to the limit of exact momentum integration.
In Section IV we present and discuss our results, starting with analysis of the band structures, phonon dispersions, and effective masses. We then assess the quality of various popular approximations, such as the neglect of the non-local pseudopotential contribution to the electron velocity, the neglect of SOC, the neglect of quadrupole correction, the use of the relaxation time approximation, the effect of the isotropic Hall factor approximation, the effect of the exchange-correlation functional, the effect of lattice parameters and pseudopotential choice. In addition, we identify the dominant scattering mechanisms responsible for limiting the carrier mobility. Finally, we conduct an in-depth comparison with available experimental mobility data, and assess the overall predictive power of the BTE. We offer our conclusions in Section V.

II. THEORY
In this Section we first present the main equations for calculating the drift mobility, and their generalization to the case of vanishing magnetic field, yielding the Hall mobility. We then present the Wannier interpolation scheme employed here, and discuss how to obtain accurate electron velocities. We also discuss the multipole expansion of the long-range part of the dynamical matrix and electron-phonon matrix elements. Finally, we describe how we use adaptive broadening in practical calculations.

A. Drift mobility
The low-field phonon-limited charge carrier mobility is calculated as [10][11][12][13][14][15][15][16][17][18][19][20][21][22][23][24][25][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45][46] : Here α, β run over the three Cartesian directions and ∂ E β f nk ≡ (∂f nk /∂E β )| E=0 is the linear variation of the electronic occupation function f nk in response to the electric field E, V uc is the unit cell volume, Ω BZ the first Brillouin zone volume, v nkα = −1 ∂ε nk /∂k α is the band velocity for the Kohn-Sham state ε nk , and n c = (1/V uc ) n (d 3 k/Ω BZ )f 0 nk is the carrier concentration. f 0 nk is the Fermi-Dirac occupation function at equilibrium (in the absence of fields). We follow the notation of Refs. 6,14. The BTE describes the detailed balance between the carriers' populations moving through phase space under the action of the driving electric field E and the carriers scattered by phonons. It can be derived from the Kadanoff-Baym equation of motion by approximating the Hartree and exchange-correlation potential, taking the diagonal Bloch state projection and assuming a spatially homogenous electric field 6 . One then obtains the quantum time-dependent BTE that can be further simplified by considering a time-independent external field (DC), static electron-one-phonon interactions and adiabatic phonons to give: where n qν is the Bose-Einstein distribution. The electron-phonon matrix elements g mnν (k, q) are the amplitude for scattering from an initial state nk to a final state mk + q via the emission or absorption of a phonon of frequency ω qν . We obtain ∂ E β f nk required for Eq. (1) by taking the field derivative of Eq. (2) at vanishing field. This yields the linearized Boltzmann transport equation (BTE) 6,14 : ∂ E β f nk =ev nkβ ∂f 0 nk ∂ε nk τ nk + 2πτ nk mν d 3 q Ω BZ |g mnν (k, q)| 2 × (n qν + 1 − f 0 nk )δ(ε nk − ε mk+q + ω qν ) with τ nk being the total scattering lifetime and the inverse τ −1 nk is the scattering rate given by: A common approximation that we refer to as the selfenergy relaxation time approximation (SERTA) consists in neglecting the second term on the right-hand side of Eq. (3). The mobility then takes the simpler form 14 : In order to analyze the role of band curvature in mobility calculations, we introduce a "mobility effective mass" as follows: This definition is obtained by setting τ nk to a constant value τ in Eq. (5), and by requiring that the resulting mobility can be expressed in the standard Drude form with the same relaxation time τ : Lastly, we introduce an approximate version of the SERTA mobility, which is useful to disentangle the convergence of electron band structures from that of the phonon energies and electron-phonon matrix elements.
To this aim we define the "constant adiabatic relaxation time approximation" (CARTA) as the mobility obtained by taking the SERTA of Eq. (5) with a constant electronphonon matrix elements and vanishing phonon frequencies: |g mnν (k, q)| = g 2 , with g an arbitrary constant with units of energy, and ω qν = 0:

B. Hall mobility
Calculation of the direct current Hall coefficient has seen a resurgence of interest 17,[49][50][51] . Experimentally, Hall mobility measurements are more common than time-offlight measurements of drift mobilities due to their superior accuracy and simplicity. This is reflected in the literature by a significantly higher number of available Hall mobility data compared to drift mobilities. However Hall measurements are performed in the presence of an external finite magnetic field, which introduces an additional Lorentz force on the carriers, thereby altering the mobility. To compare with experiment, we have to augment the BTE of Eq. (3) to account for the additional magnetic field B 6,17 : with τ nk being the total scattering lifetime defined in Eq. (4). Taking derivatives on both sides of Eq. (9) with respect to the Cartesian components of the magnetic field, B γ , at zero field, yields an equation for the linear response coefficients ∂ E β f nk and ∂ E β ∂ Bγ f nk : The Hall conductivity tensor is obtained from the second derivatives of the current density with respect to electric and magnetic field for vanishing fields: or equivalently directly from Eq. (9) as: Besides the drift and Hall conductivity and their mobility analogues, a commonly reported quantity is the dimensionless Hall tensor, which is defined as the ratio between the Hall conductivity and the drift conductivity 52,53 : where Eq. (13) with Here, x = ε/(k B T ) and we introduced the distribution function of the total decay rate: The numerical solution of the BTE requires the evaluation of a large number of electron-phonon matrix elements in order to converge the double momentum integrals (k and q). Various schemes have been developed to deal with this challenge including models with parameters computed from first principles 57 , direct evaluation of the electron-phonon matrix elements using densityfunctional perturbation theory (DFPT) 10 , linear interpolation of the scattering rates 11,58 , local orbital implementations 59 , smoothened Fourier interpolation 60,61 , Fourier interpolation of the perturbed potential 16,22,62 and interpolation based on MLWFs 63,64 . The use of MLWFs makes the calculations affordable and is the method of choice in this work, where we rely on the EPW software 63,65 . Such interpolation implies subtleties that need to be dealt with carefully, as discussed in Section II C.
C. Interpolation of the electron bands, phonon dispersions, and electron-phonon matrix elements The interpolation of the various quantities required to compute the drift mobility on ultra-dense grids relies on a discrete Fourier transform followed by a Fourier interpolation at arbitrary crystal momenta. Due to the fact that the discrete Fourier transform is in practice performed on a uniform finite grid, the quality of the interpolation depends on the localization in real space of each quantity.
In the case of the electronic Hamiltonian, we can leverage the phase freedom of the Bloch orbitals to create Wannier functions that are maximally localized in real space 48 . Since the interatomic force constants do not have such gauge freedom, we use direct Fourier interpolation. This interpolation requires some care in the case of polar materials, in order to correctly describe long-range interactions and the splitting of longitudinal-optical (LO) and transvers-optical (TO) modes 66,67 . Electronic Wannier and Bloch states are related by 48 : where R p is a lattice vector, N p is the number of unit cells in the Born-von Kárman supercell, corresponding to the number of k-points, and U nm,k is a unitary rotation matrix that transforms the Bloch wavefunctions to a Wannier gauge: In particular, the unitary transformation U can be chosen to obtain MLWFs that minimize the spatial spread functional 48 . In this context, the Hamiltonian H and dynamical matrices D can be transformed from the Bloch to Wannier representation as: where N p , N p are the number of unit cells in the Bornvon Kárman supercells for the electrons and phonons, respectively and e κα,qν are the eigendisplacement vectors corresponding to the atom κ in the Cartesian direction α for a collective phonon mode ν of momentum q. The extension of these two concepts to interpolate the electron-phonon matrix elements has been derived in Ref. 65 and yields: In practice, we compute the electron-phonon matrix elements g nmν (k, q) on a coarse momentum grid and rely on crystal symmetries to reduce the number of calculations to be performed. To carry out the Wannier transformation from coarse reciprocal space to real space, we first evaluate the matrix elements on the entire Brillouin zone using symmetries. As pointed out in a recent preprint 64 , in the case of SOC the spinor wavefunctions need to be rotated using the SU(2) unitary group. There was a bug in version 5.2 of the EPW software 63 related to this rotation. The bug has been eliminated in version 5.3, and for the systems tested so far we found a negligible impact. For example, we recalculated the mobilities of silicon, GaAs, and CsPbI 3 , and found differences smaller than 0.4% in all cases.
There are additional subtleties related to the Wannier interpolation of the various quantities required to compute the mobility. These aspects will be discussed in Section IV.

D. Band velocities
A key ingredient required in the calculation of carrier mobility is the carrier velocity as can be seen in Eq. (1). The velocity operator can be expressed in terms of the commutator between the Hamiltonian of the system and the position operator,v = (i/ )[Ĥ,r]. The matrix elements ofv in the presence of non-local pseudopotentials are [68][69][70][71] : wherep = −i ∂/∂k is the momentum operator, m e is the electron mass, andV NL is the non-local part of the pseudopotential. Neglecting the term withV NL in the case of a non-local pseudopotential is what we call the "local approximation" 14 . Within this approximation, Eq. (24) reads: (25) where u nk are the periodic part of the Bloch wavefunction. If in the second part of Eq. (25) the terms corresponding to the periodic part of the Bloch function are expanded into planewaves, u nk = G c nk (G)e iG·r , we obtain the simple velocity expression in the local approximation: We note that this expression is only valid in the case of local norm-conserving pseudopotentials. In practice, we will show that this local approximation is too crude and one needs to compute the matrix elements of the commutator in Eq. (24). We carry out this procedure in the Wannier representation, as described in Refs. 72-75. In particular, we construct the velocities interpolated at an arbitrary momentum k following Eq. (31) of Ref. 73: where A nmk is the Berry connection defined as: We construct the Berry connection in the Wannier basis [Eq. (20)] on the coarse k-point grid and compute the derivative using finite differences as 73 : where b are the vectors connecting k to its nearest neighbors and w b a weight associated with each shell of neighbours |b|. These vectors are constructed in such a way as to satisfy 76 using the smallest number of neighbours. The rotated overlap matrices (from Bloch to Wannier basis) are obtained as: where M nm (k, b) = u nk |u mk+b is the overlap matrix between the cell-periodic Bloch eigenstates at neighboring points k and k+b.  In each case, only one effect is considered, while all other settings are the same as for the reference calculation. We consider the following effects: (i) Long range part of the dynamical matrix (D) and electron-phonon matrix elements (G). The label "0" means that no long range part has been considered, and "DD", "DQ", and "QQ" mean that dipole-dipole, dipole-quadrupole and quadrupole-quadrupole contributions have been included. The labels "eD" and "eQ" mean that monopole-dipole and monopole-quadrupole interactions are included, respectively. (ii) Using the velocity computed in the local approximation via Eq.
where k belongs to the same grid employed for the Wannier interpolation of the band structure. The quantity on the left hand side decreases rapidly with R owing to the localization of MLWFs, therefore we can use these real-space quantities to interpolate back to arbitrary momenta k : We do the same for the k derivative of the Hamiltonian following Eq. (38) of Ref. 73: Finally we un-rotate from the Wannier basis to the Bloch state basis using the diagonalizers of the Hamiltonian: which gives the desired Eq. (27). We have implemented the interpolation of the velocity matrix elements including the correction for non-local potentials in the EPW software 63 . Figure 1 shows that the local approximation yields significant errors in the calculated mobilities, ranging from −80% to +40% of the correct value. For this reason, we use the exact velocities of Eq. (27) throughout this study.
To evaluate the quality of the Wannier interpolation of velocities, we perform direct DFT calculations using central finite differences, with six neighbouring points spaced by 5 · 10 −4 2π/a.

E. Long range corrections: dipoles and quadrupoles
The analytic properties of the long wavelength limit of the electron-phonon matrix elements in bulk polar materials, resulting from the Fröhlich interaction, have been know for a long time 77 , but the generalization to first-principles calculations using Wannier-Fourier interpolation 65 appeared only recently 78,79 . Following a line of thinking closely related to the Fourier interpolation of the phonon frequencies 66,67,80 , one can decompose the matrix elements into a short-(S) and long-range (L) contribution 78,79 : In Eq. (37), the long-range part is given by the infinite multipole expansion: where the first term comes from dipole contribution, then quadrupole, octopoles and higher. The analytic form of the long-range part of the dipole contribution is given as 78 : where G is a reciprocal lattice vector, τ κ is the position of atom κ, Z * κ is the Born effective charge tensor of the atom κ of mass M κ , e κqν is the vibrational eigendisplacement vector normalized in the unit cell, ε 0 is the vacuum dielectric constant and ε ∞ is the high-frequency dielectric tensor of the material. In the past, a Gaussian filter has been added to Eq. (39). Since this filter breaks the periodicity of the matrix element in reciprocal space, we recommend not to include it in future calculations. One can notice that Eq. (39) diverges as 1/|q| as we approach the zone center. This is an integrable singularity, so that upon performing integrations in reciprocal space one obtains finite values of physical properties. The quadrupolar component of the matrix element at long wavelength is given by: whereQ mnκ (k, q) is a third-rank effective quadrupole tensor defined as: and Q καβγ is the dynamic quadrupole tensor 81 . The third-rank quadrupole tensor should obey the same symmetry rules as other third-rank tensors such as the piezoelectric tensor or the second-order magnetoelectric tensor 82 . One can therefore use tools such as the Bilbao Crystallographic Server 83 to determine the non-zero components of the quadrupole tensor for a given space group, and their relations. In the case of all the polar materials investigated in this work, the quadrupole tensor takes the form: where ε αβγ is the Levi-Civita symbol and Q κ is an atomdependent scalar. In the case of silicon and diamond, the effective tensor in Eq. (41) is completely specified by a single scalar 81 : The quantity V Hxc,Eγ in Eq. (41) is the change of the Hartree and exchange-correlation potential with respect to the electric field E γ 22 . This term is null in the case of materials with vanishing Born effective charges and has been shown to be small in the case of polar materials such as GaAs and GaP 16,22 . Finally, we note that in the q + G → 0 limit, the wavefunction overlap in Eqs. (39), (41) can be written in terms of the Wannier rotation matrices U nmk as 78 : A similar approach can be used to enforce the correct behavior of the dynamical matrices at long wavelength 66,84 : where the nonanalytical, direction-dependent term, is given by the contribution of dipole-dipole, dipolequadrupole and quadrupole-quadrupole terms as: with In this expression we neglected octopoles and higherorders as well as fourth-order and higher dielectric functions. The sum over the lattice of charges is performed using the Ewald technique, with a parameter Λ. In our calculations we use Λ = 1 bohr −1 so that the real-space contribution in the Ewald summation can be neglected 66 .

F. Adaptive broadening
The numerical evaluation of Eqs. (3) and (4) requires to approximate Dirac delta functions by Lorentzian or Gaussian functions of finite broadening. This means that the resulting mobility will depend on the value of the broadening used, and careful convergence tests are required 14 . To circumvent this difficulty, various schemes have been proposed including adaptive broadening 85 and linear tetrahedron methods 22,86 .
We implemented the adaptive broadening approach in the EPW software following Eq. (18) from Ref. 85, where the state-dependent broadening is given by: where G α are the three primitive reciprocal lattice vectors, α the three crystalline directions, v nnk+q the electronic velocities defined by Eq. (27), N α the q-point grid density in the three crystalline directions. Eq. (48) can be scaled arbitrarily, and for this work we choose empirically a factor 1/2 after a number of numerical tests. In addition, the phonon velocity v qνν is obtained from the momentum derivative of the dynamical matrix, following an approach similar to Eq. (38) of Ref. 73: (49) In the case of polar materials, the momentum derivative of the long-range dipole part is obtained by taking the derivative of Eq. (46) as where we have considered only the dipole-dipole interaction term for simplicity. This choice does not change any of the results presented in this work, because it only affects the broadening parameter. This point will be discussed further in Section III F. All results described below were obtained using this adaptive broadening scheme.

III. COMPUTATIONAL METHODOLOGY
In this section we describe the software and computational parameters employed. We discuss the interpolation of electron-phonon matrix elements using the longrange dipole and quadrupole corrections. We look at the convergence rate with respect to the coarse and fine grid interpolation, and show that the drift mobility, Hall mobility, and Hall factor can be extrapolated to the limit of continuous Brillouin zone sampling. Finally, we propose in Appendix A a way to construct an optimal Wigner-Seitz cell for the double three dimensional Fourier interpolation of the electron-phonon matrix elements.

A. Density functional theory and density functional perturbation theory calculations
We use the Quantum ESPRESSO software distribution 87,88 with the relativistic Perdew-Burke-Ernzerhof (PBE) parametrization 89 of the generalized gradient approximation (GGA) to density-functional theory. All the pseudopotentials are norm-conserving, generated using the ONCVPSP code 90 , and optimized via the Pseu-doDojo initiative 91 . To assess the effect of the exchangecorrelation functional, we also regenerated the pseudopotentials using the ONCVPSP code, and changed the exchange-correlation functional from PBE to the Perdew-Zunger parametrization 92 of the local density approximation (LDA) 93 but keeping all other parameters unchanged. To assess the effect of pseudization parameters, we also used PBE pseudopotentials from the SG15 ONCV library 94 that were extended to fully relativistic pseudopotentials 95 . The electron wave functions are expanded in a planewave basis set with the kinetic energy cutoff reported in Table I, such that the total energy is converged to less than 1 mRy per atom. The Brillouin zone is sampled with a homogeneous Γ-centered grid of at least 20×20×20 k-points to converge linear-response quantities such as dielectric tensor, Born effective charges and phonon frequencies using DFPT 66,67 . The relative threshold for the self-consistent solution of the Sternheimer equations was set to 10 −17 or lower to ensure accurate first-order perturbed wavefunctions. All the DFPT data were produced in binary or XML format (as opposed to the default text format) to preserve machine precision. We include SOC in all our calculations. To obtain accurate results for properties such as the spin-orbit splitting, semicore states have been used for the pseudopotential of Ga, As, and Sb.
For definiteness we used experimental lattice parameters in all cases, as reported in Table I. The relaxed lattice parameters are also given in Table I, and tend to overestimate experiment by an average of 1.3%, which is typical for PBE calculations. As shown in Fig. 1, the lattice parameter can have a significant impact on the calculated mobility, leading to differences of up to 20% when comparing calculations with the optimized or the experimental parameter.
To ensure reproducibility and follow the FAIR principles (findable, accessible, interoperable, and reusable) 106   We indicate the initial projections and the number of iterations to reach a relative convergence of 10 −12 in the spread, the total spread, and the average spread per Wannier function for the converged kpoint grids. For the initial projection, a number next to the atom means that the initial projection is applied to only one of the two equivalent atoms in the primitive cell. The disentanglement window as well as the frozen window are given with respect to the band edge. All calculations include SOC.
When performing calculations of carrier mobility, and especially in the presence of magnetic fields, it is crucial to have symmetric and highly localized Wannier functions. In order to provide a benchmark for future work, here we describe the procedure followed to generate ML-WFs and the properties of the resulting matrix elements.
For diamond we choose 4 sp 3 orbitals (8 Wannier functions) located on each atom of the inversion-symmetric unit cell as initial projection for the valence and conduction manifold. This choice leads to a Wannier function spread σ 2 of 2.40Å 2 per function for the converged coarse k-point grid, and 0.79Å 2 for the valence bands, see Ta-   ble II. The strong localization of these MLWFs leads to a fast spatial decay of the matrix elements in the Wannier representation of the electronic Hamiltonian, the electron velocity, the interatomic force constants, and the electron-phonon matrix elements, down to nine orders of magnitude, as seen in Fig. 2. In all the tests reported here, we find as expected that Wannier functions describing the valence bands are more localized than for the conduction bands.
In the case of silicon, the Wannierization of the conduction band manifold is especially challenging. Using the wannier90 software we need 1320 iterations to converge to a relative accuracy of 10 −12Å2 in the spread (we note in passing the usefulness of the conjugate gradients algorithm such that we changed the default of resetting the search direction to steepest descents every 100 iterations, rather than 5, as default). For the conduction bands, we find that adding one d and one s orbital on one of the silicon atoms as initial projection works best, and yields a spread of 5.67Å 2 . The comparison between the electronic bandstructure of the conduction bands of silicon using this projection and an sp 3 one is shown in Appendix, Fig. 22. The challenge in obtaining carefully converged effective masses lies in the fact that the conduction band minimum of silicon does not lie at a high symmetry point and is therefore not included in the kpoint grid used for interpolation. The calculation of the conduction band manifold is accelerated in silicon and diamond by excluding the valence bands. The correspond-ing Wannier functions are shown in Fig. 2. They display a relatively complex shape that deviates from the simple chemical picture found for the other materials. The difficulty in generating Wannier functions for the conduction bands of silicon has already been discussed in earlier work 107 . For the conduction band of SrO, we find that using a similar combination of d and s orbitals on the Sr atom works better than sp 3 , but it also results in a more complicated and entangled set of Wannier functions, as shown in Fig. 2.
For GaAs we used 16 Wannier functions (8 Wannier functions times two due to spin-orbit coupling); on the Ga and As atoms, with sp 3 character. These span both the valence and conduction manifolds, and are used to calculate hole mobilities. To reduce the computational cost for the electron mobility, we used only 8 Wannier functions centered on the Ga atom, with sp 3 character, to describe the conduction manifold. Since we used semicore states for the pseudopotentials, we excluded the first 20 bands from the Wannierization in the case of valence bands, and 28 bands in the case of conduction bands. The spread of the maximally localized Wannier functions associated with Ga is 8.85Å 2 . The case of the valence band of GaAs is unique among the studied compounds. We find that Wannierizing both the valence and conduction bands yielded a smaller spread of 3.04Å 2 on the valence band manifold than using 6 or 8 Wannier functions of p or sp 3 character on the As atom with a spread of 4.12Å 2 or 3.29Å 2 , respectively. We attribute the decrease in spread with increasing number of Wannier functions in the case of GaAs to an increase in the flexibility of the Wannier function basis set and some hybridization with the conduction manifold due to the small bandgap in DFT. For SiC we used sp 3 orbitals located on the Si atom, with 8 Wannier functions for both valence and conduction, and in the conduction band calculations we excluded 8 bands. After optimization, the spreads are 4.36Å 2 and 1.21Å 2 for all orbitals, respectively.
For the other six materials, we report in Table II all the chosen initial projections, disentanglement energy windows, and frozen windows 108 . For completeness, we mention in Table III the bands that were excluded from the Wannierization. In crystals with inversion symmetry such as diamond or silicon, states with opposite spin orientations are degenerate. In zinc-blend structures however, the lack of inversion symmetry causes a spin splitting along all directions except [100] 109 . We verify in all cases that our choice of initial Wannier projections preserves the expected crystal symmetries in the spread and in the electronic band structure.
As seen in Fig. 2, the electronic Hamiltonian, velocity matrix, dynamical matrix and electron-phonon vertex in the limiting case of R = 0 or R = 0 in the Wannier representation as a function of R = |R − R | decay very rapidly. We know that the decay of the Hamiltonian for the valence band manifold must be exponential 110,111 . For the conduction band the exponential decay is not guaranteed, and for the dynamical matrix and  the phonon part of the electron-phonon matrix elements we expect a power law decay. To characterize the decay using a simple unified descriptor, we fit our data with an exponential function for all cases. The resulting decay length is reported in Table III for all materials.
In the specific case of the valence band manifold of silicon, the Hamiltonian decay length is 1.710Å, improving on earlier work which reported a decay length of 3.2Å with an empirical pseudopotential starting from four bond-centered trial functions 110 . The worst localization is found for the conduction band manifold of GaAs, GaP, and AlSb (> 2.7Å), while the best localization is achieved for the valence manifold of c-BN, diamond, and SrO (∼ 1Å). We also note that the spread and number of iterations required to reach convergence systematically increase with k-point grid density. Furthermore, it has been reported that convergence fails for ultra-dense grids (80×80×80 and above) 113 .
We also tested two aditional procedures for constructing Wannier functions, namely the selectively-localized Wannier functions (SLWFs) 114 and the symmetryadapted Wannier functions (SAWFs) 115 . With SLWFs one constrains the Wannier center r n to be located at the position r 0n using a Lagrange multiplier λ c as λ c n (r n −r 0n ) 2 , where the summation can be restricted to selected Wannier functions. In this case we consider the entire manifold. The multiplier λ c can be chosen arbitrarily; the larger the value, the more difficult the convergence, so that in practice the Wannier centers are close to but not exactly positioned at the target positions. The SAWFs guarantee that the Wannier functions respect the point group symmetry of the crystal. We note that the existing implementation in the wannier90 75 code does not currently support band disentanglement with a frozen window nor SOC. In both cases, we find that the minimum spread is larger than those reported in Table II. Therefore, we did not purse these avenues further.
The MLWFs obtained as described in this section were used to calculate interpolated band structures, phonon dispersions, and electron-phonon matrix elements in the following sections.

C. Dynamical quadrupoles
For practical calculations, we need a procedure for evaluating the quadrupole tensor introduced in Sec. II E. An efficient way to compute dynamical quadrupoles rests on perturbation theory and has been implemented in the Abinit software 81 . However, this implementation is currently limited to the LDA exchange and correlation functional, and to pseudopotentials without non-linear core corrections. To circumvent these limitations, we propose an alternative strategy for estimating the quadrupole tensor, which can be used with any exchange and correlation functional and pseudopotential type. The general idea is to determine the quadrupole tensor by matching Eqs. (38)- (41) to explicit DFPT calculations at small q.  The minimum number of calculations required is equal to the number of independent components of the tensor. In practice, to minimize numerical noise we compute the DFPT matrix elements for several points on a sphere with |q| = constant. Comparison between the deformation potential of the Γ1 valence band of silicon from the direct DFPT (dots) and the Wannier interpolation with and without quadrupoles, for different k-point grids (the q-point grid is half of that). The components included in the dynamical matrix and in the electron-phonon matrix elements are labelled as "D" and "G", respectively. The labels "DD", "DQ", and "QQ" indicate the dipole-dipole, dipole-quadrupole and quadrupole-quadrupole contributions, while "eD" and "eQ" denote the monopoledipole and monopole-quadrupole interactions, respectively. To this aim, we use a Python script with the long-range analytic electron-phonon matrix elements from Eq. (38), and the phonon frequencies, Born effective charges, dielectric constant, eigendisplacements and wavefunction overlaps from the DFPT calculations. We then perform least squares optimization with respect to the components of the quandrupole tensor, and include a modedependent offset as an additional minimization parameter to capture the O(|q| 0 ) term of the matrix elements, following Eq. (3.17) of Ref. 77. For the compounds considered in this work, we constrain the quadrupole tensor to have the form of Eq. (42) or Eq. (43) during the minimization, although the procedure holds for tensors of any symmetry.
We emphasize that this procedure for calculating the quadrupole tensor is inexpensive, as it requires the calculation of only a dozen q-points, see Fig. 23.
As a sanity check, we computed the quadrupoles of 3C-SiC using the Abinit software 62,116 . To use the Abinit implementation we considered LDA pseudopotentials without non-linear core correction and without SOC. We obtained Q Si = 7.025 e bohr and Q C = −1.416 e bohr. Using the same pseudopotentials and the Quantum Espresso software 88 we obtainedQ Si = 7.392 e bohr andQ C = −2.483 e bohr. We attribute the small difference to the fact that the implementation of Ref. 22 does not include the second term of Eq. (41), and that our procedure does not capture the q → 0 limit with the same precision as in DFPT.
Using the above procedure, we computed the effective quadrupole tensor for all ten compounds. The results are reported in Table IV. Among the compounds considered here, GaAs has the quadrupole tensor with the largest elements, while this tensor vanishes in SrO due to the cubic F m3m space group symmetry. In Table IV we also report the calculated quadrupole terms from Refs. 16,22 for comparison, as obtained using perturbation theory, LDA without non-linear core correction, and neglecting SOC. Our calculations compare well with Refs. 16,22, although close agreement is not expected due to the aforementioned differences.

D. Interpolation of the electron-phonon matrix elements
To assess the quality of the Wannier interpolation, we compare the interpolated electron-phonon matrix element with those obtained from a direct DFPT calculation. For an easier comparison, we compute the total deformation potential 79,117 : where the k = Γ point was chosen, the sum over bands is carried over the N w states of the Wannier manifold, and ρ is mass density of the crystal. The advantage of examining the deformation potential rather than the matrix elements directly is that electronic degeneracies are traced out, and the frequency that enters the zero-point amplitude is removed, so that one can concentrate on the linear variation of the potential. To investigate specific states, we also define the band-resolved deformation potential D ν mn (Γ, q) where the band summation in Eq. (51) has been removed.
In Fig. 3, we compare the deformation potential of silicon for the lowest valence band state at Γ, for the direct DFPT calculation (black dots) and the Wannier interpolation, with and without quadrupoles (lines).
In silicon there is no Fröhlich contribution to the matrix elements due to inversion symmetry. There are, however, quadrupole contributions. As can be seen in Fig. 3, the Wannier interpolation without quadrupole correction reproduces the direct DFPT calculation everywhere, except close to Γ in the ΓL and ΓK directions for the LO and TO modes. The discontinuity of the LO and TO modes at Γ poses a challenge to the interpolation without quadrupole correction, as shown by the solid lines in gray. Upon including quadrupole corrections, the deformation potential is found to be independent of grid size. We therefore only report the curve corresponding to the grid used throughout this manuscript (see Sec. III E). In contrast, the longitudinal acoustic (LA) mode of silicon shown in Fig. 3 does not exhibit any discontinuities at the zone center, and is therefore well described without quadrupole correction. The deformation potential of the transverse-acoustic modes (TA) vanishes and is not shown. These data are consistent with Ref. 15. Figure 4 shows a systematic comparison between the deformation potential computed via explicit DFPT calculations and the results of our Wannier interpolation including dipole, dipole-quadrupole and quadrupole corrections using Eq. (51). We note the sharp change in deformation potential occurring near the K point in diamond, and to a smaller extent in c-BN. This effect is due to the avoided crossing in the phonon bandstructure of diamond and c-BN at the same point in momentum space. Quantitatively, diamond and c-BN show the largest deformation potential, which is consistent with these compounds being the hardest materials in our set.
Overall, the agreement between the direct calculation and the Wannier interpolation with quadrupole correction is excellent, guaranteeing a good interpolation. The only exception is SrO which is the only material investigated here for which the quadrupole contribution vanishes by symmetry. This finding suggests that a very accurate interpolation for SrO might require the inclusion of octopole contributions. We leave the investigation of the effect of octopoles for further studies.

E. Convergece of the mobility with the coarse Brillouin zone grid
We present in Fig. 5 the convergence rate of the electron and hole carrier mobility with increasing coarse grid sizes. Overall, we find that using the same coarse k and q grid leads to the fastest convergence. On the other hand, using a coarse k-grid three times denser than the q-grid leads to similar results but at a much higher computational cost, as shown in Fig. 24 in the Appendix for c-BN. Therefore it is recommended that calculations be performed using the same kand qcoarse grids. We note that the calculations presented below were performed using a k-grid twice as dense as the q-grid for historical reasons, but in future work equal grids will be employed.
To perform the convergence study, we use either a 60 3 or a 100 3 fine grids as in indicated in Fig. 5. We verify that the cross convergence between coarse and fine grids is weak. As can be seen in Fig. 5, the interpolated mobility converges slowly with respect to the coarse grid size, despite dipole and quadrupole corrections being already included in the calculations. The origin of this slow convergence has not been investigated so far. To shed light on this effect, we compute the mobility effective mass defined in Eq. (6) and report the convergence of this purely electronic quantity with orange dots in Fig. 6. We note that the mobility effective mass reflects the quality of the Wannier interpolation of the band structures. Its value should converge to the mobility effective mass computed directly from DFT bands without interpolation. In Fig. 6 the reference DFT value is shown as horizontal dashed orange lines. Figure 6 shows that the mobility effective mass plays a significant role in the slow convergence of the mobility as a function of coarse grid size. However, the convergence of the mobility effective mass is a necessary but not sufficient condition for the overall convergence of the mobility with coarse grid size. Since computing the mobility effective mass is much cheaper than computing the mobility, this quantity can be used to estimate a minimum coarse grid size. In the cases of diamond, Si, 3C-SiC, AlP, GaP, AlAs, AlSb and SrO in Fig. 6, we see that the electron mobility effective mass does not behave consistently with the convergence of the mobility. It is also the case for the hole mobility effective mass of GaAs and AlSb. The reason for this deviation is that the effective mass does not fully capture the scattering phase space, since it neglects the energy conservation connected with the electron-phonon scattering amplitudes.
To incorporate these effects, we employ the CARTA approximation defined in Eq. (8). In this case, we have to compute the sum over q-points explicitly, making it more computationally involved, but still far cheaper than the SERTA since we do not have to interpolate the dynamical matrix and electron-phonon matrix elements. For the cases considered here the CARTA calculations are at least two orders of magnitude faster than the SERTA. We note that in CARTA calculations we cannot employ adaptive broadening, therefore we use a constant 10 meV Gaussian broadening for the Dirac deltas. The CARTA results are shown as magenta dots in Fig. 6, and show how the electronic degrees of freedom influence the convergence with respect to the coarse Brillouin zone sampling. Since these calculations do not carry information about the electron-phonon matrix element, we present the data as a deviation from the CARTA mobility computed without Wannier interpolation, by using directly DFT eigenvalues. By comparing Figs. 5 and 6, we observe that the CARTA mobility converges at a similar rate as the BTE mobility. This observation suggests to speed up the coarse grid convergence by factoring out the electronic degrees of freedom. To this aim, we rescale our BTE mobilities µ αβ using the following expression: where µ CARTA,DFT αβ is the mobility obtained in the CARTA with DFT eigenvalues and velocities. Eq. (52) allows one to correct for the small off-set observed in some materials for the effective mass mobility. Moreover, computing µ CARTA αβ is over two orders of magnitude faster than µ αβ since we only need to interpolate the eigenvalues and velocities. The results for the scaled BTE mobility are shown as orange dots in Fig. 5. Importantly, not only µ scaled αβ converges faster, but it also does so in a smoother and more systematic way.
The vertical thin gray line in Fig. 5 is the coarse grid that we consider converged and that is used in all subsequent calculations, unless stated otherwise. For the electron mobility of GaAs we do not report the mobility effective mass since we used the experimental value to describe the bands and velocities. From Fig. 5, we can note that the hole mobility converges in a smooth and systematic way while this is not necessarily the case for the electron mobility. This behavior reflects the better real-space localization of Wannier functions obtained for the valence manifold than the conduction manifold, see Fig. 2.
For the case of silicon, the electron mobility is significantly improved upon using the scaling of Eq. (52). However, we can notice that for our chosen converged grid, the calculated mobility does not converge to the correct limit without scaling (see Fig. 5). Specifically, the interpolated CARTA result overestimates the CARTA/DFT result by 1 %. We expect a similar error to occur for the BTE mobility. We anticipate that the slight difference will disappear upon using denser grids or by using a more extended Wannier functions basis set, but we have not investigated this further due to the computational cost. An overestimation of 2 % is also observed for the hole mobility of AlSb. In all other cases, the interpolated CARTA values converge smoothly to the CARTA/DFT result.
We are now in a position to answer the question on why is the convergence of the carrier mobility slow with re- spect to coarse grid size even when dipole and quadrupole corrections are included. The answer is material dependent. In the case of the electron mobility of Si, 3C-SiC, AlP, GaP, AlAs, AlSb, and SrO, the slow convergence is mostly due to purely electronic effects as demonstrated by the CARTA results. Specifically, the Wannierinterpolated electron bands converge slowly with grid size. A possible way to accelerate convergence would be to increase the number of Wannier functions in order to improve the flexibility of the basis functions. In general, for the valence bands the convergence is smoother and more systematic, and, therefore, less can be gained by using the CARTA rescaling. Nevertheless, we noticed some improvements for diamond, AlP, GaP, c-BN, AlAs, and AlSb.
In the case of diamond and c-BN, very little improvement is noticed using the CARTA rescaling. This indicates that for these two materials the slow convergence can be attributed to a slow convergence of the phonon dispersions and electron-phonon matrix elements with grid size.
In Fig. 7 we compare our coarse grid convergence results with previous work for Si and GaP in Ref. 16. Since only the SERTA mobility is reported in Ref. 16, we compare with our results at the same approximation level. Overall our results agree quite well given that there are a number of differences in the calculations. The coarse k-point grid in Ref. 16 was fixed to 18 3 , while we used a coarse k-point grid that is twice the one reported in Fig. 7 for the q-point grid. We used an equal fine kand q-point grids of 150 3 for Si and 120 3 for GaP, while Ref. 16 used a fine 72 3 k-point grid and a 144 3 q-point grid for Si and a 78 3 k-point grid and a 156 3 q-point grid for GaP. Different pseudopotentials and exchange correlation functionals were used as well. In addition, we included SOC, while Ref. 16 neglects it, but as shown in Fig. 1, the effect on electron mobility is negligible. The main difference between our results and those of Ref. 16 is that our calculations including dipole but not quadrupole corrections seem to converge rapidly, while in the earlier work there is a slow drift of in the SERTA mobility calculated with dipole correction. However, in agreement with Ref. 16, we observe that the inclusion of quadrupole corrections leads to a shift of the calculated mobility even at convergence.

F. Convergence of the mobility with the fine Brillouin zone grid
Using the converged coarse grids determined in the previous section, we can proceed to testing the convergence of the mobility with respect to the fine grid. The presence of the ∂f 0 nk /∂ε nk term in the mobility, see for example Eq. (5), implies that only those states close to the valence band or conduction band edge will contribute. To reduce the computational cost of these calculations we use the following scheme. We construct a fine k-point grid composed of only the points within a small energy window around the band edges, and rely on crystal symmetry to further reduce the number of points. For example in the case of the hole mobility of c-BN, using an energy window of 0.3 eV and an homogeneous fine grid of 250 3 points, we only have 12,390 irreducible k-points. Importantly, at no point we store the entire uniform grid in memory. This reduction allowed us to use extremely dense grids for the electron mobility of GaAs, up to 800 3 points in the first Brillouin zone.
Since all q-points can contribute (for example through inter-valley scattering), we do not use symmetry reduction, and keep all the q-points of the uniform grid that lie with the small energy window. For example, in the case of c-BN and a grid of 250 3 points, we have to compute 814,981 q-points. To this aim we generate the q-points on the fly, so that we do not have to store the full grid in memory.
Finally, we interpolate the electron-phonon matrix elements for all these k and q points and store on disk only the kand q-points identified above that contribute to the mobility. To retain only the information that will contribute to the mobility calculation, for each q-point we store the matrix elements for the nk states that fulfil the following condition: where N k , N q and N b are the total number of k-points, qpoints, and Wannierized bands, respectively. The threshold in Eq. (53) guarantees that the cumulative absolute error made by neglecting very small matrix elements will be lower than machine precision, or 10 −16 . In the c-BN example mentioned above, this procedure required the storage of a binary file of approximately 40 GB instead of multiple TBs required to store all the matrix elements. Once all the scattering rates are written to disk, we read the file in parallel using the message passing interface input-output (MPI-IO) and solve the iterative BTE. The first step of the iterative solution yields the SERTA mobility. The iterative solution of the BTE is extremely fast (a few minutes) and converges typically within a few tens of iterations. For the compounds considered in this work, we determined that the optimal energy window is between 0.2 and 0.3 eV, as shown in Table VII of the Appendix.
The results for the SERTA and BTE mobility and Hall factor are shown for all materials in Fig. 8. For similar reasons as in Ref. 47, we stress that we are showing the results as a function of inverse grid density and the fine k-point and q-point grids are always the same. For example, the label 1/100 on the horizontal axis corresponds to a 100 3 fine grid, such that the left hand side of each figure corresponds to infinitely dense sampling. Interestingly, both the mobility and the Hall factor appear to converge linearly with inverse grid density and can therefore be extrapolated. The linear interpolation range is depicted in Fig. 8 with a gray shaded region. We perform the linear extrapolation as soon as we enter the linear regime. This linearity comes from the fact that we compute the momentum integrals using the basic rectangular integration, therefore the error scales linearly with the grid spacing.
We notice that the electron and hole mobility of diamond, Si, c-BN as well as the electron mobility of GaAs and 3C-SiC show a relatively steep slope, therefore it is important to extrapolate the values for accurate results. In contrast, all the other cases considered here exhibit a milder slope, therefore the extrapolation is not necessary in these cases.
The linear behavior of the mobility is closely linked to the smearing used for the energy conserving Dirac deltas in Eq. (3). In Fig. 9 we focus on the electron mobility of c-BN since the linear trend is very pronounced, and we show how the carrier mobility depends on the grid density for various Gaussian smearing parameters. We see that, the smaller the smearing, the more pronounced are the fluctuations, and finer grids are needed before the linear regime is reached. The extrapolated value depends on the smearing used. The above discussion does not take into account the finite lifetimes of carriers in the calculations. In reality, a finite smearing in the range of 10-50 meV is to be expected as a result of the electronic linewidths. This effect could be incorporated by setting the smearing to /τ nk with τ nk given by Eq. (4), but we have not explored this direction.
Since the linear extrapolation at fixed smearing is rather tedious, in practice we proceed by using the adaptive smearing introduced in Eq. (48). As seen in Fig. 9, this approach offers a very good approximation to the zero-smearing extrapolation. All data presented in Fig. 8 have been generated using adaptive smearing. We also show a similar analysis for the case of the electron mobility of GaAs in Fig. 25 of the Appendix.
In Fig. 8 we also present the convergence of the SERTA and BTE Hall factor defined in Eq. (13). In this case, the Hall factor systematically converges linearly with grid density and can be extrapolated to infinite grid density. Also in this case we can either extrapolate the results to zero smearing, or use adaptive smearing. The results in either case are extremely close, therefore we proceed with adaptive smearing in the following.

IV. RESULTS
In this section we discuss the main results of our calculations for all the compounds considered in our work. We first describe the electronic and vibrational properties, then we discuss the Hall factor and its temperature dependence, as well as drift and Hall carrier mobilities. We also analyze in details the microscopic mechanisms responsible for limiting the carrier mobilities. Lastly, we compare our results to available experimental data.

A. Wannier interpolation of electron bands and phonon dispersions
The interpolated electronic band structures and phonon dispersions for the ten compounds considered in this work are presented in Fig. 10.
The phonon dispersions are compared to neutron scattering data. Close agreement between theory and experiment is found in all cases except for SrO, where the theory underestimates measured frequencies. The associated Born effective charges, high-frequency dielectric constants and maximum phonon frequencies are reported in Table V. The average deviation for the highest phonon frequency between calculations and experiments is 1.6%.
Since DFT/PBE systematically underestimates the electronic bandgap, the dielectric constant is systematically overestimated with respect to experiment. For the systems considered here, the average overestimation of the dielectric constant is 10%, but the largest discrepancy is as large as 30%.
As expected, the band gaps are severely underestimated. The average deviation from experiments is of 41%. On the other hand, the spin-orbit splitting is in slightly better agreement with experiment, the average deviation being 18%. Given these data, we expect that the electron-phonon matrix elements will be overscreened, and we anticipate that our calculated mobilities will tend to overestimate experimental data. We will come back to this point in Sec. IV E.
The DFT mobility effective masses obtained without Wannier interpolation are reported in Table V. The val- ues range from 0.1 to 0.8 m e . Direct comparison of such mobility effective mass with experiment is not possible, however we can compare the effective masses calculated using the standard definition with experiments. We report all the experimental and theoretical electron and hole effective masses in Table V    For AlP there are no experimental effective masses available, but we note a k · p study in relatively good agreement with our results which reports 2.68 138  Overall, the agreement between DFT effective masses and experiment is reasonable (within a factor 2-3) but there are a few large discrepancies. For the electron effective masses, the transverse mass of GaP and AlSb are overestimated with respect to the experimental values by 63 % and 82 %, respectively. In the case of silicon, the computed electron effective masses are very close to the experimental values. We note that our results are in line with previous calculations employing perturbation theory 162 . Since DFT strongly underestimates the effective mass of GaAs, with values ranging from 0.03 m e 163 to 0.053 m e 22 , we decided to use the experimental effective mass for the calculation of the electron mobility. To this aim we employed a parabolic band approximation to evaluate analytically the eigenenergies and electron velocities from the experimental effective mass. This semi-empirical procedure was used only for the electron mobility of GaAs.
In the case of hole effective masses, the largest discrepancies are found in the cases of c-BN (the mass is underestimated by almost a factor of three), GaAs, Si, and AlAs (all underestimating experimental values by 40-60%). The impact of these discrepancies on the calculated mobilities will be analyzed in Sec. IV E.

B. Impact of various approximations on the drift mobility
Having introduced all the required concepts and results, we are now able to discuss Fig. 1. For each of the approximations tested, we only change a single parameter, keeping everything else the same as in our most accurate calculations. The reference calculations have been performed using converged coarse and fine grids, velocities including non-local contributions, dipole and quadrupole corrections, SOC, BTE, adaptive broadening, and experimental lattice parameters. The following considerations also apply to the Hall factor.
We discuss the impact of the various approximations in order of importance. We do not discuss the constant relaxation time approximation (CRTA) which depends on the empirical choice of the scattering rates and whose predictive power is unclear 57 .
The most severe approximation concerns the longrange treatment of the dynamical matrix and the electron-phonon matrix element. In Fig. 1 we break down the long-range treatment into four levels: (i) no longrange treatment; (ii) including Fröhlich dipoles; (iii) including dipole and quadrupole corrections in the dynamical matrix; and (iv) including dipole and quadrupole corrections in the electron-phonon matrix elements.
Completely neglecting the long-range behavior has been recognized in 2015 to be a severe limitation 78,79 but its impact on the mobility has not yet been quantified. The consequence of neglecting long-range contributions during Fourier transformation of the matrix elements is that the Fröhlich divergence when approaching the zone centre is spuriously suppressed. As a result, the overall electron-phonon coupling, and therefore the scattering rates, will be strongly reduced. This leads to a strong overestimation of the carrier mobility, up to 220% in the case of the electron mobility of SrO. Surprisingly, the hole mobility of 3C-SiC is only slightly underestimated (8%) when ignoring long-range couplings. This error cancellation is due to incorrect eigendisplacement vectors and will be discussed shortly.
Next in order of significance, with a standard deviation of 30%, is the neglect of dynamical quadrupoles during the interpolation of the dynamical matrix and the electron-phonon matrix element. Neglecting quadrupoles typically leads to an underestimation of the mobility by an average of -5.6%, and up to -45% in the case of 3C-SiC. A notable exception is the case of c-BN, for which the neglect of quadrupole corrections leads to an overestimation of the mobility by 71%. This important effect has been discovered recently, and to our knowledge it has been taken into account only in a handful of works 15,16,22,25 . In the light of these developments it will be necessary to revisit prior mobility calculations.
In the case of silicon, neglecting quadrupoles for the LO and TO modes in the Γ − L and Γ − K directions shown in Fig. 3 produces spurious oscillations, leading to a 4%/5% overestimation of the electron/hole mobility.  11: (a) Valence band deformation potential corresponding to c-BN from 0 to 2 eV/bohr. The label "D" refers to the dynamical matrix, the label "G" refers to electronphonon matrix elements. The labels "DD", "DQ" and "QQ" indicate that dipole-dipole, dipole-quadrupole, or quadrupolequadrupole contributions are included. The labels "eD" and "eQ" indicate monopole-dipole and monopole-quadrupole terms. The black circles are the reference DFPT calculations. (b) Same as in (a), but for the Γ1 valence band of 3C-SiC.
In the case of diamond, the neglect of quadrupole corrections yields an overestimation of less than 1% of the mobilities. Data for all compounds are collected in Fig. 1, and their magnitude is in line with prior work 22,164 . The effect of quadrupole correction can be substantial in some compounds. For example, as seen in Fig. 1, the hole mobility of c-BN is strongly overestimated. To better understand why this is the case, in Fig. 11(a) we compare the deformation potential obtained by including or neglecting the quadrupole term. We can see that the high-energy modes are unaffected by quadrupoles in this case. Conversely, the quadrupole correction becomes important for low-energy acoustic modes. The red lines in Fig. 11(a) shows the result of the interpolation performed by including only the dipole long-range correction. This result clearly deviates from the direct DFPT calculations, shown as the black dots close to the zone center. In particular, the deformation potential corresponding to the LA mode in the L − Γ and K − Γ directions is strongly underestimated, leading to the observed large overestimation of the hole mobility of c-BN (in Fig. 1).
Interestingly, we found that including quadrupole corrections to the electron-phonon matrix elements only was not sufficient to correctly describe the deformation potential, see the orange lines in Fig. 11. Indeed, also in this case the deformation potential of the LA mode is overestimated, albeit the overestimation is less significant, and the hole mobility of c-BN is underestimated accordingly (23%, see fourth column of Fig. 1). This result is counter-intuitive since the phonon frequencies are extremely well described by including dipole-only terms in the dynamical matrix, see Fig. 20 for example (in certain classes of materials such as piezoelectrics this is often not the case 84 ). A careful investigation shows that the root cause which explains the difference between the orange and green curves in Fig. 11 are the eigendisplacement vectors obtained by diagonalization of the dynamical matrix. The eigendisplacement vectors e κα,qν enter in Eq. (23) during the Fourier transformation of the electron-phonon matrix elements and within the longrange term in Eq. (40). It is therefore crucial to include at least dipole-dipole, dipole-quadrupole and quadrupolequadrupole in the dynamical matrix, see Eq. (46) as well as dipole-dipole and quadrupole-quadrupole effects in the electron-phonon matrix elements, see Eq. (40) to achieve a reliable interpolation of the deformation potential. We also emphasize that this correction must be applied both to the dynamical matrix and the scattering matrix elements: a partial correction of the dynamical matrix only is not sufficient, as shown by the purple line in Fig. 11 and the third column of Fig. 1.
Another interesting situation is found for 3C-SiC. We show the deformation potential for the lowest energy band in Fig. 11(b). Again we can see that the deformation potential close to Γ along the K − Γ direction has a spurious finite value at the zone center. Including quadrupoles in the dynamical matrix and electronphonon matrix elements fixes the problem. Among all materials, we note that including quadrupole correction in the electron-phonon matrix element but not in the dynamical matrix, yields an average underestimation of the mobility by 5%.
Beyond long-range dipole and quadrupole corrections, the most significant source of error in mobility calculations is the local velocity approximation, whereby the non-local part of the pseudopotential is neglected. This approximation yields a large standard deviation of 28% with values ranging from -71% in c-BN to +37% for 3C-SiC. We note that the local approximation was used by some of us for the calculation of the mobility of silicon 14 .
We are now able to state that this approximation yields an overestimation of the mobility by +17% and +7% for the electron and hole mobility of silicon, respectively. It is worth mentioning that the velocity including non-local contributions was implemented in the EPW software shortly after the publication of Ref. 14, so that subsequent work is not impacted by this approximation.
The next source of error in order of importance is the neglect of SOC. Without SOC the mobility is underestimated by up to 62%. For the compounds investigated here, this effect is almost negligible for electron mobility, while it is significant for hole mobility. The strongest underestimation occurs in AlSb with a -62% decrease in mobility. These results are in line with earlier findings 13,14 . The reason for this underestimation is that SOC only affects valence bands in the semiconductors investigated here, because the valence band maximum is predominantly of p character. The effect of SOC is to lower the spin-split band, thereby reducing the phase space available for scattering, which increases the mobility. For systems with p states at the conduction band bottom we expect a similar effect, for example in halide perovskites 46 . Less significant than the above but still important is the error that one makes by using the popular SERTA method. This approximation systematically underestimates the mobility within the range 1% to 40%. At an intuitive level, the fact that the SERTA mobility underestimates the BTE mobility can be understood as follows: in the SERTA all scattering processes reduce the electron current in the same way, while in the BTE one takes into account the fact that forward scattering does not reduce the current as strongly as backward scattering. Therefore the BTE is expected to yield a higher mobility.
The last three effects shown in Fig. 1 have a smaller impact on mobility, with errors up to 20%. We go over these effects in order of appearance in the figure. First, we investigated the sensitivity of the mobility to the lattice parameter. To this aim we re-calculated all mobilities using the relaxed lattice parameters. As discussed in Sec. III A, the PBE lattice parameter is overestimated with respect to experiment by an average of +1.3%. This effect yields an average decrease of mobility by -6.5%. We then considered the effect of the exchange-correlation functional by using LDA instead of PBE but keeping all other parameters the same. In this case we find only small changes, with an average of -1.5%. Lastly, we explored the effect of pseudization parameters. To this aim we repeated the calculations using the SG15 ONCV library 94 that was extended to fully relativistic potentials 95 instead of the PseudoDojo library 91 . In this case we find an average deviation of the order of 2%. The small impact of these effects on the calculated mobilities highlights the robustness of our approach.
We note that many other effects can impact the calculation of carrier mobilities, including electron-twophonon interactions 20 , thermal lattice expansion 14 , and band structure renormalization due to electron-phonon coupling 47,165 . We expect these effects to be small in most materials, but it will be important in the future to perform a detailed analysis. The effect of quasiparticle corrections 14,24 and many-body renormalization of the electron-phonon matrix elements 166,167 will be discussed in Sec. IV E in relation with experimental data.

C. Dominant scattering mechanisms
To gain a better understanding of the microscopic phenomena responsible for limiting the intrinsic transport properties of the compounds investigated here, we show in the Appendix Fig. 26 the electron and hole scattering rates τ nk computed via Eq. (4) at 150 K, 300 K, and 500 K. Without surprises, the higher the temperature, the higher the scattering rates and the lower the mobility.
As expected, in most cases we find a mild increase of the scattering rates as we move away from the band edges, which corresponds to acoustic-phonon scattering. At the onset for the emission or absorption of optical phonons, the rates show a marked increase. For example in diamond this sharp increase occurs at around 160 meV, when optical phonon emission/absorption processes become allowed. We note a couple of exceptions to this general trend: the scattering rates of electrons in GaAs and SrO, and both the electron and hole scattering rates  of c-BN first decrease when moving away from the band edges which is attributed to LO phonon absorption and was investigated in Ref. 19.
In Fig. 12 we present the mode decomposition of the room-temperature scattering rates, as well as the total scattering rates. Starting from the holes: For diamond, the onset of scattering is mostly related to transverse acoustic mode scattering, while for silicon the longitudinal acoustic mode also plays an important role. In the case of GaAs, AlP, GaP, AlAs, AlSb, and SrO, the hole scattering near the band edges is dominated by the longitudinal optical mode. For the electrons: The scattering of electrons near the band edge in GaAs, AlP, GaP, AlAs, AlSb, and SrO is dominated by longitudinal optical modes, while diamond, silicon, 3C-SiC and c-BN show significant contributions from LA and TA modes.
In 1966, Berman, Lax and Loudon 168 have derived intervalley-scattering selection rules for III-V semiconductors. For example, in the case of electrons in GaP, X-X intervalley scattering is expected. The rule is that, if the mass of the group-V element is heavier than that of the group-III element, then LO scattering should be allowed by selection rules, otherwise the LA scattering would be allowed. Therefore, for GaP one would expect LA-phonon scattering, while for AlP, AlSb, and AlAs one would expect LO-phonon scattering to dominate the intervalley scattering. Looking at Fig. 12 confirms the dominance of LO scattering in AlP, AlSb, and AlAs. For the case of GaP, it seems the scattering at an energy 3/2k B T away from the band edge has about equal contribution from LA and LO phonon modes.
A complementary viewpoint on the scattering rates is obtained by considering only carriers at a given energy, and resolving the rates by phonon energy. This analysis is shown in Fig. 13. The calculations correspond to carriers at an energy of 3k B T /2 away from the band edge and room temperature, following previous work 46 . From this analysis, we can see that the scattering rate is overwhelmingly dominated by acoustic scattering in diamond and c-BN. In contrast, optical scattering dominates for electrons in GaAs, SrO, AlAs, AlP, and for holes in GaAs, AlAs, AlSb, GaP, and AlP. We also mention a few counter examples to the common wisdom whereby if a material is dominated by deformation potential acoustic scattering, then the temperature exponent of the mobility is (T −3/2 ) 169 . Diamond, silicon, c-BN and 3C-SiC are all dominated by acousticphonon scattering, but the temperature exponents evaluated on the 150 K-500 K range are -1.81, -2.05, -1.33, and -2.86 for the electrons and -2.06, 2.93, -1.51, and -1.97 for the holes, respectively. These data suggest that one should exert some caution when applying simplified models based on the electron gas to real materials.
For completeness we also examine the predictive power of the classic Drude formula [Eq. (7)]. To this aim we compare our BTE calculations for all the compounds considered in this work with the Drude prediction. In the Drude formula we employ the mobility effective masses reported in Table V and the total scattering rates evaluated as in Fig. 13. The resulting data are shown in Fig. 14. The Drude formula reproduces the trend of ab initio BTE calculations very closely, although it tends to overestimate mobilities by a factor of 1.1 to 3.1. This favourable comparison implies that the materials descriptors that we introduced, namely the mobility effective mass and the scattering rates of carriers at 3k B T /2, provide a meaningful description of carrier transport in the systems considered here.

D. Temperature dependence of the Hall factor
In this section we investigate the temperature dependence of the Hall factor. We use Eqs. (9)-(13) to compute the Hall factors within both the SERTA and full BTE. For each temperature we extrapolate the fine grids as discussed in Sec. III F, and we present the results in Fig. 15. We find that the Hall factor can vary substantially with temperature, ranging from 0.7 to 1.9 at room temperature. These results suggest that caution should be used when comparing ab initio calculations of carrier mobilities to experiments, since most reported experimental data refer to Hall mobilities, not drift mobilities. As a technical note, we mention that below 150 K we could not converge the Brillouin zone grids due to excessive computational cost. Accordingly these data are not shown in Fig. 15.
Our calculations indicate that the Hall factor exhibits rather complex temperature trends, which reflects the variety of phonon energy scales involved.  Comparison with experimental data for the Hall factor is challenging given the scarcity of available data. We were only able to find experimental values for the hole Hall factor of diamond and for the electron Hall factor of GaAs, which we report in Fig. 15. In both cases we find reasonable agreement between our calculations and experiments, with the largest deviation being of the order of 30% (note the magnified scales in Fig. 15).
In a recent preprint 51 , the authors calculated the electron Hall factor of silicon to be 1.15 at room temperature. These results are consistent with our Fig. 8 for the same fine grid where we have an electron Hall factor of 1.12. The authors of Ref. 51 also computed the electron Hall factor of GaAs, finding a BTE value ranging from 1.22 to 1.14 in the 250-400 K temperature range. Their results compare well with our calculations, yielding a Hall factor of 1.28 to 1.13 in the same temperature range.
In Fig. 15 we also show the Hall factor as computed using the standard isotropic approximation given by Eqs. (15)-(16) (gray dashed lines). For these calculations we evaluated the Dirac deltas using a Gaussian of width 1 meV. This approximation reproduces the BTE results quite well for diamond, GaAs, AlP(e), GaP(h), c-BN(e), AlAs, AlSb, and SrO. However, it fails even qualitatively in the case of silicon, AlP(h), and GaP(e). One important limitation of the isotropic approximation is that it cannot lead to a Hall factor smaller than unity. In this section we examine the accuracy of the various methods presented here in reproducing experimental carrier mobility. In Fig. 16 we show the dependence of the BTE drift and Hall mobility on temperature, and compare our calculations to experimental data whenever available. The calculations are presented in a log-log scale. We indicate with dashed lines the linear fits used to determine the temperature exponents in Table VIII of the Appendix. Our results span a large range of power law decays, from T −1.18 for c-BN to T −3.15 for AlP. In the temperature range considered here (150-500 K), we find that a single exponent captures fairly well the calculated trends for all materials. Some deviations are found to occur above 400 K in the case of GaAs, AlP, c-BN, AlAs, and SrO. These deviations are indicative of highfrequency optical phonons becoming increasingly important in carrier scattering at these temperatures. We note that this fitting procedure is done in order to extract a single exponent which can be compared to experimental data fitted with the same procedure and in the same temperature range, whenever available.
In Fig. 17(a) we compare the calculated and measured temperature exponents. We find that in general the measured exponents are well reproduced by the calculations,  Table VIII and obtained on the 150 K to 500 K range. The references for the experimental data are as follows: diamond 146,173-177 , silicon [178][179][180][181][182] , GaAs 171,183-187 , 3C-SiC [188][189][190][191]AlP 192,193 , , and AlSb 201-203 . with an error of less than 20%. The only exceptions are found for SiC and AlSb, for which very few experimental studies exist. Given the overall good agreement found for the other compounds, we expect that the agreement between theory and experiment for SiC and AlSb will improve as new experimental data will be produced in the future.
In Fig. 17(b) we show the comparison between calculated and measured mobilities, in a logarithmic scale. In all cases the experimental data lay below the predicted phonon-limited Hall mobility. The difference can be as severe as 100% in some cases. The calculated drift mobility is found to be closer to the experimental data, but we have to emphasize that most experimental data reported in this figure are for the Hall mobility, as indicated in Table VIII. The SERTA drift mobility appears to be closer to experiments, and in some cases below the experimental data (for example in diamond, AlSb and GaAs). This behavior is probably due to a cancellation of errors, and should not be taken as an endorsement of the SERTA method. Whenever possible we recommend to perform complete BTE calculations. We note that, in calculations of phonon-limited carrier mobilities, it is not meaningful to validate computational methods by direct comparison with experiments, since experimental data incorporate additional scattering mechanisms.

F. Corrections to band structures and matrix elements from experimental data
In this section we investigate the role of many-body renormalization of the band structures and electronphonon matrix elements on the mobility. Since it is not currently possible to perform complete many-body calculations of the BTE, we estimate these effects by using the experimental effective masses and by rescaling the electron-phonon matrix elements via the experimental dielectric permittivities.
To this aim, starting from the experimental masses reported in Table V, For the holes we average the band contribution with the following weighted average formula which assumes band Using these experimentally-derived quantities, we rescale our calculated electron and hole Hall mobilities using:  The relative change in mobility obtained by this procedure is reported in Table VI. In all materials considered here this rescaling decreases the hole mobility, with the exception of SiC and BN. For the electron mobility, in all cases but diamond the mobility increases upon rescaling. Along the same line, we also proceed to renormalize the electron-phonon matrix elements by using the experimental dielectric screening. This strategy allows us to capture effects similar to the many-body renormalization of the electron-phonon coupling 5,166 . Since the mobility contains the square moduli of the matrix elements, we rescale the mobility using: In this operation we employ the dielectric constants reported in Table V.
In Fig. 18 we show the impact of rescaling the matrix elements (black triangles) on the mobility together with the rescaling of the matrix elements and effective masses (blue crosses). For all cases we also report the mean absolute error (MAE) and the mean absolute relative error (MARE) between the computed and experimental mobilities.
Our best first-principles BTE Hall mobilities have a large absolute error around 1000 cm 2 /Vs due to the strong overestimation of the high electron mobility of GaAs and with a MARE of 64 %. If we remove GaAs from the dataset, which is the compound with the highest mobility and hence it weighs disproportionately in the calculation of the average, we obtain a maximum absolute error of 335 cm 2 /Vs, and a MARE of 62%. When we proceed to renormalize the electron-phonon matrix ele- ments using the experimental dielectric constant, the results improve drastically with a MAE below 205 cm 2 /Vs and MARE of 35 %. However, this improvement is partially cancelled when we also use the experimental effective masses, leading to a MARE of 55 %.
These results indicate that improving the accuracy of the band structures and electron-phonon matrix elements, for example by using hybrid functionals 205 or many-body perturbation theory 5,166 will be important to achieve agreement with experiment. Whether these improvements will lead to close agreement with experiment, or else a more powerful theory beyond the BTE formalism will be needed in the future, remains an important open question.
Finally, we note that regardless of the level of approximation or re-scaling applied, the electron mobility of 3C-SiC and AlAs, and the hole mobility of GaP are significantly higher than experiment. This suggests that higher mobilities could be achieved in these materials, paving the way for improved electronic devices.

V. CONCLUSION
In this work we have carried out an extensive investigation of the carrier transport properties of 10 simple semiconductors including diamond, silicon, GaAs, 3C-SiC, AlP, GaP, c-BN, AlAs, AlSb, and SrO. Exploiting a precise interpolation approach and an efficient procedure for numerical convergence, we have been able to compute from first-principles transport coefficients with high accuracy. This key step has allowed us to directly quantify the impact of a wide range of approximations often used in transport calculations. In addition, we have been able to assess the role of a number of physical factors such as the contribution of long-range dipole and quadrupole interactions, the SOC, effective masses and lattice parameters as well as the effect of temperature and magnetic field. One of the most interesting result of this extensive study is that our most accurate Hall mobility calculations seem to systematically overestimate the experimental data. In a few cases the overestimation is as large as a factor of two. This is a quite promising outcome as in our calculation we only included electronphonon scattering and our result should be regarded as an upper limit to the Hall mobility achievable in highly pure samples, suggesting the possibility of further improvement in the transport properties of semiconductors of technological relevance. Here it is also important to stress that our work has clearly shown that the Hall factor exhibits significant variations across materials and as a function of temperature, ranging between 0.7 and 2. This is a clear indication that the common practice of comparing computed drift mobilities with experimental Hall mobilities should be done with particular care. Finally, it is important to stress that in our analysis we have also showed that part of the difference between theory and experiment may be connected with the inaccuracy in the DFT effective masses and electron-phonon matrix elements, but further work using explicit many-body perturbation theory calculations will be needed to confirm this point.
Our work demonstrates that first-principles calculations of carrier transport are making considerable progress. We believe we are reaching a point where detailed comparison with high-quality experimental data on pure samples will be increasingly important. For example, we highlight the need for new accurate measurements of the carrier mobility of SiC, AlAs, and GaP, and we hope to reignite interest in the experimental community to generate modern transport datasets for theorists to compare with. This manuscript describes a comprehensive protocol to perform such delicate calculations, and as a blueprint for reporting first-principles transport data in a way that is accessible, reliable, and reproducible. In this spirit, we encourage the adoption of the FAIR 106 data and software principles in the study first-principles transport. The interpolation of the Hamiltonian, dynamical matrix, and electron-phonon matrix elements from real space to reciprocal space via Eqs. (21), (22) and (23) is usually performed by using a radial cutoff on the direct lattice vectors 63,206 . While this is certainly the simplest possible choice, it is not the most efficient, especially in the case of unit cells with high aspect ratios.
In order to improve the accuracy of the interpolation, it is convenient to select the direct lattice vectors so that the distance between two Wannier centers or the distance between two atoms are within a cutoff radius. For example, in the interpolation of the Hamiltonian we construct a shifted Wigner-Seitz cell of lattice vectors R p with the requirement that this cell be centered on one of the two Wannier functions, say r n . The radial cutoff is then imposed with respect to the distance |R p + r m − r n |, where r m denotes the center of the second Wannier function.
In practice, we determine the set of Wigner-Seitz vectors to describes the electronic Hamiltonian: (A1) such that R p + r m − r n ∈ (i p × j p × o p ) where r m is the origin vector pointing to the position of the mth Wannier center. A similar approach was recently implemented in the Wannier90 software for the interpolation of the Hamiltonian 75 .
Similarly, in the case of the dynamical matrix we consider pairs of atoms τ κ and τ κ . For each such pair, we construct a Wigner-Seitz cell centered at τ κ , and we impose a truncation based on the distance |R p + τ κ − τ κ |. For the optimal choice of a Wigner-Seitz vector for a discretized Brillouin zone grid i p ×j p ×o p , the short-range part of the dynamical matrix takes the form 66 : for all direct lattice vectors R p such that R p +τ κ −τ κ ∈ (i p ×j p ×o p ) and 0 otherwise. In Eq. (A2) we have made a rigid shift to the central primitive cell such that it only depends on one lattice vector R p . As a result one needs to store a set of atom-dependent Wigner-Seitz cells. In first-principles software such as Quantum ESPRESSO 88 , a Fourier-transform grid is constructed with weights centered on pairs of atoms and zeros elsewhere. In EPW, we slightly optimize the computational cost of the interpolation by only retaining the union of nonzero elements between all the Wigner-Seitz cells.
In the case of the electron-phonon matrix elements, we exploit the same ideas as above, this time for both the Wannier functions and the atomic positions. The Wigner-Seitz cells for electrons and phonons are chosen so that both the |R p + r m − r n | and |R p + τ κ − r n | distances are minimized: 2M κ ω qν e † κα,qν U † mm k+q g m n ν (k, q)U n nk , (A3) where the lattice vectors are obtained such that R p + τ κ −r n ∈ (i p ×j p ×o p ) and R p +r m −r n ∈ (i p ×j p ×o p ). We note that the set of lattice vectors R p used for the interpolation of electron-phonon matrix elements can in general be different from the set employed for the interpolation of the dynamical matrix. This procedure used to determine the optimal set of vectors is sketched in Fig. 19. The set of direct lattice vectors R p and R p are selected such that the trial point t = R p + r m − r n , constructed from the distance between two Wannier centers, or t = R p + τ κ − r n constructed from a Wannier center and an atomic position in the central cell and any other cell, is closer to the supercell center 0 than any other supercell centers S. In the case of a trial point at equal distance between 0 and S (i.e., on the surface of the supercell Wigner-Seitz cell) a weight proportional to the number of degeneracies is given. The optimal Wigner-Seitz construction is activated with the new use ws input variable in EPW. As a sanity check, we verified that the mobility of 3C-SiC calculated on a 8×8×8 k-point and 4×4×4 q-point grid was the same as the one obtained from a 2×1×1 supercell with a 4×8×8 k-point and 2×4×4 q-point grid.
In Fig. 20 we report the calculated phonon dispersions of the three acoustic modes of 3C-SiC. The dashed gray line represents the interpolation from a coarse qpoint grid using a Γ-centered fixed Wigner-Seitz cell with dipole-dipole correction. The red line is for an optimal Wigner-Seitz cell with dipole-dipole correction. The blue line is for an optimal Wigner-Seitz cell with dipoledipole, dipole-quadrupole and quadrupole-quadrupole corrections. We can see that the Γ-centered Wigner-Seitz cell construction yields spurious soft phonon modes along the Γ−L direction. The improved Wigner-Seitz construction eliminates this artifact. Deformation potential of 3C-SiC without spinorbit coupling of the first electronic band (m = n = 1) on a coarse 3×3×3 k-point and q-point grid including dipole and quadrupole corrections. The Wigner-Seitz cell is either Γcentered or constructed using Rp + rm − rn ∈ (ip × jp × op) for the electrons and R p + τ κ − rn ∈ (i p × j p × o p ) for the phonons, respectively. The black dots are reference datapoints obtained via direct DFPT calculations.
We also show in Fig. 21 a comparison between the deformation potential of 3C-SiC computed with a Γcentered Wigner-Seitz cell or with the improved construction. The new construction produces results closer to those obtained directly with DFPT everywhere. We note that the results shown in Fig. 21 are for an unconverged small coarse grid of 3×3×3 for both kand q-points. When using denser coarse grids, these difference rapidly disappear (at least for the present case of cubic semiconductors). Given that we use much denser coarse grids in this work, for simplicity we employed the basic Γ-centered construction in all calculations reported in this work.     Effective quadrupole termQ βγ κα obtained from the least-mean square fit minimization of direct DFPT calculations for 3C-SiC. In these calculations we used the LDA exchange and correlation functional and without nonlinear core correction, so that we can compare with Ref. 22. The two non-equivalent values of the quadrupole tensor are obtained by using the phonon frequencies, Born-effective charges, dielectric constant, eigendisplacement vectors and wavefunction overlaps from explicit DFPT calculations, and by optimizing the quadrupole tensor using Eq. (40). We repeat the fitting for several sets of q-points on spherical surfaces of varying radii R. The radii are in units of 2π/a. The horizontal axis indicates the coarse q-point grid, and the lines correspond to different choices of the coarse k-point grid: same as the q-point grid (blue), twice as dense as the q-point grid (green), and three times as dense as the q-point grid (red). In all cases the fine grids were set to 60 3 points for both k and q. We also show with dashed lines the relative computational time needed in each case (normalized to the most time-consuming calculation). These tests show that the results obtained by using equal coarse grids for electrons and phonons are of the same accuracy as those with denser electron grids, but much cheaper computationally. Therefore we recommend to always use equal coarse grids.