EPIC STAR: a reliable and efficient approach for phonon- and impurity-limited charge transport calculations

A computationally efficient first-principles approach to predict intrinsic semiconductor charge transport properties is proposed. By using a generalized Eliashberg function for short-range electron–phonon scattering and analytical expressions for long-range electron–phonon and electron–impurity scattering, fast and reliable prediction of carrier mobility and electronic thermoelectric properties is realized without empirical parameters. This method, which is christened “Energy-dependent Phonon- and Impurity-limited Carrier Scattering Time AppRoximation (EPIC STAR)” approach, is validated by comparing with experimental measurements and other theoretical approaches for several representative semiconductors, from which quantitative agreement for both polar and non-polar, isotropic and anisotropic materials is achieved. The efficiency and robustness of this approach facilitate automated and unsupervised predictions, allowing high-throughput screening and materials discovery of semiconductor materials for conducting, thermoelectric, and other electronic applications.


INTRODUCTION
Computational screening of electronic materials, empowered by new theories, computational tools and high computing power, is an efficient and economical alternative of purely experimental exploration. For instance, these tools have assisted the discovery of many new high performance thermoelectric (TE) candidates, such as TmAgTe 2 1 , NbFeSb 2 , NbCoSn 3 , and CdIn 2 Te 4 4 , to name a few. High mobility oxides 5 have also been explored to uncover new materials as transparent conductors with potential application in display, light-emitting diodes, and solar cells. Some computationally identified candidates, such as Ba 2 BiTaO 6 6 and TaIrGe 7 , have been experimentally verified. Many of the previous works are based on the prediction of carrier mobility μ, electrical conductivity σ, and Seebeck coefficient S. Although the state-ofthe-art tools and theories based on Boltzmann transport equation (BTE) are able to accurately predict both electron 8,9 and phonon 10,11 transport properties, the high complexity and computational cost make it a challenge for high-throughput screening of possible materials. To overcome this, various approximations have been established to make the calculations more tractable.
To determine μ, σ, and S, BTE solution within scattering time approximation (STA) is commonly adopted 12 which yields good agreement with fully numerical solution of the equation 13 . However, to determine the phonon-limited intrinsic transport properties, very dense sampling of both electron and phonon Brillouin zones are required for the integration. 10 5 -10 6 k and q points are needed for a material as simple as Si 9,13,14 and even more for polar semiconductors such as GaAs [14][15][16] . This is computationally complex and expensive even with the use of linear 13 or Wannier 8,17 interpolation techniques. Therefore, constant scattering time approximation (CSTA) is commonly employed in tools such as BoltzTraP 18,19 and BoltzWann 20 where only electronic structure calculation from density functional theory (DFT) 21 is needed, and it indeed provides a lot of insights into the role of the electron band structure. For example, high band degeneracy was found to be related to high zT, whereas band convergence, which increases such degeneracy, has become an effective method for enhancing TE performance 1,22,23 . Quantitative descriptors calculated in CSTA such as Fermi surface complexity factor 24 and electronic fitness function 25 have also been proposed as screening criteria for TE materials.
While CSTA could predict the Seebeck coefficient, which is independent of scattering time τ in CSTA, the carrier mobility and electrical conductivity inevitably vary as the magnitude of τ differs from one system to another. It is also found that when τ depends on the electron energy E, this dependency would also affect the value of S [26][27][28][29] . As such, this hinders the predictive power of CSTA, implying the need for evaluation of τ from first-principles by including important scattering mechanisms such as electron-phonon interaction and charged impurity scattering. Deformation potential theory (DPT), pioneered by Bardeen and Shockley 30 , describes the scattering from long-wavelength acoustic phonons, which is an important scattering mechanism in semiconductors. It has been widely used in predicting TE properties of various materials [31][32][33][34][35][36] and recently in highthroughput search for new TE compounds 4 . The electron-phonon averaged (EPA) approximation 37 takes into account all electron-phonon interactions and uses the average coupling to evaluate τ. It achieves quantitative agreement with full and accurate calculations using electron-phonon Wannier (EPW) method 8 for some half-Heusler compounds in the heavily doped region. However, the absence of polar optical phonon (POP) scattering impairs the accuracy of τ predictions, whereby the scattering time and performance may be overestimated for certain systems. For example, in lightly doped half-Heusler compounds, optical phonons dominate the scattering of charge carriers near band edges as a consequence of weak electron-acoustic-phonon coupling 38 . Therefore the intrinsic carrier mobility may be overestimated. In fact, we will show that EPA predicts an intrinsic hole mobility of 149 cm 2 V −1 s −1 for NbFeSb which is twice the value from EPW calculation 38 (77 cm 2 V −1 s −1 ). Direct averaging of electron-phonon matrix elements and frequencies may not be appropriate for acoustic phonons and POP due to large variations of their magnitudes in Brillouin zone 39 . Moreover, using one single frequency for a whole phonon branch is not accurate especially for acoustic phonons, where long wavelength electron-phonon scattering is almost elastic. Therefore, a more general approximation is needed to incorporate the appropriate treatment for all phonon branches. The short-range electron-phonon interaction in metals, which is important for the prediction of conventional superconductor and metal resistivity 40 , was described by Eliashberg electron-phonon spectral function α 2 FðωÞ which was used in Allen-Dynes formula for transition temperature 41,42 and Ziman's resistivity formula 12,43 . Its dependence on phonon frequency ω could reflect some details that are missing in DPT and EPA methods, and first-principles calculation of α 2 FðωÞ has been demonstrated for systems like Nb 44,45 , Al 13 , Pb 8,45 , and B-doped diamond 46 . Inspired by these works, the extension of Eliashberg function to non-metallic systems could potentially be a feasible approach to realize fast and accurate prediction of short-range electron-phonon interactions and charge transport properties. And the long-range electron-phonon scattering can be treated analytically due to their diverging behavior for small momentum phonons 47 . The charged impurity scattering could also be treated in a similar way as a Coulomb potential from point charges.
In this work, we propose a new first-principles approach for "Energy-dependent Phonon-and Impurity-limited Carrier Scattering Time AppRoximation (EPIC STAR)" from density functional perturbation theory (DFPT) calculations, which overcomes the above-mentioned difficulties with reliable prediction for electron scattering from all phonon branches. The EPIC STAR calculation after obtaining DFPT phonon and electron-phonon information is very fast and computational cost is only limited by DFPT itself. This new approach treats the short-and long-range scattering mechanisms separately to reproduce their different behavior. Specifically, it employs analytical expressions to treat the longrange POP scattering and impurity scattering, using ab initio information including Born effective charges and dielectric tensor; and the remaining short-range electron-phonon interactions are described numerically using a generalized Eliashberg function, α 2 FðE; ωÞ from first-principles calculations which depends on both electron energy E and phonon frequency ω. The long-range POP scattering formulism is based on recent works for first-principles treatment of Fröhlich interactions in polar materials 47,48 , and impurity scattering is simulated with Brooks-Herring model 39,49,50 . α 2 FðE; ωÞ is a generalization of the electron-phonon spectral function where E-dependence and energy conservation are explicitly included. The new approach is validated for Si, GaAs, Mg 2 Si, and NbFeSb, which represent high mobility and high TE performance materials with different extents of polarization and band structure anisotropy. EPIC STAR calculation also reveals a new TE material, NaInSe 2 whose performance has not been studied before. The agreement with EPW prediction and experimental measurements is investigated and high predictive capability is achieved at costs lower by at least 1 order of magnitude. The computational cost is not just similar to those of phonon calculations and the EPA method, but also reliable for more general family of materials with improved accuracy. Moreover, the DFPT calculations for EPIC STAR will not be wasted as they can be directly fed into accurate calculations like EPW if necessary, while the phonon-related information can be utilized in other relative studies such as thermal, elastic and piezoelectric properties. This approach is therefore suitable for automated charge transport property calculation and high-throughput screening of electronic materials with different functionalities arising from non-equilibrium electron scattering times, such as thermoelectrics, high power electronics and photovoltaics.

RESULTS AND DISCUSSIONS
Energy-dependent scattering time approximation The electron energy scattering time limited by electron-phonon interactions is given by 40 where g nml k; q ð Þ is the electron-phonon coupling matrix element among electrons in state nk, state mk þ q, and phonon in branch l with momentum q. n 0 lq is the Bose-Einstein distribution function of phonons with frequency ω lq while f 0 mkþq is the Fermi-Dirac distribution function of electrons with energy ε mkþq . The phonon frequencies and electron-phonon coupling matrix can be obtained from density functional perturbation theory (DFPT) calculation as implemented in PHonon code in QUANTUM ESPRESSO 51,52 .
While direct evaluation of Eq. (1) becomes available using interpolation techniques for various materials 14,[53][54][55] , it is still computationally expensive. To simplify the computation in EPIC STAR, we again reduce the complex nk-dependent scattering rate to an energy-dependent averaged scattering rate τ À1 E ð Þ. This is exact in the case of an isotropic band dispersion 39 . For anisotropic system, it is still a reasonable approximation which will be verified in this work. In this approximation, the transport distribution function simply becomes In the EPA approximation 37 , the energy-dependent scattering time is obtained by replacing electron-phonon matrix element squared g nml k; q ð Þ j j 2 with its average value at energy pair for ε nk and ε mkþq . This can be achieved by energy-bin-based averaging 37 or moving least squares averaging (MLS-EPA) 56 . This is also well justified for non-polar optical phonon scattering and inter-valley scattering, where the matrix element varies slowly and can be approximated by the corresponding optical or intervalley deformation potential constants 39 . However, we noticed that such treatment may not be good enough in cases where g nml k; q ð Þ j j 2 varies rapidly. This is particularly important for intra-valley acoustic phonon scattering, where g nml k; q ð Þ j j 2 / q j j, and POP scattering where g nml k; q ð Þ j j 2 / q j j À2 . In the first case, while g nml k; q ð Þ j j 2 is linear for long-wavelength acoustic phonons, its Bose-Einstein distribution is not exactly q j j À1 if temperature is not very high. Therefore, separated averaging of g nml k; q ð Þ j j 2 and ω lq may lead to different results in such case. For POP scattering, the scattering matrix diverges as q j j À2 which decays rapidly as q j j increases, so the average value cannot represent the scattering strength.
In accordance with the DPT 30 , both the electron-acousticphonon matrix and the phonon dispersion are linear in q j j around Γ point 38 . Therefore, it can be expected that g nml k; q ð Þ j j 2 =ω lq varies slowly for long wavelength and could be a better candidate for averaging. This also holds for non-polar optical phonon scattering and inter-valley scattering as both g nml k; q ð Þ j j 2 and ω lq vary slowly for these cases. With this in mind, we define a generalized Eliashberg function for short-range electron-phonon interaction 43 Here, the energy and momentum conservations during electron-phonon scattering are preserved through Dirac delta functions, while the positive and negative ω stand for phonon absorption and emission, respectively. The matrix g S nml k; q ð Þ 2 is the remainder of g nml k; q ð Þ j j 2 after removing the long-range contribution g L nml k; q ð Þ 2 47 . By setting E at Fermi energy E F and considering only small ω, it reduces to the conventional isotropic Eliashberg function for metals. Since g nml k; q ð Þ j j 2 =ω lq varies slowly as discussed above, this expression converges fast with increasing Brillouin zone sample density. In the actual first-principles calculations, the delta functions are approximated as Gaussian with broadening of 0.2 eV for electron energy and 0.005 eV for phonon energy, according to our convergence test as given in Supplementary Fig. 1 where converged results were found to be insensitive to smearing parameters. We have also considered the smearing-free tetrahedron method for the integration, however we found that when band edge is not directly sampled by the k-mesh, α 2 l F E; ω ð Þ near band edge is significantly underestimated using tetrahedron integration, as shown in Supplementary Fig. 2. So, it is not well-suited for the target of this study and we use the smearing-based method.
The evaluation of α 2 l F E; ω ð Þ is very fast after obtaining electron-phonon matrix elements using density functional perturbation theory. The computational cost roughly scales as N EPIC $ 3N at N q N k N 2 bnd N E N ω , where N at , N q , N k , N bnd , N E , and N ω are the numbers of atoms, irreducible q points, irreducible k points, relevant bands, energy sample points, and frequency sample points. N E N ω is typically in the order of 10 5 . For comparison, EPW interpolation roughly scales as 3N at N q N WF N k ð Þ 2 N fine q N fine k N bnd where N WF is the number of Wannier functions and N fine q N fine k is the total number of (k, q) pairs to be interpolated. As shown in previous works, N fine q N fine k can be as large as over 10 10 for silicon 9 and even higher for GaAs [14][15][16] . Comparison of computational costs by EPIC STAR and EPW are given in Supplementary Information.
The energy-dependent scattering rate from phonon branch l simply becomes For polar semiconductors like III-V compounds, the electron-POP coupling matrix diverges around Γ point, as described by Fröhlich model 40,57 . The analytical ab initio description based on Born effective charge, dielectric tensor, and phonon dispersion from DFPT calculations has been proposed recently 47,48 . Such strong electron-POP interaction often dominates the scattering as demonstrated using state-of-the-art electron-phonon coupling method for GaAs [14][15][16]48 , anatase TiO 2 47 , two-dimensional transition metal chalcogenides 58 , etc.
The piezoelectric scattering from polar acoustic phonons usually contributes at low temperatures 39 , so we will focus on POP scattering here. Such electron-POP coupling matrix divergence significantly increases the difficulty of scattering rate evaluation as very dense q mesh needs to be sampled [14][15][16] . Therefore, we take advantage of the analytical expression of such POP scattering matrix elements to simplify the computation. Here, we take the average electron-POP coupling matrix iC POP;l q j j for each phonon branch l where C 2 POP;l is the angular averaged Fröhlich coupling strength given by an integration over a solid angle Here g L l;q is the ab initio long-range contribution to the intraband electron-phonon interaction as defined in ref. 47 and q is unit vector along q. For acoustic phonons, this term vanishes due to translational symmetry. For transverse optical phonons, it also vanishes as polarization is orthogonal to momentum q. This leaves only LO branches with non-zero C 2 POP;l . The POP contribution from phonon branch l is then 39 where m b E ð Þ is the average energy-dependent effective mass from band curvature calculation with the energy E relative to the band edge, ω l is the average frequency of phonon branch l near Γ point, while μ and T are chemical potential and absolute temperature, respectively. After subtracting iC POP;l q j j 2 from the squared electron-phonon coupling matrix g nml k; q ð Þ j j 2 , we use the remainder to evaluate α 2 l F E; ω ð Þ and corresponding scattering rate τ À1 l . All scattering rate contributions will be added together to calculate the scattering time τ E ð Þ. This information including C POP;l and α 2 l F E; ω ð Þ may also be used to construct a full Boltzmann equation for direct iterative solution 9,14 , but we will focus on scattering time approximation in this work.
In the presence of free carriers, electrostatic screening effect could be important 39 which reduces the impurity and POP scattering strength. Taking this effect into account, we include it as screening length L which is defined as 38,39 The relative dielectric constant ϵ is obtained from the same DFPT calculation for phonons 59 . It will be used in the dielectric function ϵ q ð Þ ¼ 1 þ L À2 q j j À2 . For heavily doped semiconductors, it is also necessary to consider charged impurities induced by doping. Here we consider the simple yet powerful Brooks-Herring 39,49,50 model of charged impurity scattering with free carrier screening. The corresponding scattering rate is Validations and discussions Silicon is an extensively studied and widely used high mobility semiconductor. Its room temperature electron intrinsic mobility from experimental measurements ranges from about 1400 to 1800 cm 2 V −1 s −1 60-63 . The phonon-limited electron mobility of silicon has been modeled using BTE with electron-phonon scattering from DPT 30 , linear interpolation 13 , and Wannier interpolation 9,14,64 calculations. These full electron-phonon coupling calculations reported intrinsic electron mobility from 1366 to 1915 cm 2 V −1 s −1 , depending on the accuracy of band structure, lattice constant, and Brillouin zone sampling density. The temperature dependence of mobility is also accurately reproduced by phonon-limited Boltzmann transport for a wide range of temperatures, while the effect of doping can also be well described with inclusion of impurity scattering in Brooks-Herring model 9 . The general agreement between experimental measurement and theoretical prediction demonstrates the predictive power of BTE. Therefore, it would be useful to compare the carrier mobility of silicon calculated using EPIC STAR with those in the literature for validation of this fast and simplified method. First, we calculated the phonon-limited scattering time for silicon using both EPIC STAR and MLS-EPA, and compared them with results from literature as shown in Fig. 1a, b. As expected, EPIC STAR shows good agreement especially for low energy carriers, where optical phonon emission is forbidden by energy conservation and acoustic phonon scattering is important. This accordance is attributed to the explicit ω-dependence of α 2 l F E; ω ð Þ which should be more accurate for the rapidly varying acoustic phonon frequency, instead of a single ω for the whole acoustic phonon branch in MLS-EPA. For carriers with higher energy, both EPIC STAR and MLS-EPA agree well with results from full electron-phonon coupling calculations using EPW and from other published works 13,14 .
Then, we calculated the intrinsic electron (hole) mobility defined as μ eðhÞ ¼ σ=en eðhÞ , with n e h ð Þ being n-(p-)type carrier concentration which is lower than 10 14 cm −3 , such that impurity effect is negligible. The calculated intrinsic mobilities are compared with experimental measurement for high purity silicon with temperatures from 200 to 500 K [60][61][62][63]65 , and theoretical predictions using phonon-limited BTE as shown in Fig. 1c,  d 9,13,14,64 . The electron mobility calculated using EPIC STAR is in  67 are illustrated in gray. Doping dependences for e n-type electron mobility and f p-type hole mobility in lightly to heavily doped silicon are compared with experimental measurement and theoretical predictions. Upward-pointing triangles are from ref. 68 while downwardpointing triangles are from ref. 69 . Dashed lines are first-principles predictions 9,64 . In heavily doped region, a better agreement is achieved by including ionized impurity scattering (el−ph+imp), which clearly shows the important role of impurities in silicon when doping concentration is increased. Phonon-limited BTE prediction in literatures produced electron mobility from 1366 to 1915 cm 2 V −1 s −1 at 300 K, while EPIC STAR predicts it to be 1637 cm 2 V −1 s −1 . The hole mobility is overestimated in all theoretical predictions as compared with experimental values, which could be ascribed mainly to the significant underestimation of heavy hole effective mass. The heavy hole effective mass along [100] direction extracted from experiments is 0.46m e in ref. 66 and 0.55m e in ref. 65 while firstprinciples predictions, including both LDA/GGA and GW calculations 9,14 , are only half of the experimental values. Detailed analyses of the role of density functionals, GW corrections, spin-orbit couplings and electron-phonon renormalization are presented in ref. 9 where the hole mobility at 300 K using band structure from first-principles varies from 658 to 820 cm 2 V −1 s −1 .
In that work, a better agreement was achieved using experimental effective mass, resulting in hole mobility of 502 cm 2 V −1 s −1 as compared with experimental value of 450 cm 2 V −1 s −1 67 . EPIC STAR calculation using LDA without spin-orbit coupling predicts hole mobility of 840 cm 2 V −1 s −1 , achieving agreement with full electron-phonon coupling calculation at a much lower computational cost. And the difference with experimental value mainly comes from hole effective mass and dielectric constant discrepancy between DFT band structure and experimental value 9 . By using experimental hole effective mass, we obtain a hole mobility of 660 cm 2 V −1 s −1 at 300 K, and after rescaling the electron-phonon coupling using experimental dielectric constant following ref. 9 it further reduces to 520 cm 2 V −1 s −1 which is very close to the experimental value. This observation suggests that an accurate DFT band structure is important for reliably predicting the charge transport properties. The carrier mobilities in heavily doped cases are also investigated, as shown in Fig. 1e, f. At high doping level, inclusion of impurity scattering is necessary to reproduce experimental trend 68,69 where mobility decreases significantly as more dopants are introduced.
The full electron-phonon BTE calculations require very dense sampling in both electron and phonon Brillouin zone. Previous works using linear interpolation 13 and Wannier interpolation 14 show that converged results are achieved using uniform k-and qgrids with 10 6 points in the electron and phonon Brillouin zone. With symmetry considerations and random sampling, a smaller number of inequivalent points, i.e., 8.5 × 10 4 k-points and 2 × 10 5 q-points would be sufficient 9 . However, this is still a large number as every (k, q) pair needs to be sampled. To estimate the computational cost, we investigated the total number of CPU hours consumed by EPIC STAR and EPW respectively, using the same computer with 2 Intel Xeon E5-2690 v3 CPUs. DFPT calculation using QUANTUM ESPRESSO, which is needed for both EPIC STAR and EPW, took 6 CPU hours with 6 3 q-points. The evaluation of electron-phonon matrix in Wannier basis took 112 CPU hours. And a subsequent Wannier interpolation from 6 3 qpoints and 12 3 k-points to 40 3 q-mesh and 3322 irreducible kpoints (reduced from 40 3 uniform grid) using EPW 8 was performed for speed test. With very small interpolated sampling in Brillouin zone, it still took 322 CPU hours which is much longer than the DFPT calculation itself. With a denser 8.5 × 10 4 k-points and 2 × 10 5 q-points, the hours could be further increased by orders of magnitude. For comparison, EPIC STAR calculation based on the same DFPT calculation took less than 2 CPU hours to obtain α 2 l F, ω l , v 2 αβ E ð Þ, and ρ E ð Þ, and the subsequent scattering time and BTE calculations required only a few seconds. So, the total CPU time used for EPW calculation on a small fine grid, which is 440 CPU hours, is much longer than the total time of 8 CPU hours for EPIC STAR calculation. And the difference would be even larger to reach convergence for EPW. In addition, if the Wannierization is not done correctly, interpolation may result in dissimilar band structures especially in cases with band entanglements 70,71 , so it may not be trivial to automate the process for arbitrary materials. EPIC STAR, on the other hand, uses Fourier interpolation as implemented in BoltzTraP which does not require such information, and can be interpolated from a much denser initial k-grid at low cost. In combination with order-of-magnitude lower computational cost, EPIC STAR enables fast and reliable prediction of phonon-and impurity-limited charge transport properties.
Polar optical phonon scattering, which is important in polar materials such as III-V compounds like GaAs [14][15][16]48 , is also treated in EPIC STAR. Due to the diverging behavior of electron-POP coupling matrix elements, direct Brillouin zone integration methods like EPW requires an even denser q-point sampling for such materials as compared to non-polar materials. A uniform qgrid of 10 6 -10 7 points is needed for electron scattering rate in GaAs, while extremely fine k-points over 10 7 are needed for determining charge transport properties like carrier mobility [14][15][16] . First-principles transport properties prediction of polar semiconductors is therefore difficult. Via EPIC STAR, we are able to provide a workaround by calculating the energy-dependent scattering rates using analytical expression. DFPT and EPIC STAR calculation using a coarse q-grid of 6 3 takes 50 CPU hours. Using EPW to interpolate from the same coarse grid to 80 3 q-points and 23,842 irreducible k-points (reduced from an 80 3 uniform mesh), it took 767 CPU hours which is longer by over 1 order of magnitude. Since full EPW calculation for GaAs requires much more k-and qpoints 14 , the actual difference should be even larger.
As shown in Fig. 2a, while EPIC STAR gives quantitative agreement with full electron-phonon calculations in literature [14][15][16] , the MLS-EPA results deviate significantly from them. This is a natural consequence of the absence of POP scattering, which dominates the electron scattering in GaAs and is greater than the scattering from acoustic and non-polar optical phonons by almost an order of magnitude near the band edge. While the experimentally measured electron mobility varies from 0.7 × 10 4 to 1 × 10 4 cm 2 V −1 s −1 at 300 K 72-76 , their theoretical results range from 0.7 × 10 4 to about 1.2 × 10 4 cm 2 V −1 s −1 14-16 , with energy scattering time approximation results being lower than the iterative solution of BTE. Thus, EPIC STAR prediction of 9770 cm 2 V −1 s −1 , at a much smaller cost, shows general agreement with both experimental measurements and previous theoretical works which used accurate yet more expensive methods by direct integration of the Brillouin zone.
In addition to carrier mobility, EPIC STAR can also predict other transport coefficients from BTE, including electrical conductivity, Seebeck coefficient, and electronic contribution to thermal conductivity. Here, we use EPIC STAR to study the thermoelectric performance of Mg 2 Si. Mg 2 X (X = Si, Ge, Sn) is a family of high performance thermoelectric materials which do not comprise expensive and less abundant elements like Bi or Te, or toxic elements like Pb 77 . Simple Mg 2 X such as n-type Mg 2 Si is able to achieve zT of 0.86 78  The electron mobilities of doped Mg 2 Si were first studied and compared with those in literature 80 . The calculation was performed with the same doping concentration as the experiments, from which good agreement has been achieved between EPIC STAR and experimental measurement for a wide temperature range as shown in Fig. 3a.
Bi-doped Mg 2 Si with different doping levels have been experimentally studied with Bi:Mg 2 Si ratio from 0.1 to 2% 78 . The highest zT of 0.86 was realized at 862 K with 2% Bi doping. Here, we performed EPIC STAR calculations for each doping level studied experimentally and compared the resulting electrical conductivity, Seebeck coefficient, and figure of merit zT with the experimental measurements. The carrier concentrations are fixed to the experimental value at 300 K from Hall measurements 78 .
Since the PBE GGA band gap of 0.22 eV is an underestimation which will affect the predicted value of Seebeck coefficient for low doping levels, we instead used the TB-mBJ meta-GGA 81 band gap of 0.6 eV from ref. 82 by performing a rigid shift of conduction band. The lattice contribution to thermal conductivity κ L , which is beyond the scope of this study, is taken from experimental data which was estimated based on Wiedemann-Franz law where Lorenz number was assumed to be 2.45 × 10 −8 V 2 K −2 78 . However, the actual Lorenz number may differ from this value and depends on both doping and temperature 37 . Therefore, the figure of merit zT here is just an estimation. The calculated zT, electrical conductivity, and Seebeck coefficient are compared with experimental measurements in ref. 83 as shown in Fig. 3b-d. The predicted intrinsic mobility and thermoelectric performance agrees with experimental measurements along a wide temperature range. This demonstrates the predictive power of our approach for both intrinsic and doped semiconductors.
Half-Heusler (HH) alloys constitute an important class of high power factor TE materials 84,85 . Among many HH compounds, NbFeSb has been proposed as one of the promising candidates in  a theoretical study in 2008 86 using BTE in CSTA. Later in 2014 it has been confirmed as a good TE material with zT~1 at 700°C in a combined experimental and theoretical study using EPA 2 . A record-high room temperature power factor of 106 μWcm À1 K À2 was achieved with Ti-doped NbFeSb, and it has been suggested that electron-phonon interaction is the dominant charge carrier scattering mechanism 87 . A very recent theoretical study using accurate EPW method successfully reproduced NbFeSb TE performance for a wide range of temperatures, and indicates that the high power factor of HH alloys mainly arises from weak electron-acoustic-phonon scattering which comes from the nonbonding feature of electron orbitals 38 . The extensive experimental and theoretical study of NbFeSb provides another good candidate for EPIC STAR validation, and our approach indeed achieves excellent agreement with both experimental measurement and EPW prediction from 300 to 1000 K. As shown in Fig. 4, EPIC STAR predicts a power factor of 99 μWcm À1 K À2 , comparable to experimentally measured 106 μWcm À1 K À2 and EPW prediction of around 90 μWcm À1 K À2 38 . The Seebeck coefficient is slightly higher than literature while conductivity is lower. This EPIC STAR calculation used a coarse 6 3 q-mesh, which is the same as that in EPW before Wannier interpolation. With EPW, interpolation to 40 3 k-and q-points needs 497 CPU hours, while 80 3 is needed to reach converged results 38 , which requires at least thousands of CPU hours based on previous estimation. However, DFPT calculation took 145 CPU hours and the subsequent EPIC-STAR calculation only took less than 5 CPU hours to obtain comparable results, reducing time by at least 1 order of magnitude. And the intrinsic hole mobility from EPIC STAR calculation, which is 64 cm 2 V −1 s −1 , is comparable to 77 cm 2 V −1 s −1 from EPW calculation 38 while EPA gives a significant overestimation of 149 cm 2 V −1 s −1 which comes from the absence of POP scattering. The results from EPA on p-type NbFeSb is given in Supplementary Information. This again emphasizes the efficiency and reliability of EPIC STAR which is expected to be generally applicable for a broader range of materials.
In addition to the validation of the above theoretical and experimental results in the literature, we also apply this method to predict the thermoelectric performance of NaInSe 2 which has not been reported yet. It has been synthesized experimentally whereby a layered anisotropic crystal structure was found, as shown in Supplementary Fig. 2 88 . It possesses a similar crystal structure to NaCoO 2 , of which its thermal conductivity is greatly suppressed due to rattling modes in the presence of Na vacancy 89 . Although Na ions can be rattling in the presence of vacancy, its ionic contribution to charge transport is very small in similar structures 90 , so electronic transport is expected to dominate the electrical properties. Reduced thermal conductivity is desirable for thermoelectric applications, and therefore, its electronic part of thermoelectric performance is explored in this work.
First, we performed EPW calculation to validate our calculation of scattering times. The case for NaInSe 2 is complicated because of its complex electronic band structure and phonon dispersion as shown in Fig. 5a, b. The lowest conduction valley at Γ point is energetically very close to the second lowest triply degenerate F valley which is only slightly more positive by 0.1 eV, and they possess very different effective masses. The valence band has a strong effective mass anisotropy which is ten times heavier along the z direction than the in-plane directions (see Supplementary  Fig. 2 for the definitions). The phonon dispersion is also different along in-plane and out-of-plane directions as shown in Fig. 5b. The electron-phonon scattering as illustrated in Fig. 5c is dominated by polar longitudinal optical phonons even for high energy electrons at 0.5 eV above conduction band edge. This again emphasizes the importance of correct treatment for POP scattering in electron-phonon scattering time prediction. By comparing EPIC STAR scattering time with averaged EPW scattering time for carriers with the same energy, we found good agreement even for such a complicated situation, as shown in Fig. 5d.
With such a complex electronic band structure, the Fermi  24 should be large and be correlated with good thermoelectric performance. Here, we use a slightly modified definition for anisotropic material by distinguishing the in-plane (m Ã k ) and out-of-plane transport effective mass (m Ã ? ). The N Ã v K Ã value for ntype NaInSe 2 is found to be 2.8 for in-plane and 1.8 for out-ofplane transport, while for p-type it is 7.6 for in-plane and 0.7 for out-of-plane. The peak n-type power factor predicted using EPIC STAR at 300 K is PF k ¼ 13:4 μWcm À1 K À2 and PF ? ¼ 8:5 μWcm À1 K À2 , while for p-type PF k ¼ 28:4 μWcm À1 K À2 and PF ? ¼ 2:4 μWcm À1 K À2 , respectively. At low temperature, ptype NaInSe 2 achieves higher PF due to large power factor from complex band structure. However, as temperature increases, ntype performance increases rapidly and outperforms p-type after 600 K. This is because of increased effective valley degeneracy due to electrons in F-valleys which have a 3-fold degeneracy. F-valleys are 0.1 eV higher than the single Γ-valley so they are important only when temperature is high enough. The relative trend of PF is in line with N Ã v K Ã which is expected as the high anisotropy and valley degeneracy lead to enhanced Seebeck coefficients. However, the CSTA result using a single scattering time τ CSTA ¼ 4:5 fs, which albeit fits the peak p-type PF k , results in a serious underestimation of a factor of 3 in n-type power factors. Moreover, the optimal p-doping level from CSTA prediction is also different from EPIC STAR result. It is also observed that CSTA PF increases monotonically as temperature rises, as shown in Supplementary  Fig. 3, which differs significantly from EPIC STAR calculation where PF decreases at high temperature due to stronger electron-phonon scattering (Fig. 5f). This again emphasizes the importance of reliable calculation of scattering times.
In conclusion, we proposed a swift and automation-friendly approach for intrinsic and impurity-limited charge transport property prediction from first-principles. By introducing generalized Eliashberg function and adding polar optical phonon contribution, impurity scattering, and free carrier screening, the new approach achieves high fidelity especially for both non-polar and polar semiconductors with very small computational cost after DFPT calculations. We demonstrated the importance of polar optical phonon scattering, which suggests that care should be taken when the electronic properties of polar semiconductors are studied without polar optical phonon scattering. We verified our approach by comparing with previous experimental and theoretical results for Si, GaAs, Mg 2 Si, and NbFeSb, and also revealed NaInSe 2 as a potential new TE material. As high-throughput DFPT computations have been demonstrated recently 91,92 , our methodology can be widely applied for reliable high-throughput screening of high mobility semiconductors and highperformance thermoelectric and photovoltaic materials.

Transport coefficients from Boltzmann transport equation
Semi-classical Boltzmann transport equation (BTE) could be approximately solved in scattering time approximation, leading to the electronic transport integral tensor of order p as 12 Here, for each state nk, v α nk is the α-direction component of group velocity, τ nk is the scattering time, ε nk is the energy, and f 0 E; μ; T ð Þis the Fermi-Dirac distribution with given chemical potential μ and temperature T. Spin degeneracy is included in the band index n. Ω BZ ¼ 2π ð Þ 3 Ω is the Brillouin zone volume where Ω is the unit cell volume.
The electrical conductivity tensor σ, Seebeck coefficient S, and electronic contribution to thermal conductivity κ e can be evaluated using K p of different orders through 12,93 κ e ¼ 1 TΩ Here, e is the elementary charge such that electron charge is −e. The carrier mobility can be defined simply as μ e h ð Þ ¼ σ=en e h ð Þ with n e h ð Þ being electron (hole) concentration. The only system-dependent ingredient in this formulation is the so-called transport distribution function 93 which Fig. 5 Band structures and properties of NaInSe 2 . a Electronic band structure of NaInSe 2 . b Phonon dispersion of NaInSe 2 . The observed longitudinal and transverse optical phonon (LO-TO) splitting and non-analytical behavior demonstrates the polar semiconductor nature of NaInSe 2 . c Mode-resolved electron-phonon scattering rate contribution from longitudinal and transverse acoustic (LA + TA), longitudinal optical (LO), and transverse optical (TO) phonons. The polar LO branches contribute more than 80% of the scattering rate in a wide energy range. d Electron and hole scattering time calculated using EPW (symbols) and their average (black lines), MLS-EPA (blue lines), and EPIC STAR (red lines). The EPIC STAR is in agreement with averaged EPW results from the band edge to the high energy region while the EPA results show larger discrepancy mainly due to the absence of POP scattering. e The predicted doping-dependent power factors at 300 K compared with CSTA results with τ CSTA ¼ 4:5 fs. A single scattering time (in the case for CSTA) is clearly not enough for the prediction of both electron and hole, and the peak power factor position is also different from phonon-limited study. f The predicted temperature-dependent power factor at fixed hole (2 × 10 20 cm −3 ) and electron (1 × 10 20 cm −3 ) concentration. p-Type power factor is higher at relatively low temperatures, while ntype power factor peaks around 700 K.