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

Abstract

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₂¹, NbFeSb², NbCoSn³, and CdIn₂Te₄⁴, to name a few. High mobility oxides⁵ 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₂BiTaO₆⁶ and TaIrGe⁷, have been experimentally verified. Many of the previous works are based on the prediction of carrier mobility $\mu$, electrical conductivity $\sigma$, and Seebeck coefficient $S$. Although the state-of-the-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 $\mu$, $\sigma$, and $S$, BTE solution within scattering time approximation (STA) is commonly adopted¹² which yields good agreement with fully numerical solution of the equation¹³. 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⁵–10⁶ 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¹³ or Wannier^8,17 interpolation techniques. Therefore, constant scattering time approximation (CSTA) is commonly employed in tools such as BoltzTraP^18,19 and BoltzWann²⁰ where only electronic structure calculation from density functional theory (DFT)²¹ 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²⁴ and electronic fitness function²⁵ 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³⁰, 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 high-throughput search for new TE compounds⁴. The electron–phonon averaged (EPA) approximation³⁷ 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⁸ 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³⁸. Therefore the intrinsic carrier mobility may be overestimated. In fact, we will show that EPA predicts an intrinsic hole mobility of 149 cm² V⁻¹ s⁻¹ for NbFeSb which is twice the value from EPW calculation³⁸ (77 cm² V⁻¹ s⁻¹). 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³⁹. 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⁴⁰, was described by Eliashberg electron–phonon spectral function $\alpha ^2F(\omega )$ which was used in Allen–Dynes formula for transition temperature^41,42 and Ziman’s resistivity formula^12,43. Its dependence on phonon frequency $\omega$ could reflect some details that are missing in DPT and EPA methods, and first-principles calculation of $\alpha ^2F(\omega )$ has been demonstrated for systems like Nb^44,45, Al¹³, Pb^8,45, and B-doped diamond⁴⁶. 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⁴⁷. 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 long-range 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, $\alpha ^2F(E,\omega )$ from first-principles calculations which depends on both electron energy $E$ and phonon frequency $\omega$. 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. $\alpha ^2F(E,\omega )$ 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₂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₂ 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⁴⁰

$$\tau _{n{\mathbf{k}}}^{ - 1} = \mathop {\sum }\limits_l \tau _{ln{\mathbf{k}}}^{ - 1} = \mathop {\sum }\limits_{lm} \frac{{2\pi }}{\hbar } {\int \nolimits_{BZ}} \left| {g_{nml}\left( {{\mathbf{k}},{\mathbf{q}}} \right)} \right|^2\left[ {\left( {1 - f_{m{\mathbf{k}} + {\mathbf{q}}}^0 + n_{l{\mathbf{q}}}^0} \right)\delta \left( {\varepsilon _{n{\mathbf{k}}} - \hbar \omega _{l{\mathbf{q}}} - \varepsilon _{m{\mathbf{k}} + {\mathbf{q}}}} \right)} \right. + \left. {\left( {f_{m{\mathbf{k}} + {\mathbf{q}}}^0 + n_{l{\mathbf{q}}}^0} \right)\delta \left( {\varepsilon _{n{\mathbf{k}}} + \hbar \omega _{l{\mathbf{q}}} - \varepsilon _{m{\mathbf{k}} + {\mathbf{q}}}} \right)} \right]\frac{{d{\mathbf{q}}}}{{{\mathrm{\Omega }}_{{\mathrm{BZ}}}}}$$

(1)

where $g_{nml}\left( {{\mathbf{k}},{\mathbf{q}}} \right)$ is the electron–phonon coupling matrix element among electrons in state $n{\mathbf{k}}$, state $m{\mathbf{k}} + {\mathbf{q}}$, and phonon in branch $l$ with momentum ${\mathbf{q}}$. $n_{l{\mathbf{q}}}^0$ is the Bose–Einstein distribution function of phonons with frequency $\omega _{l{\mathbf{q}}}$ while $f_{m{\mathbf{k}} + {\mathbf{q}}}^0$ is the Fermi–Dirac distribution function of electrons with energy $\varepsilon _{m{\mathbf{k}} + {\mathbf{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 $n{\mathbf{k}}$-dependent scattering rate to an energy-dependent averaged scattering rate $\tau ^{ - 1}\left( E \right)$. This is exact in the case of an isotropic band dispersion³⁹. 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

$$\left[ {{\mathbf{\Sigma }}\left( E \right)} \right]_{\alpha \beta } \approx \overline {v_{\alpha \beta }^2} \left( E \right)\rho \left( E \right)\tau \left( E \right)$$

(2)

where the density of states $\rho \left( E \right)$, averaged squared velocity tensor $\overline {v_{\alpha \beta }^2}$ and averaged scattering rate $\tau ^{ - 1}$ can be calculated from band structure as^18,37

$$\rho \left( E \right) = \mathop {\sum }\limits_n {\int \nolimits_{BZ}} \delta \left( {E - \varepsilon _{n{\mathbf{k}}}} \right)\frac{{d{\mathbf{k}}}}{{{\mathrm{\Omega }}_{BZ}}}$$

(3)

$$\overline {v_{\alpha \beta }^2} \left( E \right) = \frac{1}{{\rho \left( E \right)}}\mathop {\sum }\limits_n{\int\nolimits_{BZ}} v_{n{\mathbf{k}}}^\alpha v_{n{\mathbf{k}}}^\beta \delta \left( {E - \varepsilon _{n{\mathbf{k}}}} \right)\frac{{d{\mathbf{k}}}}{{{\mathrm{\Omega }}_{BZ}}}$$

(4)

$$\tau ^{ - 1}\left( E \right) = \mathop {\sum }\limits_l \tau _l^{ - 1}\left( E \right) = \mathop {\sum }\limits_l \frac{1}{{\rho \left( E \right)}}\mathop {\sum }\limits_n {\int\nolimits_{BZ}} \tau _{ln{\mathbf{k}}}^{ - 1}\delta \left( {E - \varepsilon _{n{\mathbf{k}}}} \right)\frac{{d{\mathbf{k}}}}{{{\mathrm{\Omega }}_{BZ}}}$$

(5)

In the EPA approximation³⁷, the energy-dependent scattering time is obtained by replacing electron–phonon matrix element squared $\left| {g_{nml }\left( {{\mathbf{k}},{\mathbf{q}}} \right)} \right|^2$ with its average value at energy pair for $\varepsilon _{n{\mathbf{k}}}$ and $\varepsilon _{m{\mathbf{k}} + {\mathbf{q}}}$. This can be achieved by energy-bin-based averaging³⁷ or moving least squares averaging (MLS-EPA)⁵⁶. 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³⁹. However, we noticed that such treatment may not be good enough in cases where $\left| {g_{nml}\left( {{\mathbf{k}},{\mathbf{q}}} \right)} \right|^2$ varies rapidly. This is particularly important for intra-valley acoustic phonon scattering, where $\left| {g_{nml}\left( {{\mathbf{k}},{\mathbf{q}}} \right)} \right|^2 \propto \left| {\mathbf{q}} \right|$, and POP scattering where $\left| {g_{nml}\left( {{\mathbf{k}},{\mathbf{q}}} \right)} \right|^2 \propto \left| {\mathbf{q}} \right|^{ - 2}$. In the first case, while $\left| {g_{nml}\left( {{\mathbf{k}},{\mathbf{q}}} \right)} \right|^2$ is linear for long-wavelength acoustic phonons, its Bose–Einstein distribution is not exactly $\left| {\mathbf{q}} \right|^{ - 1}$ if temperature is not very high. Therefore, separated averaging of $\left| {g_{nml}\left( {{\mathbf{k}},{\mathbf{q}}} \right)} \right|^2$ and $\omega _{l{\mathbf{q}}}$ may lead to different results in such case. For POP scattering, the scattering matrix diverges as $\left| {\mathbf{q}} \right|^{ - 2}$ which decays rapidly as $\left| {\mathbf{q}} \right|$ increases, so the average value cannot represent the scattering strength.

In accordance with the DPT³⁰, both the electron–acoustic-phonon matrix and the phonon dispersion are linear in $\left| {\mathbf{q}} \right|$ around Γ point³⁸. Therefore, it can be expected that $\left| {g_{nml}\left( {{\mathbf{k}},{\mathbf{q}}} \right)} \right|^2/\omega _{l{\mathbf{q}}}$ 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 $\left| {g_{nml }\left( {{\mathbf{k}},{\mathbf{q}}} \right)} \right|^2$ and $\omega _{l {\mathbf{q}}}$ vary slowly for these cases. With this in mind, we define a generalized Eliashberg function for short-range electron–phonon interaction⁴³

$$\begin{array}{l} \alpha _l^2F\left( {E,\omega } \right) = \frac{\omega }{{\rho \left( E \right)}}\mathop {\sum }\limits_{mn} {\int\nolimits_{BZ}} \frac{{d{\mathbf{k}}}}{{{\mathrm{\Omega }}_{BZ}}}\frac{{d{\mathbf{q}}}}{{{\mathrm{\Omega }}_{BZ}}}\frac{{\left| {g_{nml}^S\left( {{\mathbf{k}},{\mathbf{q}}} \right)} \right|^2}}{{\omega _{l{\mathbf{q}}}}}\delta \left( {E - \varepsilon _{n{\mathbf{k}}}} \right)\\\times\,\delta \left( {\hbar \omega + E - \varepsilon _{m{\mathbf{k}} + {\mathbf{q}}}} \right)\left[ {\delta \left( {\omega - \omega _{l{\mathbf{q}}}} \right) + \delta \left( {\omega + \omega _{l{\mathbf{q}}}} \right)} \right]\end{array}$$

(6)

Here, the energy and momentum conservations during electron–phonon scattering are preserved through Dirac delta functions, while the positive and negative $\omega$ stand for phonon absorption and emission, respectively. The matrix $\left| {g_{nml}^S\left( {{\mathbf{k}},{\mathbf{q}}} \right)} \right|^2$ is the remainder of $\left| {g_{nml }\left( {{\mathbf{k}},{\mathbf{q}}} \right)} \right|^2$ after removing the long-range contribution $\left| {g_{nml}^L\left( {{\mathbf{k}},{\mathbf{q}}} \right)} \right|^2$ ⁴⁷. By setting $E$ at Fermi energy $E_F$ and considering only small $\omega$, it reduces to the conventional isotropic Eliashberg function for metals. Since $\left| {g_{nml}\left( {{\mathbf{k}},{\mathbf{q}}} \right)} \right|^2/\omega _{l{\mathbf{q}}}$ 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, $\alpha _l^2F\left( {E,\omega } \right)$ 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 $\alpha _l^2F\left( {E,\omega } \right)$ is very fast after obtaining electron–phonon matrix elements using density functional perturbation theory. The computational cost roughly scales as $N_{EPIC}\sim 3N_{at}N_qN_kN_{bnd}^2N_EN_\omega$, where $N_{at}$, $N_q$, $N_k$, $N_{bnd}$, $N_E$, and $N_\omega$ are the numbers of atoms, irreducible q points, irreducible k points, relevant bands, energy sample points, and frequency sample points. $N_EN_\omega$ is typically in the order of 10⁵. For comparison, EPW interpolation roughly scales as $3N_{at}N_q\left( {N_{WF}N_k} \right)^2N_q^{fine}N_k^{fine}N_{bnd}$ where $N_{WF}$ is the number of Wannier functions and $N_q^{fine}N_k^{fine}$ is the total number of (k, q) pairs to be interpolated. As shown in previous works, $N_q^{fine}N_k^{fine}$ can be as large as over 10¹⁰ for silicon⁹ 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

$$\tau _l^{ - 1}\left( {E,\mu ,T} \right) = 2\pi \left\{ {{\int \nolimits_0^{+ \infty}} \alpha _l^2F\left( {E,\omega } \right)\left[ {f\left( {E + \hbar \omega - \mu ,T} \right) + n\left( {\hbar \omega ,T} \right)} \right]d\omega + {\int \nolimits_{- \infty }^0} \alpha _l^2F\left( {E,\omega } \right)\left[ {1 - f\left( {E + \hbar \omega - \mu ,T} \right) + n\left( {\hbar \omega ,T} \right)} \right]d\omega } \right\}$$

(7)

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₂⁴⁷, two-dimensional transition metal chalcogenides⁵⁸, etc. The piezoelectric scattering from polar acoustic phonons usually contributes at low temperatures³⁹, 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 ${\mathbf{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 $\frac{{iC_{{\mathrm{POP}},l}}}{{\left| {\mathbf{q}} \right|}}$ for each phonon branch $l$ where $C_{{\mathrm{POP}},l}^2$ is the angular averaged Fröhlich coupling strength given by an integration over a solid angle

$$C_{{\mathrm{POP}},l}^2 = \mathop {{\lim }}\limits_{{\mathrm{q}} \to 0} {\int} {\left| {\mathbf{q}} \right|^2\left| {g_{l,{\mathbf{q}}}^L} \right|^2\frac{{d\hat q}}{{4\pi }}}$$

(8)

Here $g_{l,{\mathbf{q}}}^L$ is the ab initio long-range contribution to the intraband electron–phonon interaction as defined in ref. ⁴⁷ and $\hat q$ is unit vector along ${\mathbf{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_{{\mathrm{POP}},l}^2$. The POP contribution from phonon branch $l$ is then³⁹

$$\begin{array}{l}\tau _{{\rm{POP}},l}^{ - 1}(E) = C_{{\rm{POP}},l}^2\frac{\Omega }{{\pi {\hbar ^2}}}\sqrt {\frac{{{m_b}(E)}}{{2E}}} \left\{ {{{\sinh }^{ - 1}}\sqrt {\frac{E}{{\hbar {\omega _l}}}} \left[ {n(\hbar {\omega _l},\,T) + f(E + \hbar {\omega _l} - \mu ,\,T)} \right] +\, {{\sinh }^{ - 1}}\sqrt {\frac{E}{{\hbar {\omega _l}}} - 1} \left[ {n\left( {\hbar {\omega _{l,}}\,T} \right) + 1 - f\left( {E - \hbar {\omega _l} - \mu ,\,T} \right)} \right]} \right\}\end{array}$$

(9)

where $m_b\left( E \right)$ is the average energy-dependent effective mass from band curvature calculation with the energy $E$ relative to the band edge, $\omega _l$ is the average frequency of phonon branch $l$ near Γ point, while $\mu$ and $T$ are chemical potential and absolute temperature, respectively. After subtracting $\left| {\frac{{iC_{{\mathrm{POP}},l}}}{{\left| {\mathbf{q}} \right|}}} \right|^2$ from the squared electron–phonon coupling matrix $\left| {g_{nml}\left( {{\mathbf{k}},{\mathbf{q}}} \right)} \right|^2$, we use the remainder to evaluate $\alpha _l^2F\left( {E,\omega } \right)$ and corresponding scattering rate $\tau _l^{ - 1}$. All scattering rate contributions will be added together to calculate the scattering time $\tau \left( E \right)$. This information including $C_{{\mathrm{POP}},l}$ and $\alpha _l^2F\left( {E,\omega } \right)$ 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³⁹ 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

$$L^{ - 2} = \frac{{e^2}}{{{\it{\epsilon }}\varepsilon _0}}{\int} {\left( { - \frac{{\partial f}}{{\partial E}}} \right)} \rho \left( E \right)dE$$

(10)

The relative dielectric constant ${\it{\epsilon }}$ is obtained from the same DFPT calculation for phonons⁵⁹. It will be used in the dielectric function ${\it{\epsilon }}\left( {\mathbf{q}} \right) = 1 + L^{ - 2}\left| {\mathbf{q}} \right|^{ - 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

$$\tau _{{\mathrm{imp}}}^{ - 1}\left( E \right) \approx \frac{{n_{{\mathrm{imp}}}e^4}}{{16\sqrt {2m_b\left( E \right)} \pi {\it{\epsilon }}^2\varepsilon _0^2}}\frac{1}{{\sqrt {E^3} }}\left[ {\ln \left( {1 + \frac{{8m_b\left( E \right)EL^2}}{{\hbar ^2}}} \right) - \frac{{8m_b\left( E \right)EL^2}}{{\hbar ^2 + 8m_b\left( E \right)EL^2}}} \right]$$

(11)

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² V⁻¹ s^{−1 60,61,62,63}. The phonon-limited electron mobility of silicon has been modeled using BTE with electron–phonon scattering from DPT³⁰, linear interpolation¹³, and Wannier interpolation^9,14,64 calculations. These full electron–phonon coupling calculations reported intrinsic electron mobility from 1366 to 1915 cm² V⁻¹ s⁻¹, 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⁹. 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 $\omega$-dependence of $\alpha _l^2F\left( {E,\omega } \right)$ which should be more accurate for the rapidly varying acoustic phonon frequency, instead of a single $\omega$ 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.

Fig. 1: Electron and hole scattering times and mobilities of Si.

Energy dependence of scattering time in silicon for a electrons and b holes, calculated using EPIC STAR and MLS-EPA as compared with literature¹⁴ using full electron–phonon calculations. Temperature-dependent intrinsic mobilities of c electrons and d holes, compared with experimental measurements (Exp.) and theoretical predictions (Theo.) in literature. Intrinsic mobilities are calculated with carrier concentration lower than 10¹⁴ cm⁻³ to avoid impurity effect. Experimental results are indicated as symbols with crosses⁶⁰, stars⁶¹, squares⁶², circles⁶³, and triangles⁶⁵, Dashed lines are phonon-limited BTE prediction^9,13,14,64. The region with variation of ±30% of the power law fitting of temperature dependence with T^−2.42 and T^−2.2 for electrons and holes, as defined in ref. ⁶⁷ 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. ⁶⁸ while downward-pointing triangles are from ref. ⁶⁹. 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.

Then, we calculated the intrinsic electron (hole) mobility defined as ${\boldsymbol{\mu}}_{e(h)} = {\boldsymbol{\sigma }}/en_{e( h)}$, with $n_{e\left( h \right)}$ being n-(p-)type carrier concentration which is lower than 10¹⁴ cm⁻³, 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 excellent agreement with both experimental and theoretical results in literature, while MLS-EPA overestimates the mobility. Phonon-limited BTE prediction in literatures produced electron mobility from 1366 to 1915 cm² V⁻¹ s⁻¹ at 300 K, while EPIC STAR predicts it to be 1637 cm² V⁻¹ s⁻¹. 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. ⁶⁶ and 0.55m_e in ref. ⁶⁵ while first-principles 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. ⁹ where the hole mobility at 300 K using band structure from first-principles varies from 658 to 820 cm² V⁻¹ s⁻¹. In that work, a better agreement was achieved using experimental effective mass, resulting in hole mobility of 502 cm² V⁻¹ s⁻¹ as compared with experimental value of 450 cm² V⁻¹ s^−1 67. EPIC STAR calculation using LDA without spin–orbit coupling predicts hole mobility of 840 cm² V⁻¹ s⁻¹, 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⁹. By using experimental hole effective mass, we obtain a hole mobility of 660 cm² V⁻¹ s⁻¹ at 300 K, and after rescaling the electron–phonon coupling using experimental dielectric constant following ref. ⁹ it further reduces to 520 cm² V⁻¹ s⁻¹ 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¹³ and Wannier interpolation¹⁴ show that converged results are achieved using uniform k- and q-grids with 10⁶ 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⁴ k-points and 2 × 10⁵ q-points would be sufficient⁹. 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³ q-points. The evaluation of electron–phonon matrix in Wannier basis took 112 CPU hours. And a subsequent Wannier interpolation from 6³ q-points and 12³ k-points to 40³ q-mesh and 3322 irreducible k-points (reduced from 40³ uniform grid) using EPW⁸ 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⁴ k-points and 2 × 10⁵ 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 $\alpha _l^2F$, $\omega _l$, $v_{\alpha \beta }^2\left( E \right)$, and $\rho \left( E \right)$, 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 q-grid of 10⁶–10⁷ points is needed for electron scattering rate in GaAs, while extremely fine k-points over 10⁷ 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³ takes 50 CPU hours. Using EPW to interpolate from the same coarse grid to 80³ q-points and 23,842 irreducible k-points (reduced from an 80³ 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 q-points¹⁴, 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⁴ to 1 × 10⁴ cm² V⁻¹ s⁻¹ at 300 K^{72,73,74,75,76}, their theoretical results range from 0.7 × 10⁴ to about 1.2 × 10⁴ cm² V⁻¹ s^{−1 14,15,16}, with energy scattering time approximation results being lower than the iterative solution of BTE. Thus, EPIC STAR prediction of 9770 cm² V⁻¹ s⁻¹, 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.

Fig. 2: Electron scattering time and mobility of GaAs.

a Energy dependence of scattering time for electrons in GaAs, calculated using EPIC STAR and MLS-EPA, as compared with refs. ^14,15,16 using full electron–phonon calculation which is much more computationally expensive. EPIC STAR achieves comparable results for both scattering rate and intrinsic mobility compared to full electron–phonon calculation, with at least 1 order of magnitude faster computation speed. Due to the absence of POP scattering, the MLS-EPA results (blue dashed line) deviate substantially from the full electron–phonon calculations. Solid red line is purely first-principles EPIC STAR results. The sharp increase in scattering rate near 0.04 eV marks the onset of POP emission process. b Temperature-dependent intrinsic mobility of electrons in GaAs, compared with experimental measurements (Exp.) and full first-principles theoretical predictions (Theo.) in literature. The experimental results are indicated as crosses⁷², stars⁷³, squares⁷⁴, circles⁷⁵, and triangles⁷⁶. Thin dashed lines are phonon-limited BTE prediction with energy scattering time approximation (STA) and with iterative solution of BTE^14,15,16.

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₂Si. Mg₂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⁷⁷. Simple Mg₂X such as n-type Mg₂Si is able to achieve zT of 0.86⁷⁸, while higher zT has been realized using solid solutions of different X through band convergence. For example, zT of Mg_2.15Si_0.28Sn_0.71Sb_0.006 can reach up to 1.3 at 700 K²², and that of Mg_2.16(Si_0.3Ge_0.05Sn_0.65)_0.98Sb_0.02 with 1.45 at 750 K⁷⁹.

The electron mobilities of doped Mg₂Si were first studied and compared with those in literature⁸⁰. 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.

Fig. 3: Transport properties of Mg₂Si.

a Electron mobility calculated at different temperature, as compared with experimental measurement by Morris et al.⁸⁰. Predicted values for b figure of merit zT, c electrical conductivity, and d Seebeck coefficients of Bi-doped Mg₂Si are compared with previous experimental measurement⁷⁸.

Bi-doped Mg₂Si with different doping levels have been experimentally studied with Bi:Mg₂Si ratio from 0.1 to 2%⁷⁸. 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⁷⁸. 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⁸¹ band gap of 0.6 eV from ref. ⁸² by performing a rigid shift of conduction band. The lattice contribution to thermal conductivity $\kappa _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⁻⁸ V² K^−2 78. However, the actual Lorenz number may differ from this value and depends on both doping and temperature³⁷. 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. ⁸³ 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⁸⁶ 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². A record-high room temperature power factor of 106 ${\mathrm{\upmu Wcm}}^{ - 1}{\mathrm{K}}^{ - 2}$ was achieved with Ti-doped NbFeSb, and it has been suggested that electron–phonon interaction is the dominant charge carrier scattering mechanism⁸⁷. 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 non-bonding feature of electron orbitals³⁸. 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 ${\mathrm{\upmu Wcm}}^{ - 1}{\mathrm{K}}^{ - 2}$, comparable to experimentally measured 106 ${\mathrm{\upmu Wcm}}^{ - 1}{\mathrm{K}}^{ - 2}$ and EPW prediction of around 90 ${\mathrm{\upmu Wcm}}^{ - 1}{\mathrm{K}}^{ - 2}$ ³⁸. The Seebeck coefficient is slightly higher than literature while conductivity is lower. This EPIC STAR calculation used a coarse 6³ q-mesh, which is the same as that in EPW before Wannier interpolation. With EPW, interpolation to 40³ k- and q-points needs 497 CPU hours, while 80³ is needed to reach converged results³⁸, 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² V⁻¹ s⁻¹, is comparable to 77 cm² V⁻¹ s⁻¹ from EPW calculation³⁸ while EPA gives a significant overestimation of 149 cm² V⁻¹ s⁻¹ 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.

Fig. 4: p-type NbFeSb thermoelectric performance.

Thermoelectric performance of p-type NbFeSb, including a power factor, b electrical conductivity, and c Seebeck coefficient from EPIC STAR prediction from 300 to 1000 K, compared with experimental measurement of Ti-doped NbFeSb⁸⁷ and EPW prediction in literature³⁸. Calculations were performed at a fixed hole-doping concentration of 2.0 × 10²⁰ cm⁻³ for all temperatures.

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₂ 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⁸⁸. It possesses a similar crystal structure to NaCoO₂, of which its thermal conductivity is greatly suppressed due to rattling modes in the presence of Na vacancy⁸⁹. Although Na ions can be rattling in the presence of vacancy, its ionic contribution to charge transport is very small in similar structures⁹⁰, 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₂ 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.

Fig. 5: Band structures and properties of NaInSe₂.

a Electronic band structure of NaInSe₂. b Phonon dispersion of NaInSe₂. The observed longitudinal and transverse optical phonon (LO–TO) splitting and non-analytical behavior demonstrates the polar semiconductor nature of NaInSe₂. 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 ${\uptau}_{{\mathrm{CSTA}}} = 4.5\,{\mathrm{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²⁰ cm⁻³) and electron (1 × 10²⁰ cm⁻³) concentration. p-Type power factor is higher at relatively low temperatures, while n-type power factor peaks around 700 K.

With such a complex electronic band structure, the Fermi Surface Complexity Factor $N_v^\ast K^\ast = \left( {\frac{{m_S^\ast }}{{m_c^\ast }}} \right)^{\frac{3}{2}}$ proposed by Gibbs et al.²⁴ 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_\parallel ^\ast$) and out-of-plane transport effective mass ($m_ \bot ^\ast$). The $N_v^\ast K^\ast$ value for n-type NaInSe₂ is found to be 2.8 for in-plane and 1.8 for out-of-plane 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 ${\mathrm{PF}}_\parallel = 13.4\,{\mathrm{\upmu Wcm}}^{ - 1}{\mathrm{K}}^{ - 2}$ and ${\mathrm{PF}}_ \bot = 8.5\,{\mathrm{\upmu Wcm}}^{ - 1}{\mathrm{K}}^{ - 2}$, while for p-type ${\mathrm{PF}}_\parallel = 28.4\,{\mathrm{\upmu Wcm}}^{ - 1}{\mathrm{K}}^{ - 2}$ and ${\mathrm{PF}}_ \bot = 2.4\,{\mathrm{\upmu Wcm}}^{ - 1}{\mathrm{K}}^{ - 2}$, respectively. At low temperature, p-type NaInSe₂ achieves higher PF due to large power factor from complex band structure. However, as temperature increases, n-type 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^\ast K^\ast$ which is expected as the high anisotropy and valley degeneracy lead to enhanced Seebeck coefficients. However, the CSTA result using a single scattering time $\tau _{{\mathrm{CSTA}}} = 4.5\,{\mathrm{fs}}$, which albeit fits the peak p-type ${\mathrm{PF}}_\parallel$, 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₂Si, and NbFeSb, and also revealed NaInSe₂ 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 high-performance thermoelectric and photovoltaic materials.

Methods

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¹²

$$\left[ {{{\mathbf{\cal{K}}}}^p} \right]_{\alpha \beta } = {\int \nolimits_{ - \infty }^{ + \infty }} \left[ {\mathop {\sum }\limits_n{\int \nolimits_{BZ}} v_{n{\mathbf{k}}}^\alpha v_{n{\mathbf{k}}}^\beta \tau _{n{\mathbf{k}}}\delta \left( {E - \varepsilon _{n{\mathbf{k}}}} \right)\frac{{d{\mathbf{k}}}}{{{\mathrm{\Omega }}_{BZ}}}} \right]\left( {E - \mu } \right)^p\left[ { - \frac{{\partial f^0\left( {E,\mu ,T} \right)}}{{\partial E}}} \right]dE$$

(12)

Here, for each state $n{\mathbf{k}}$, $v_{n{\mathbf{k}}}^\alpha$ is the $\alpha$-direction component of group velocity, $\tau _{n{\mathbf{k}}}$ is the scattering time, $\varepsilon _{n{\mathbf{k}}}$ is the energy, and $f^0\left( {E,\mu ,T} \right)$ is the Fermi–Dirac distribution with given chemical potential $\mu$ and temperature $T$. Spin degeneracy is included in the band index $n$. ${\mathrm{\Omega }}_{{\mathrm{BZ}}} = \frac{{\left( {2\pi } \right)^3}}{{\mathrm{\Omega }}}$ is the Brillouin zone volume where ${\mathrm{\Omega }}$ is the unit cell volume.

The electrical conductivity tensor ${\boldsymbol{\sigma }}$, Seebeck coefficient ${\boldsymbol{S}}$, and electronic contribution to thermal conductivity ${\boldsymbol{\kappa }}_e$ can be evaluated using ${\mathbf{\cal{K}}}^p$ of different orders through^12,93

$${\boldsymbol{\sigma }} = \frac{1}{\Omega}e^2{\bf{\cal{K}}}^0$$

(13)

$${\mathbf{S}} = \frac{1}{- eT}\left[ {\mathbf{\cal{K}}}^0 \right]^{ - 1}{\mathbf{\cal{K}}}^1$$

(14)

$${\boldsymbol{\kappa }}_e = \frac{1}{T\Omega}\left[ {{\mathbf{\cal{K}}}^2 - {\mathbf{\cal{K}}}^1\left[ {{\mathbf{\cal{K}}}^0} \right]^{ - 1}{\mathbf{\cal{K}}}^1} \right]$$

(15)

Here, $e$ is the elementary charge such that electron charge is −e. The carrier mobility can be defined simply as ${\boldsymbol{\mu }}_{e\left( h \right)} = {\boldsymbol{\sigma }}/en_{e\left( h \right)}$ with $n_{e\left( h \right)}$ being electron (hole) concentration. The only system-dependent ingredient in this formulation is the so-called transport distribution function⁹³ which projects the complex electronic state $n{\mathbf{k}}$-dependence into simple energy $E$-dependent formula

$$\left[ {{\mathbf{\Sigma }}\left( E \right)} \right]_{\alpha \beta } = \mathop {\sum }\limits_n {\int\nolimits_{BZ}} v_{n{\mathbf{k}}}^\alpha v_{n{\mathbf{k}}}^\beta \tau _{n{\mathbf{k}}}\delta \left( {E - \varepsilon _{n{\mathbf{k}}}} \right)\frac{{d{\mathbf{k}}}}{{{\mathrm{\Omega }}_{BZ}}}$$

(16)

DFT calculations

All first-principles calculations were performed using QUANTUM ESPRESSO^51,52 including structure optimization, ground state calculation, and band structure calculations using density functional theory. Phonon and electron–phonon coupling matrix were obtained from DFPT^59,94 calculation implemented in QUANTUM ESPRESSO. Non-self-consistent-field band structure calculations were performed on a uniform 48³ k-mesh and Fourier interpolation into a much denser 240³ mesh was performed using BoltzTraP in order to obtain the density of states, squared velocity tensor, and energy-dependent effective mass. Local density approximation (LDA)⁹⁵ was used for Si and GaAs for comparison with previous theoretical literature, while generalized-gradient approximation (GGA) with exchange-correlation function given by Perdew, Burke, and Ernzerhof (PBE)⁹⁶ was used for Mg₂Si, NbFeSb, and NaInSe₂. Norm-conserving pseudopotentials^97,98 were used to replace core electrons. Although corrections such as spin–orbit coupling and Hubbard U were not considered in this work, they can be added for all calculations having been implemented and tested for band structure and DFPT calculations in QUANTUM ESPRESSO⁵².