Towards predictive many-body calculations of phonon-limited carrier mobilities in semiconductors

# Towards predictive many-body calculations of phonon-limited carrier mobilities in semiconductors

## Abstract

We probe the accuracy limit of ab initio calculations of carrier mobilities in semiconductors, within the framework of the Boltzmann transport equation. By focusing on the paradigmatic case of silicon, we show that fully predictive calculations of electron and hole mobilities require many-body quasiparticle corrections to band structures and electron-phonon matrix elements, the inclusion of spin-orbit coupling, and an extremely fine sampling of inelastic scattering processes in momentum space. By considering all these factors we obtain excellent agreement with experiment, and we identify the band effective masses as the most critical parameters to achieve predictive accuracy. Our findings set a blueprint for future calculations of carrier mobilities, and pave the way to engineering transport properties in semiconductors by design.

During the last decade, materials design guided by first-principles calculations has emerged as a powerful research strategy. Nowadays it is often possible to accurately predict ground-state properties of new materials in silico. This information can be used to screen for promising new materials Curtarolo et al. (2013); Jain et al. (2016). At variance with ground-state properties, the prediction and screening of materials properties involving electronic excitations is still in its infancy. For example charge and heat transport coefficients are typically evaluated using a combination of ab initio and semi-empirical approaches Wang et al. (2011); Chen et al. (2013); Hautier et al. (2013); Xu and Verstraete (2014); Krishnaswamy et al. (2017). The reasons for this lag are that the evaluation of transport coefficients is considerably more challenging than total energies, the computational infrastructure is not yet fully developed, and the lack of a clear set of reference data for validation and verification Lejaeghere et al. (2016).

In this work, we focus on phonon-limited carrier mobilities in semiconductors. The theoretical framework for calculating mobilities is well established, and is rooted in the Boltzmann transport equation (BTE), as described in Refs. Ziman (1960); Grimvall (1981); Mahan (2000). The BTE is a semiclassical, quasiparticle theory of electron transport, which can be rigorously derived from a many-body quantum-field theoretic framework by neglecting two-particle correlations Kadanoff and Baym (1962). The key ingredients are the electronic band structures, the phonon dispersion relations, and the electron-phonon matrix elements. The calculations of these quantities have reached maturity Giustino (2017), therefore there should be no fundamental obstacles towards predicting mobilities. However, already for the most studied semiconductor, silicon, one finds that (i) calculations of carrier mobilities are scarce, (ii) there is considerable scatter in the calculated data, and (iii) reproducing measured mobilities remains a challenge. For example, Refs. Zhou et al., 2015; Fiorentini and Bonini, 2016; Li, 2015; Restrepo et al., 2009 calculate intrinsic electron mobilities at room temperature , 1750, 1860, and 1970 cm/Vs, respectively, while experiments are in the range 1300-1450 cm/Vs Canali et al. (1975); Norton et al. (1973); Sze and Ng (2007).

Motivated by these considerations, here we set to clarify the accuracy limit and the predictive power of ab initio mobility calculations based on the BTE. We show that in order to correctly reproduce experimental data we need to take into account GW quasiparticle corrections to the band structures and the electron-phonon matrix elements, to include the spin-orbit splitting of the valence bands, and to properly converge the integrals over the Brillouin-zone. We also find that accurate band effective masses are absolutely critical to reproduce measured mobilities. By considering all these aspects, we succeed in reproducing measured data with high accuracy, thus establishing unambiguously the predictive power of the ab initio BTE.

In a semiconductor the steady-state electric current is related to the driving electric field via the mobility tensors as: , where Greek indices denote Cartesian coordinates. In this expression , and , are the mobility and particle density of electrons and holes, respectively. Within Boltzmann’s transport formalism Ziman (1960) the current density is expressed as , where and are the volume of the crystalline unit cell and the first Brillouin zone, respectively. The occupation factor plays the role of a statistical distribution function, and reduces to the Fermi-Dirac distribution in the absence of the electric field. The band velocity is given by , where is the single-particle electron eigenvalue for the state .

Using these definitions, the electron mobility is obtained via the derivative of the current with respect to the electric field: . Here the summations are restricted to the conduction bands, and is short for . An analogous expression holds for holes. From this expression we see that in order to calculate mobilities we need to evaluate , that is the linear response of the distribution function to the electric field . This quantity can be computed starting from the BTE Ziman (1960):

 (−e)E⋅1ℏ∂fnk∂k=2πℏ∑mν∫dqΩBZ|gmnν(k,q)|2 ×{(1−fnk)fmk+qδ(εnk−εmk+q+ℏωqν)(1+nqν) +(1−fnk)fmk+qδ(εnk−εmk+q−ℏωqν)nqν −fnk(1−fmk+q)δ(εnk−εmk+q−ℏωqν)(1+nqν) −fnk(1−fmk+q)δ(εnk−εmk+q+ℏωqν)nqν}. (1)

The left-hand side of Eq. (1) represents the collisionless term of Boltzmann’s equation for a uniform and constant electric field, in the absence of temperature gradients and magnetic fields; the right-hand side represents the modification of the distribution function arising from electron-phonon scattering in and out of the state , via emission or absorption of phonons with frequency , wavevector , and branch index  Grimvall (1981). is the Bose-Einstein distribution function. The matrix elements in Eq. (1) are the probability amplitude for scattering from an initial electronic state into a final state via a phonon , as obtained from density-functional perturbation theory Baroni et al. (2001); Giustino (2017). By taking derivatives of Eq. (1) with respect to we obtain an explicit expression for the variation :

 ∂Eβfnk=e∂f0nk∂εnkvnk,βτ0nk+2πτ0nkℏ∑mν∫dqΩBZ|gmnν(k,q)|2 ×[(1+nqν−f0nk)δ(εnk−εmk+q+ℏωqν) +(nqν+f0nk)δ(εnk−εmk+q−ℏωqν)]∂Eβfmk+q, (2)

having defined the relaxation time:

 1τ0nk = 2πℏ∑mν∫dqΩBZ|gmnν(k,q)|2 (3) × [(1−f0mk+q+nqν)δ(εnk−εmk+q−ℏωqν) + (f0mk+q+nqν)δ(εnk−εmk+q+ℏωqν)].

Equation (Towards predictive many-body calculations of phonon-limited carrier mobilities in semiconductors) is the linearized BTE and is valid under the assumption that the energy gained by a carrier accelerated by the electric field over the mean free path is much smaller than the thermal energy, ; this assumption is verified in most semiconductors under standard operating conditions. This equation needs to be solved self-consistently for , and is also referred to as the iterative BTE (IBTE). A simpler approach consists in neglecting the integral on the r.h.s. of Eq. (Towards predictive many-body calculations of phonon-limited carrier mobilities in semiconductors). In this case we obtain the variation without solving iteratively. It can be shown that the relaxation time is related to the imaginary part of the Fan-Migdal electron self-energy Giustino (2017) via . Based on this analogy, in the following we refer to the approximation of neglecting the integral in Eq. (Towards predictive many-body calculations of phonon-limited carrier mobilities in semiconductors) as the ‘self-energy relaxation time approximation’ (SERTA). In this approximation the mobility takes the simple form:

 μe,αβ=−eneΩ∑n∈CB∫dkΩBZ∂f0nk∂εnkvnk,αvnk,βτ0nk. (4)

We perform calculations within density-functional theory (DFT), planewaves, and pseudopotentials using the EPW code Poncé et al. (2016) of the Quantum ESPRESSO distribution Giannozzi et al. (2017), in conjunction with the wannier90 library Mostofi et al. (2014). This approach employs a generalized Wannier-Fourier interpolation technique Giustino et al. (2007) in order to obtain electron eigenvalues, phonon eigenfrequencies, and electron-phonon matrix elements on dense Brillouin zone grids by means of maximally localized Wannier functions Marzari et al. (2012). A fine sampling of the Brillouin zone is required because, at finite temperature, the Fermi level lies within the band gap, therefore we need to sample scattering processes taking place in the tails of the Fermi-Dirac distribution. In our calculations the Fermi level is determined in such a way that the net charge density at a given temperature, , equals the doping level ( for an intrinsic material). We now analyze in turn the key ingredients when calculating mobilities. We consider the paradigmatic case of silicon, for which extensive experimental data are available.

### .1 Brillouin-zone sampling

We find that in order to obtain reliable intrinsic mobilities it is necessary to employ extremely fine quasi-random grids, with a densified sampling around the band extrema. Convergence of mobility values to within 0.5% is reached when using grids with 85K inequivalent -points and 200K inequivalent -points [white dot in Fig. 4(a)]. Subsequent calculations in this article are performed using these grids. In Appendix Fig. 4(b) we compare calculations of the intrinsic mobility of silicon within the SERTA and the IBTE approaches. We find that the iterative solution of Eq. (Towards predictive many-body calculations of phonon-limited carrier mobilities in semiconductors) leads to converged values which are 6% higher than the SERTA result for electrons, and 1% lower for holes. Since the IBTE is drastically more expensive because it requires homogeneous and commensurate grids Li (2015); Fiorentini and Bonini (2016), in the following discussion we focus on SERTA calculations. We use a finite broadening of 5 meV to evaluate the Dirac delta function in Eq. (3). The sensitivity of the results to the broadening parameter is analyzed in Fig. 5 of the Appendix.

### .2 Exchange and correlation

In order to investigate the effect of the DFT exchange and correlation we perform calculations within both the local density approximation (LDA) Ceperley and Alder (1980); Perdew and Zunger (1981) and the generalized gradient approximation of Perdew, Burke, and Ernzerhof (PBE) Perdew et al. (1996), using scalar-relativistic pseudopotentials Hamann (2013). Figure 1 shows that the intrinsic mobilities at 300 K differ by 16% between LDA and PBE for electrons, and by 3% for holes. Closer inspection shows that these differences arise primarily from the optimized lattice parameters obtained with these functionals ( Å in LDA and 5.47 Å in PBE). In fact, when using the experimental lattice parameter ( Å) the deviation between LDA and PBE mobilities reduces to 0.4% for electrons and 2% for holes (Fig. 1). These results indicate that the choice of exchange and correlation is not critical so long as accurate lattice parameters are employed.

### .3 Spin-orbit coupling

Spin-orbit interactions in silicon are very weak Yu and Cardona (2010), therefore relativistic effects are usually neglected. However, here we find that spin-orbit coupling is important for predictive calculations, yielding hole mobilities 9% higher than non-relativistic calculations (Fig. 1). This effect can be understood by considering the band structures in Fig. 2(b). The spin-orbit interaction splits the six-fold degenerate states at the top of the valence bands, leading to the formation of two doubly-degenerate light-hole and heavy-hole bands, and one doubly-degenerate split-off hole band. As a result the effective mass of the light hole decreases (see Appendix Table 1), leading to a higher mobility. On the other hand, Fig. 2(a) shows that the conduction band bottom is relatively unaffected by spin-orbit coupling, and correspondingly the effect on the electron mobility is less pronounced (2.7%).

### .4 Many-body quasiparticle corrections

Given the sensitivity of the calculated mobilities to the band extrema, we investigate the effect of many-body correlations within the GW quasiparticle approximation. To obtain quasiparticle energies we use the Yambo code Marini et al. (2009); the values calculated on a 121212 uniform grid are then interpolated using the EPW code. Figure 2 shows the modification to the band extrema resulting from quasiparticle corrections. In the case of the valence bands, quasiparticle corrections increase the mass of the light holes (see Table 1 in Appendix); as a result the hole mobility decreases by 3%, as shown in Fig. 1. The opposite effect is observed for the conduction bands where the electron mobility is increased by 5%.

### .5 Corrections to the DFT screening

Another source of error in the DFT calculations of carrier mobilities is the overscreening of the electron-phonon matrix elements associated with the DFT band gap problem Giustino (2017). In fact, in the case of silicon DFT yields a static dielectric constant , which is higher than the measured value  Sze and Ng (2007). In order to overcome this issue it is necessary to modify the screening in the calculation of phonon dispersion relations. Since this is computationally prohibitive, here we take a simpler approach and renormalize the matrix elements as follows: . Here is meant to be the most accurate description of the screening that we can afford, and we are neglecting local-field effects which should yield an error on the order of a few percent Hybertsen and Louie (1987). For practical purposes we replace the dielectric functions by an analytic expression Resta (1977), where the only input parameter is the head of the dielectric matrix. The validity of this procedure is demonstrated in Appendix Fig. 6 using explicit calculations in the random-phase approximation. This correction to the matrix elements leads to a decrease of the electron and hole mobilities by 8.8% and 12.4%, respectively, as shown in Fig. 1.

### .6 Thermal expansion and electron-phonon renormalization

We computed the effect of thermal lattice expansion on the DFT eigenenergies using the thermo_pw code http://people.sissa.it/dalcorso/thermo_pw_dist.html () within the quasi-harmonic approximation and concluded that this effect is negligible, see Appendix Figs. 7-8. We also determined the electron-phonon renormalization of the effective masses using data from Ref. Poncé et al., 2015. This effect increases the masses by 3%, and results into a decrease of the mobilities by 5%.

After considering all the effects discussed so far, and after accounting for the corrections to the SERTA results arising from the solution of the complete IBTE, our most accurate theoretical mobilities at 300 K are  cm/Vs and  cm/Vs. These values are to be compared to the measured drift mobilities -1450 cm/Vs Ludwig and Watters (1956); Cronemeyer (1957); Li and Thurber (1977); Jacoboni et al. (1977) and -510 cm/Vs Dorkel and Leturcq (1981); Jacoboni et al. (1977); Ludwig and Watters (1956); Cronemeyer (1957) (Fig. 1). From the comparison with experiment we see that by pushing the theory to its limits we can obtain electron mobilities in very good agreement with experiment. On the contrary, the hole mobility are still approximately 30% above the measured range. This discrepancy can be traced back to the underestimation of the [100] heavy hole effective masses within the GW approximation. In fact, by repeating the calculation using the experimental hole effective mass instead of the GW mass, we obtain a hole mobility  cm/Vs, this time in very good agreement with experiment as shown in Fig. 1. This result leads us to conclude that the effective mass plays an absolutely critical role in mobility calculations. Our finding can be understood by considering that the mobility varies with the effective mass as with being a coefficient between 1 and 2.5 Bardeen and Shockley (1950); Blatt (1957); Keyes (1959); as a result a 20% error in the effective mass leads to an error in the mobility of up to 60%. This finding highlights the critical role of accurate calculations of quasiparticle band structures, and raises the question on whether the standard GW method and pseudopotential calculations (see Table 2) are sufficient for delivering predictive mobilities.

Using the best possible computational setup we can now compare our calculations with experiment over a range of temperatures and doping levels. Figure 3(a) shows the intrinsic electron and hole mobilities of silicon between 100 K and 500 K. In the case of the hole mobilities we show both our best ab initio results (solid line), as well as those re-calculated using the experimental effective masses (dashed line). Overall, the agreement between our calculations and experiment is very good throughout the entire temperature range. Figure 3(b) shows a comparison between calculated and measured mobilities at 300 K, as a function of carrier concentration between and  cm. In this case, in addition to the ab initio electron-phonon scattering, we used the semi-empirical model of Brooks and Herring with the Long-Norton correction Brooks (1951); Li and Thurber (1977) to account for impurity scattering (see Appendix for details). Also in this case we find very good agreement with experiment, although the contribution of impurity scattering is evaluated semi-empirically.

In conclusion, we pushed the accuracy of transport calculations within the BTE formalism to its limits, and we demonstrated that this approach can deliver predictive accuracy for a prototypical semiconductor. Our findings raise two important questions for future work on transport in semiconductors: (i) the present formalism yields results which fall within the experimental uncertainty. In order to enable further progress in this area it will be important to produce a high-quality experimental data from single-crystal samples. (ii) An unexpected challenge that we faced is to perform accurate ab initio calculations of effective masses. Going forward it will be important to establish whether the GW method and pseudopotential calculations can provide effective masses with the accuracy required for predictive mobility calculations. Meanwhile, the present work opens the way to predictive calculations of mobilities and lays the groundwork for the ab initio design of semiconductor devices.

After submission of this work, a related calculation for Si was reported, where the authors found a significant increase in Si hole mobility with SOC and no effect from SOC on the electron mobility in line with our results Ma et al. (2018).

###### Acknowledgements.
We acknowledge fruitful discussions with C. Verdi, M. Schlipf and W. Li. This work was supported by the Leverhulme Trust (Grant RL-2012-001), the UK EPSRC Research Council (grants No. EP/J009857/1 and EP/M020517/1), the EU H2020 programme under grant No. 696656 GrapheneCore1, the University of Oxford Advanced Research Computing (ARC) facility (http://dx.doi.org/810.5281/zenodo.22558), the ARCHER UK National Supercomputing Service under the AMSEC and CTOA projects, PRACE DECI-13 resource Cartesius at SURFsara, and the PRACE DECI-14 resource Abel at UiO. E.R.M. acknowledges the NSF support (Award No. OAC-1740263).

## I Appendix

### i.1 Computational Methods

In this work we use norm-conserving pseudopotentials with planewave kinetic energy cutoffs of 45 Ry and 35 Ry for LDA and PBE calculations, respectively. The phonon dispersion relations are evaluated using density-functional perturbation theory Baroni et al. (2001), starting from a 666 uniform grid of -points. An 181818 uniform grid of -points was required to correctly obtain vanishing Born effective charges. Representative phonon dispersion relations obtained within this setup can be found in Ref. Poncé et al. (2016). The coarse grids for the electron-phonon interpolation required 121212 -points and 666 -points. Such a dense -grid was needed to obtain a good Wannier interpolation of the conduction bands, since the minimum is along the line (approx. 0.85X) and does not fall on a high-symmetry point.

For the self-consistent, iterative solution of the Boltzmann transport equation (IBTE) we employ uniform Brillouin-zone grids, and the -point sums are restricted to the irreducible wedge of the Brillouin zone using crystal symmetry operations. The IBTE is solved using homogeneous and commensurate - and -point grids since the variations at the -th iteration require the knowledge of the variations at the -th iteration, see Eq. (Towards predictive many-body calculations of phonon-limited carrier mobilities in semiconductors).

For the direct solution of the BTE within the self-energy relaxation time approximation (SERTA) the Brillouin zone grids do not need to be commensurate. In this case, in order to improve the sampling accuracy, we employ quasi-random Sobol sequences of - and -points. Following recommended practice, we skip the first 1000 elements of a sequence and we retain one element every 100 of the remainder Bratley and Fox (1988); furthermore we employ a linear scramble and shift of the resulting sequence, using standard routines from Matlab R2015a Hong and Hickernell (2003). As a further refinement we replace the homogeneous Sobol weights using a Voronoi triangulation with the code Voro++ Rycroft (2009). In the Voronoi triangulation we take into account the periodicity of the Brillouin zone by building periodic replicas of the random grid in neighboring reciprocal unit cells. For the -point grid we also densify the distribution around the band extrema, in order to capture the fine features of the scattering near the band edges. This is achieved by generating additional random points with the Lorentz distribution and by recomputing the Voronoi weights of the resulting grid. Here indicates the location of the band extrema and  Å.

Figure 4(a) shows the convergence of the intrinsic mobility of silicon at 300 K with respect to the number of electron and phonon wavevectors in the Brillouin zone within the SERTA approximation. Figure 4(b) shows the comparison between calculations of the intrinsic mobility of silicon within the SERTA and the IBTE approaches.

The GW calculations are performed starting from the PBE band structure and using the experimental lattice parameters on a 121212 -point grid. To obtain direct and indirect band gaps converged to within 5 meV we use 120 bands and a planwaves cutoff of 15 Ry for the dielectric matrix. The renormalization of the band velocity is evaluated as in Ref. Rohlfing and Louie (2000): , where indicates the momentum operator. When the previous expression is replaced by .

For completeness the effective masses computed within scalar-relativistic DFT, fully-relativistic DFT, and including GW quasiparticle corrections are reported in Table 1. We also show in Table 2 the effective masses calculated without SOC at the experimental lattice parameter with two different types of pseudization (norm-conserving and ultrasoft), and two exchange and correlation functionals (LDA and PBE).

In order to calculate mobilities using band structures as close as possible to experiments (i.e. the lowermost bars in Fig. 1), we repeated the calculations using the low-energy dispersion relations parametrized in Refs. Dresselhaus et al., 1955; Yu and Cardona, 2010 starting from the measured effective masses:

 εcb= ℏ2(kx−k0,x)22m||+ℏ2(ky−k0,y)22m⊥ +ℏ2(kz−k0,z)22m⊥+εc, (5) εhh= Ak2+[B2k4+C2(k2xk2y+k2yk2z+k2zk2x)]1/2, (6) εlh= Ak2−[B2k4+C2(k2xk2y+k2yk2z+k2zk2x)]1/2, (7) εso= −k2ℏ22mso−εso, (8)

where ( is the free electron mass), , , denotes the wavevectors of the conduction band minima, and is the conduction band bottom. The coefficients are , and  Dresselhaus et al. (1955); Yu and Cardona (2010) and = 48 meV.

### i.2 Broadening of Dirac delta functions

The numerical evaluation of phonon-limited mobilities using Eqs. (Towards predictive many-body calculations of phonon-limited carrier mobilities in semiconductors)-(4) requires one to replace the Dirac delta functions in Eqs. (Towards predictive many-body calculations of phonon-limited carrier mobilities in semiconductors)-(3) by Lorentzian functions with finite broadening : . This procedure makes the calculated mobility dependent on the broadening parameter, hence it is important to check how sensitive are the results to the choice of .

Figure 5(a) shows the intrinsic electron mobility of silicon at 0 K, evaluated as a function of . From this figure we see that the mobility tends to diverge towards as . This trend can be rationalized by noting that the mobility is directly proportional to the relaxation time [cf. Eq. (4)], and the relaxation time due to acoustic phonon scattering in a non-polar semiconductor is inversely proportional to the temperature Bardeen and Shockley (1950). As a result, we expect that the phonon-limited mobility will increase indefinitely as becomes smaller and the Lorentzian approaches the Dirac delta function. This observation is in agreement with the explicit calculations in Fig. 5(a).

This behavior poses a problem when one has to decide which broadening parameter to use in the calculations. As a general rule here we set to the smallest possible value where the curve vs.  is relatively flat, so that our results are insensitive to this choice. Based on Fig. 5(a), we use  meV in all calculations presented in the article. This choice is consistent with the notion that real quasiparticles do not have an infinite lifetime as it is assumed in the BTE formalism, but have a finite lifetime due to electron-electron and electron-phonon interactions. In Fig. 5(b) we show our calculated quasiparticle broadening from electron-phonon interactions at 0 K and 300 K. It can be seen that at 300 K the broadening reaches values up to 4-5 meV for quasiparticle energies located one phonon energy away from the band bottom (the highest phonon energy in silicon is 63 meV). These values are consistent with our choice of broadening parameter.

### i.3 Screening of the electron-phonon matrix elements

The strategy that we used to correct for the DFT overscreening of the electron-phonon matrix elements consists of un-screening the matrix elements via the DFT dielectric function, so as to obtain the bare matrix elements, and then screening the bare matrix elements using the best possible dielectric function. The dielectric function can be factored out of the integral in the matrix element if we neglect local field effects. Since local field effects are known to decrease the head of the dielectric function of silicon by 10% and the body of the dielectric function is typically one or two orders of magnitude smaller than the head Hybertsen and Louie (1987), we expect to make an error on the order of a few percent.

In order to perform this operation for a large number of phonon wavevectors we use the Thomas-Fermi model dielectric function of Ref. Resta, 1977:

 ϵ(q)=k20+q2k20sin(qR)/(qRϵ0)+q2, (9)

where , is the macroscopic (electronic) dielectric constant. are are obtained from the valence electron density as and . The only free parameter of the model is .

We test the validity of this model by computing the dielectric matrix within the random phase approximation (RPA), using a 121212 unshifted grid (corresponding to 72 inequivalent wavevectors). Figure 6 shows a comparison between the model dielectric function of Eq. (9) and the RPA calculation, after matching to the head of the RPA dielectric matrix. We see that Eq. (9) reproduces well the RPA screening, therefore it is sensible to use it in the renormalization of the electron-phonon matrix elements. In this work we renormalized the matrix elements by setting to the experimental dielectric constant of silicon (11.94).

### i.4 Brooks-Herring model for impurity scattering

In order to account for impurity scattering in Fig. 3(b), we use the semi-empirical model developed by Brooks and Herring Brooks (1951); Li and Thurber (1977). In this model the mobility is evaluated analytically by taking into account quantum-mechanical scattering rates, spherical energy surfaces, negligible electron-electron interactions, and complete ionization of the impurities. The explicit expression of the hole mobility is:

 (10)

where , , and . Here is the density-of-state effective mass for the holes Balkanski and Wallis (2000), and are the hole densities and the density of ionized impurities [impurity concentration in Fig. 3(b)], respectively, is the dielectric constant, is the permittivity of vacuum, and is Planck’s constant. In the above expressions, the concentrations are expressed in cm, and the temperature is in K.

In the case of silicon, Eq. (10) cannot be used for the electron mobility because the electron mass is highly anisotropic and leads to incorrect results. To account for the electron mass anisotropy, we instead used the Long-Norton mobility expression Norton et al. (1973); Li and Thurber (1977):

 μLNi=7.3⋅1017T3/2niG(b)[cm2Vs], (11)

where the electron density-of-state effective mass is  Sze and Ng (2007).

Finally, the mobility including phonon () and impurity () scattering can be computed using the mixed-scattering formula Li and Thurber (1977):

 μ=μl[1+X2{ci(X)cos(X)+sin(X)(si(X)−π2)}], (12)

where and ci(X) and si(X) are the cosine and sine integrals.

### i.5 Effect of thermal lattice expansion

Within the quasi-harmonic approximation Baroni et al. (2010), the Helmholtz free energy of a cubic crystal is given by Palumbo and Corso (2017):

 F(T,V)=U(V)+Fvib(T,V)+Fel(T,V), (13)

where is the static energy at 0 K, is the contribution due to lattice vibration and the energy due to electronic thermal excitations. We rely on the adiabatic approximation to treat each term independently. The vibrational Helmholtz free energy per cell is given in the harmonic approximation by Palumbo and Corso (2017):

 Fvib(T,V)= 12N∑q,νℏωq,ν(V) + kBTN∑q,νln[1−exp(−ℏωq,ν(V)kBT)], (14)

where is the number of -points, the first term is the contribution to the zero-point energy and the second term is the phonon contribution at finite temperature. can be neglected as the band gap is much larger than thermal energies.

The energy minimum of at a given temperature corresponds to zero pressure and gives the variation of volume with temperature due to thermal expansion. To perform those calculations we used the thermo_pw code http://people.sissa.it/dalcorso/thermo_pw_dist.html (); Corso (2016). The phonon frequencies were computed using the same LDA and PBE pseudopotentials as in the manuscript, without spin-orbit coupling, at nine different volumes. The resulting energies were fitted using the Murnaghan equation of state Murnaghan (1944). We used a 181818 -point grid for the electron and a 666 -point grid for the phonons. The obtained volume variation is given in Fig. 7.

The change of eigenenergies due to thermal expansion is given by Lautenschlager et al. (1985):

 Δεnk(T)=−∂εnk∂P∣∣T∫T0dT′3α(T′)B(T′), (15)

where is the bulk modulus and is the thermal expansion coefficient. and are obtained via numerical differentiation starting from the volumes calculated in the above figure.

We note that in Eq. 15 we carried the term out of the temperature integral. This common temperature-independent approximation is valid in the elastic regime. To make sure that this approximation is valid, we compute at 4 K and 300 K by numerical derivation around the equilibrium volume for that temperature. From Table 3 we see that indeed the temperature dependence is negligible and we therefore use the value at 4 K in Eq. 15.

The bulk modulus, the thermal expansion and eigenstates renormalization with temperature of silicon computed with the LDA and PBE exchange-correlation functionals are presented in Fig. 8. From the bottom panel we see that thermal lattice expansion leads to a slight increase of the band gap of silicon. For PBE, the valence band top and conduction band bottom change by 9.5 meV and 7.8 meV from 0 K to 300 K, therefore the net increase of the band gap is 1.7 meV. This variation is much smaller than the gap renormalization arising from electron-phonon interactions (as discussed next), therefore in the present case this effect can safely be neglected when calculating carrier mobilities.

### i.6 Electron-phonon renormalization of the bandstructures and free carrier screening

The electron-phonon renormalization of the bandstructure has been discussed for the case of silicon in considerable detail in Ref. Poncé et al., 2015. The calculated zero-point renormalization of the fundamental gap is 56.2 meV within the non-adiabatic Rayleigh-Schrödinger perturbation theory. This change corresponds to 5% of the band gap. We extracted the effective masses using data from that paper, and the results are shown in Table 4 for specific directions.

The electron-phonon renormalization of the bands leads to an increase of both electron and hole effective masses between 1% and 6%. Therefore, we can reasonably estimate that similar changes in effective masses will occurs in our calculations. Given the dependence of the mobility on effective mass, we estimate a 5% reduction in mobility due to this effect.

The effect of free carrier screening can be included via a Lindhard dielectric function using ab initio parameters as described in Ref. Verdi et al., 2017. In this case the screening only affects phonons with energy below the plasma energy of the doped carriers. For intrinsic silicon, which is the main focus of our work, the carrier concentration is below 10 cm. Using data from Ref. Caruso and Giustino, 2016, we estimate that the plasma energy in this case would be well below 0.1 meV, therefore the free carrier screening would be ineffective for nearly all phonons. We also mention that the renormalization of the band structure arising from the free carriers is negligible in this case Caruso and Giustino (2016).

### References

1. S. Curtarolo, G. L. W. Hart, M. B. Nardelli, N. Mingo, S. Sanvito,  and O. Levy, Nature Materials 12, 191 (2013).
2. A. Jain, Y. Shin,  and K. A. Persson, Nature Reviews Materials 1, 1 (2016).
3. S. Wang, Z. Wang, W. Setyawan, N. Mingo,  and S. Curtarolo, Phys. Rev. X 1, 021012 (2011).
4. X. Chen, D. Parker,  and D. J. Singh, Scientific Reports 3, 3168 (2013).
5. G. Hautier, A. Miglio, G. Ceder, G.-M. Rignanese,  and X. Gonze, Nature Communications 4, 2292 (2013).
6. B. Xu and M. J. Verstraete, Phys. Rev. Lett. 112, 196603 (2014).
7. K. Krishnaswamy, B. Himmetoglu, Y. Kang, A. Janotti,  and C. G. Van de Walle, Phys. Rev. B 95, 205202 (2017).
8. K. Lejaeghere, G. Bihlmayer, T. Björkman, P. Blaha, S. Blügel, V. Blum, D. Caliste, I. E. Castelli, S. J. Clark, A. Dal Corso, S. de Gironcoli, T. Deutsch, J. K. Dewhurst, I. Di Marco, C. Draxl, M. Dułak, O. Eriksson, J. A. Flores-Livas, K. F. Garrity, L. Genovese, P. Giannozzi, M. Giantomassi, S. Goedecker, X. Gonze, O. Grånäs, E. K. U. Gross, A. Gulans, F. Gygi, D. R. Hamann, P. J. Hasnip, N. A. W. Holzwarth, D. Iuşan, D. B. Jochym, F. Jollet, D. Jones, G. Kresse, K. Koepernik, E. Küçükbenli, Y. O. Kvashnin, I. L. M. Locht, S. Lubeck, M. Marsman, N. Marzari, U. Nitzsche, L. Nordström, T. Ozaki, L. Paulatto, C. J. Pickard, W. Poelmans, M. I. J. Probert, K. Refson, M. Richter, G.-M. Rignanese, S. Saha, M. Scheffler, M. Schlipf, K. Schwarz, S. Sharma, F. Tavazza, P. Thunström, A. Tkatchenko, M. Torrent, D. Vanderbilt, M. J. van Setten, V. Van Speybroeck, J. M. Wills, J. R. Yates, G.-X. Zhang,  and S. Cottenier, Science 351 (2016).
9. J. Ziman, Electrons and Phonons, edited by N. Mott, B. E.C.,  and W. D.H. (Oxford University Press, 1960).
10. G. Grimvall, The electron-phonon interaction in metals (North-Holland Publishing Company, 1981).
11. G. D. Mahan, Many-Particle Physics (Springer, 2000).
12. L. P. Kadanoff and G. Baym, Quantum Statistical Mechanics (Benjamin, 1962).
13. F. Giustino, Rev. Mod. Phys. 89, 015003 (2017).
14. J. Zhou, B. Liao, B. Qiu, S. Huberman, K. Esfarjani, M. S. Dresselhaus,  and G. Chen, PNAS 112, 14777 (2015).
15. M. Fiorentini and N. Bonini, Phys. Rev. B 94, 085204 (2016).
16. W. Li, Phys. Rev. B 92, 075405 (2015).
17. O. D. Restrepo, K. Varga,  and S. T. Pantelides, Appl. Phys. Lett. 94, 212103 (2009).
18. C. Canali, C. Jacoboni, F. Nava, G. Ottaviani,  and A. Alberigi-Quaranta, Phys. Rev. B 12, 2265 (1975).
19. P. Norton, T. Braggins,  and H. Levinstein, Phys. Rev. B 8, 5632 (1973).
20. S. Sze and K. K. Ng, Physics of semiconductor Devices - Third Edition (Wiley, 2007).
21. S. Baroni, S. de Gironcoli, A. Dal Corso,  and P. Giannozzi, Rev. Mod. Phys. 73, 515 (2001).
22. S. Poncé, E. R. Margine, C. Verdi,  and F. Giustino, Computer Physics Communications 209, 116 (2016).
23. P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. de Gironcoli, P. Delugas, R. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. KüÃ§zükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O. de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu,  and S. Baroni, Journal of Physics: Condensed Matter 29, 465901 (2017).
24. A. A. Mostofi, J. R. Yates, G. Pizzi, Y.-S. Lee, I. Souza, D. Vanderbilt,  and N. Marzari, Computer Physics Communications 185, 2309 (2014).
25. F. Giustino, M. L. Cohen,  and S. G. Louie, Phys. Rev. B 76, 165108 (2007).
26. N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza,  and D. Vanderbilt, Rev. Mod. Phys. 84, 1419 (2012).
27. D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 566 (1980).
28. J. P. Perdew and A. Zunger, Phys. Rev. B 23 (1981).
29. J. P. Perdew, K. Burke,  and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
30. D. R. Hamann, Phys. Rev. B 88, 085117 (2013).
31. P. Y. Yu and M. Cardona, Fundamental of Semiconductors, edited by H. E. Stanley and W. T. Rhodes (Springer, 2010).
32. A. Marini, C. Hogan, M. GrÃ¼ning,  and D. Varsano, Computer Physics Communications 180, 1392 (2009).
33. F. J. Morin and J. P. Maita, Phys. Rev. 96, 28 (1954).
34. R. A. Logan and A. J. Peters, Journal of Applied Physics 31, 122 (1960).
35. G. W. Ludwig and R. L. Watters, Phys. Rev. 101, 1699 (1956).
36. C. Jacoboni, C. Canali, G. Ottaviani,  and A. A. Quaranta, Solid-State Electronics 20, 77 (1977).
37. H. Brooks, Phys. Rev. 83, 879 (1951).
38. S. S. Li and W. R. Thurber, Solid-State Electronics 20, 609 (1977).
39. M. S. Hybertsen and S. G. Louie, Phys. Rev. B 35, 5585 (1987).
40. R. Resta, Phys. Rev. B 16, 2717 (1977).
41. http://people.sissa.it/dalcorso/thermo_pw_dist.html,  .
42. S. Poncé, Y. Gillet, J. Laflamme Janssen, A. Marini, M. Verstraete,  and X. Gonze, J. Chem. Phys. 143, 102813 (2015).
43. D. C. Cronemeyer, Phys. Rev. 105, 522 (1957).
44. J. Dorkel and P. Leturcq, Solid-State Electronics 24, 821 (1981).
45. J. Bardeen and W. Shockley, Phys. Rev. 80, 72 (1950).
46. F. J. Blatt, Solid State Physics 4, 199 (1957).
47. R. W. Keyes, Journal of Applied Physics 30, 454 (1959).
48. J. Ma, A. S. Nissimagoudar,  and W. Li, Phys. Rev. B 97, 045201 (2018).
49. P. Bratley and B. L. Fox, ACM Trans. Math. Softw. 14, 88 (1988).
50. H. S. Hong and F. J. Hickernell, ACM Trans. Math. Softw. 29, 95 (2003).
51. C. H. Rycroft, Chaos 19 (2009).
52. M. Rohlfing and S. G. Louie, Phys. Rev. B 62, 4927 (2000).
53. G. Dresselhaus, A. F. Kip,  and C. Kittel, Phys. Rev. 98, 368 (1955).
54. R. N. Dexter and B. Lax, Phys. Rev. 96, 223 (1954).
55. M. Balkanski and R. F. Wallis, Semiconductor Physics and Applications (Oxford University Press, 2000).
56. S. Baroni, P. Giannozzi,  and E. Isaev, Rev. Miner. Geochem. 71, 39 (2010).
57. M. Palumbo and A. D. Corso, Journal of Physics: Condensed Matter 29, 395401 (2017).
58. A. D. Corso, Journal of Physics: Condensed Matter 28, 075401 (2016).
59. F. D. Murnaghan, Proc. Natl. Acad. Sci. USA 30, 244 (1944).
60. P. Lautenschlager, P. B. Allen,  and M. Cardona, Phys. Rev. B 31, 2163 (1985).
61. H. J. McSkimin and P. A. Jr., Journal of Applied Physics 35, 2161 (1964).
62. Y. Okada and Y. Tokumaru, Journal of Applied Physics 56, 314 (1984).
63. C. Verdi, F. Caruso,  and F. Giustino, Nature Communications 8, 15769 (2017).
64. F. Caruso and F. Giustino, Phys. Rev. B 94, 115208 (2016).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters