Computationally efficient matrix elements

The scattering rate of an electron from an initial state nk, where n is a band index and k is a wave vector, to final state mk + q is described by Fermi’s golden rule as

Historically, this challenge has been avoided. In the constant relaxation time approximation (CRTA), Eq. (1) is simplified to a single constant. An alternative is to employ model matrix elements formulated for isotropic band structures based on intrinsic materials parameters. For example, the treatment of deformation potential scattering due to long-wavelength acoustic phonons proposed by Bardeen and Shockley⁸ depends only on an averaged elastic constant and band edge deformation potential; it ignores perturbations from transverse phonon modes and anisotropy in the deformation response. This simple approach has been employed widely in computations of acoustic phonon scattering but is unreliable and does not generalise to complex systems or metals^49–51. An alternative approach, developed by Khan and Allen⁴⁹, can reproduce the fully first principles electron–phonon scattering rate if the strain tensor caused by the phonon and an additional velocity term are included. The resulting matrix element is given by

In the present work, we combine the simplicity of the Bardeen-Shockley approach with the accuracy of the Khan-Allen matrix element by exploiting the acoustoelastic properties of materials. The dispersion relations for acoustic waves are contained in the Christoffel equation⁵²

Scattering by acoustic phonons through the piezoelectric interaction (pi) occurs in non-centrosymmetric systems and can dominate at low temperatures (≲50 K). We have applied a similar treatment to extend the isotropic matrix element of Meijer and Polder¹¹, Harrison⁵⁵, and Zook⁵³, to include the full piezoelectric stress tensor h and scattering from all three acoustic modes. The resulting matrix element is given by

We treat polar optical phonon scattering (po) by extending the Frölich model¹² to include quantum mechanical wave function overlaps and anisotropic permittivity. Here, electrons in a dielectric medium are perturbed by a dispersionless longitudinal optical phonon mode with frequency ω_po. Our electron–phonon matrix element takes the form

Following the classic treatment of Brooks and Herring^9,58 we consider the scattering from fully ionized impurities (ii) modeled as screened Coulomb potentials, with the matrix element given by

The final k-dependent scattering rates are obtained by integrating Eq. (1) over all phonon wave vectors (q) in the first Brillouin zone. Elastic scattering processes are well described by the MRTA to the BTE due to the requirement that τ_nk→mk+q = τ_mk+q→nk⁴⁰. As this condition does not hold for inelastic processes, we adopt the self-energy relaxation time approximation (SERTA) to obtain the final polar phonon coupling rates³². Further justification for this approach is detailed in the Supplementary Methods. Electronic eigenvalues and group velocities needed to calculate scattering and transport properties are Fourier interpolated onto dense Brillouin zone grids using the BoltzTraP2 software²² (as detailed in the Supplementary Methods). Electron mobility and Seebeck coefficient are calculated using the linearized BTE via the Onsager transport coefficients^22,61 (Supplementary Eqs. 13 to 15). We also employ a custom procedure for selecting the most important k-points at which to calculate scattering to further reduce the computational expense (detailed in the Supplementary Methods).

Unlike other state-of-the-art approaches in which a computationally expensive DFPT calculation is required to obtain g(k, q), in our method all matrix elements depend only on common materials parameters (ω_po, ϵ_s, ϵ_∞, etc.) that can be calculated relatively inexpensively. Crucially, many of these properties are already tabulated in databases such as the Materials Project⁶² or can be obtained through relatively cheap ab initio calculations. Furthermore, the matrix elements can be evaluated in a fixed time regardless of the number of atoms in the system, and multiple temperatures and carrier concentrations can be calculated simultaneously with only a modest increase in the computational time. Full-timing information for the calculation of all first-principles inputs required to compute the transport properties of the materials discussed in this work and the scaling performance of each code routine is given in the Supplementary Methods.

Analysis of scattering rates and electron mobility

In Fig. 1, we compare mode-dependent scattering rates for n-Si and n-GaAs calculated by our method against fully first principles calculations (DFPT + Wannier) at 300 K obtained using the EPW and PERTURBO softwares^32,35. The scattering of electrons in Si is dominated by acoustic phonons whereas polar optical phonon scattering dominates in GaAs, as revealed by the mobility analysis in Supplementary Fig. 9. Excellent agreement is seen for both systems, with the onset of polar optical emission scattering in GaAs well described by our calculations. The high degree of agreement for Si is somewhat surprising considering our calculations do not include optical phonons that are known to contribute to scattering at room temperature. However, as we do not consider symmetry selection rules in our calculations, the acoustic phonon scattering rate will be slightly overestimated and is expected to partially pick up some of the missing optical deformation scatterings. Whereas in Si, group theoretical rules dictate that intervalley g-phonon transitions are only possible through longitudinal optical phonons, Rode⁶³ and others⁶⁴ have demonstrated the temperature-dependent transport behavior can reproduced considering only acoustic phonons if selection processes are ignored. Additional comparisons against DFPT + Wannier scattering rates for 3C-SiC and p-SnSe are provided in Supplementary Fig. 13. In both cases, the shape and magnitude of the scattering rates is well reproduced, particularly at low energies, despite the simpler approach that does not involve an expensive DFPT calculation to obtain the matrix elements.

Fig. 1

Calculated scattering rates.

Comparison of the calculated (pink) scattering rates, $$ {\tau }_{n\mathbf{k}}^{-1}$$, against those obtained using density functional perturbation theory combined with Wannier interpolation (DFPT + Wannier, light teal) for (a) n-GaAs⁴⁴ and (b) n-Si³² at 300 K.

In Fig. 2, we compare the time taken to compute the transport properties of NbFeSb and Ba₂BiAu (the full-timing breakdown is tabulated in Supplementary Tables 2 and 3). Taking into account the time required to compute all first-principles inputs and the electron mobility at a single temperature and carrier concentration, our method offers over a 2 order of magnitude speed up compared to DFPT + Wannier (an average of 29 core hours versus 8350 core hours). Considering only the time needed to obtain the scattering rates and transport properties (i.e., presuming all inputs have already been tabulated), our approach offers a 4 order of magnitude speed up (Fig. 2a). This can be exploited when performing calculations at multiple temperatures and carrier concentrations. For example, calculating the mobility of Ba₂BiAu for 10 temperatures requires ~32,000 core hours using DFPT+Wannier compared to <35 core hours with our approach (95% of which is required to calculate the first principles inputs). Furthermore, we expect the relative cost advantage of our method to increase with system size as unlike in DFPT + Wannier the computational expense of the matrix elements does not depend on the number of atoms. This reduction in computational time, combined with similar accuracy to DFPT + Wannier [within 10%, see Fig. 2b], makes our approach amenable to the large-scale calculation of electronic transport properties.

Fig. 2

Comparison of speed against accuracy for transport calculations.

Existing methods for calculating electron transport properties are either computationally efficient but inaccurate (constant relaxation time, CRT, orange) or accurate but highly computationally demanding (density functional perturbation theory combined with Wannier interpolation, DFPT + Wannier, teal). The approach outlined in this work (pink) demonstrates accuracy comparable to state-of-the-art methods at ~1/500th of the computational cost. a The time required to obtain electron mobility for each method is broken down by the time spent computing first-principles inputs and performing the scattering and transport calculations. b The mean absolute percentage error in the calculated mobility at 300 K is compared to the total computational time (including the time to obtain all first-principles inputs). Results are averaged for NbFeSb (p-type, n = 2 × 10²⁰ cm⁻³, DFPT + Wannier^39,90) and Ba₂BiAu (n-type, n = 1 × 10¹⁴ cm⁻³, DFPT + Wannier⁹¹). In (b), the mobility error is referenced with respect to state-of-the-art DFPT + Wannier calculations as high-quality experimental data were not available. The full timing breakdown for each material is provided in Supplementary Tables 2 and 3. Constant relaxation time calculations were performed with τ = 10 fs.

Figure 3a plots the calculated mobility of GaN against experimental measurements, indicating very close agreement from 150 to 500 K. As each scattering mechanism is treated with a separate matrix element, this allows the impact of individual scattering processes to be assessed. At low temperatures, the mobility of GaN is limited by impurity scattering, with polar optical phonon scattering dominating above 300 K, as illustrated by the dashed lines in Fig. 3a. The total mobility taking into account all scattering mechanisms reproduces the experimental mobility with very high agreement. Further insight into the competing nature of the scattering mechanisms is provided by the energy dependence of the electron lifetimes and the resulting spectral conductivity, Σ(ε) = v(ε)²τ(ε)N(ε) where N is the density of states and v is the group velocity, computed at 300 K and an electron concentration of 5.5 × 10¹⁶ cm⁻³ (Fig. 3b, c). Impurity scattering dominates at the conduction band edge but diminishes quickly as energy increases. At energies above ω_po of the band minimum (above the phonon emission threshold), polar-optical interactions are two orders of magnitude stronger than any other competing mechanism and act as the primary limiting factor for electron mobility, in agreement with the experimental findings of ref. ⁶⁵ and DFPT + Wannier calculations⁴². In contrast, the mobility calculated using a constant relaxation time of τ = 10 fs—a value on the higher end of that typically employed in screening studies^25,26,28—underestimates the mobility by a factor of 2–10 depending on the temperature, as shown in Fig. 3a. More fundamentally, the CRTA does not reproduce the correct shape of temperature dependence as depicted in Fig. 3c. The ability of our method to reproduce the qualitative temperature dependence of transport properties, as well as make good approximations of quantitative behavior (often closely in-line with more detailed theoretical methods), thus represents a major advance for improving the accuracy of high-throughput methods.

Fig. 3

Temperature-dependent transport properties of GaN and SnSe.

a Comparison of the electron mobility, μ, of GaN against experiment (black triangles,⁹²). Mobility limited by ionized impurity (teal, ii), acoustic deformation potential (orange, ad), and polar optical phonon scattering (pink, po) is indicated in dashed lines. Total mobility taking into account all scattering mechanisms ($$ 1/{\tau }_{n\mathbf{k}}^{\mathrm{ii}}+1/{\tau }_{n\mathbf{k}}^{\mathrm{ad}}+1/{\tau }_{n\mathbf{k}}^{\mathrm{po}}$$) is given by the black solid line. Constant relaxation time (CRT) calculations with τ = 10 fs is given by dotted gray line. b Electron lifetimes, τ_nk and (c) spectral conductivity, Σ, arising from different scattering processes in GaN at 300 K. The valence band maximum is set to zero eV. In (b), the vertical dotted gray line indicates the energy of the effective polar phonon frequency, ω_po. d Comparison of the direction-dependent mobility of SnSe against experiments—a (orange), b (teal), c (pink) points from ref. ⁶⁶, b–c (black) points from ref. ⁹³.

A primary goal of the present approach is to extend well-established scattering matrix elements that were formulated for isotropic materials properties to be compatible with highly anisotropic materials. To that end, we have calculated the direction-dependent hole mobilities of Pnma structured SnSe at a carrier concentration of 3 × 10¹⁷ cm⁻³, with the results compared to Hall measurements in Fig. 3d. Single-crystal SnSe has recently attracted significant attention as a thermoelectric material. Due to its layered structure, SnSe exhibits anisotropic transport properties, with the highest thermoelectric performance observed along the b axis⁶⁶. Our calculations reproduce the strong directional dependence in transport measurements, in which the mobility parallel to the layers (along b and c) is almost an order of magnitude larger than that perpendicular to the layers (along a). Our mobility results agree remarkably well with the considerably more computationally expensive electron–phonon calculations performed using DFPT+Wannier and G₀W₀ band structures⁴¹ (Supplementary Fig. 8). We note that additional anisotropy in the mobility between the b and c directions has been observed in high-temperature experimental measurements⁶⁶. In both our calculations and DFPT + Wannier, however, the mobility along b and c are almost the same for temperatures above 300 K⁴¹. The discrepancy against the experiment is thought to derive from the use of a Hall factor r_H of unity when extracting the carrier concentrations needed to compute mobility⁴¹ (further analysis is provided in the Supplementary Discussion). As we demonstrate in Supplementary Fig. 7, access to band and k-dependent lifetimes can further be used to calculate electron linewidths that are qualitatively comparable to those measured through techniques such as angle-resolved photoemission spectroscopy (ARPES)⁶⁷.

Electron mobility and Seebeck coefficient across many systems

To demonstrate the generality of our approach, we investigate the transport properties of 21 semiconductors ranging from 2 to 48 atoms in their primitive unit cells. To highlight the compatibility of the method with high-throughput computations, all inputs (eigenvalues, wave functions, materials parameters) are obtained from density functional theory (DFT) using low-cost exchange-correlation functionals (see the “Methods” section). All such materials parameters are listed in Supplementary Table 4. Results are compared to transport measurements on high purity single-crystalline samples to minimize the effects of grain boundaries and crystallographic defects. Further details on the calculation methodology and selection of reference data are provided in the Supplementary Methods. The materials span multiple chemistries, doping polarities, and band structure types including anisotropic and multiband systems, and comprise: (i) conventional semiconductors, Si, GaN, GaP, GaAs, InP, ZnS, ZnSe, CdS, CdSe, and SiC; (ii) the thermoelectric candidates SnS, SnSe, PbTe, Bi₂Te₃, and BiCuOSe; (iv) photovoltaic absorbers PbS and CdTe; and (iii) transparent conductors, SnO₂, ZnO, and CuAlO₂. Our dataset also includes the relatively complex CH₃NH₃PbI₃ hybrid perovskite containing 48-atoms. In Fig. 4a we compare calculated mobility against experimental measurements for all 21 materials in our dataset. Calculations were performed using the experimentally determined carrier concentrations at a temperature of 300 K. Results regarding the temperature and carrier concentration dependence of mobility for all materials (calculated, experimental, and comparison with CRTA) is provided in Supplementary Figs. 8 and 9 and include the breakdown of mobility by scattering type. These plots represent a comprehensive test of our approach, across many materials, not only at the single condition plotted in Fig. 4a but when conditions are varied.

Fig. 4

Calculated mobility and Seebeck coefficient against experiment.

a Comparison of carrier mobilities, μ, at 300 K between calculations and experiments, with points colored by the conductivity effective mass $$ m_c^{*}$$. b Comparison of Seebeck coefficients, α, at 300 K between calculations and experiments, with points colored by the majority carrier concentration n. For Si and CdS we compare directly to the diffusive component of Seebeck coefficient only. CH₃NH₃PbI₃ has been abbreviated as MAPbI₃. In (a) and (b), orange crosses indicate results computed using a constant relaxation time of 10 fs. Detailed temperature and carrier concentration results for each material are provided in Supplementary Figs. 8 and 11.

The calculated mobilities agree closely with experiment across all materials, covering several orders of magnitude from ZnO (180 cm²/Vs) to n-type GaAs ($$ {\mu }_{\exp }=$$2.1 × 10⁴ cm²/Vs). Notably, the calculated mobility (Spearman rank coefficient against experiment r_s = 0.93) improves significantly on results obtained using a constant relaxation time of τ = 10 fs (r_s = 0.52). As we demonstrate in Supplementary Fig. 14, our method also improves the dependence of mobility on temperature, obtaining a mean squared error (MSE) of 0.20 consistent with DFPT + Wannier (MSE = 0.19) and dramatically more accurate than a constant relaxation time (MSE = 3.2). We note that more reasonable results can be obtained in the constant relaxation time approximation by fitting the lifetime to experimental measurements. However, in most cases, experimental measurements will not be available and thus this approach cannot easily be applied to new or uncharacterized systems. One alternative is to linearly scale the constant lifetime by the temperature according to $$ {\tau }_T=\tau \frac{300}{T}$$ ⁶⁸. As demonstrated in Supplementary Figure 9, this results in a considerable improvement in the temperature dependence of mobility, reducing the mean squared error from 3.2 to 0.9. By definition, however, this does not impact the mobility at 300 K and thus does not change the constant lifetime results presented in Fig. 4a.

Our approach demonstratesgreater deviation from experiment in materials with smaller mobilities such as p-CuAlO₂ (3 cm²/Vs in a-b plane), where a local hopping mechanism is proposed to compete with band transport⁶⁹, and p-CdTe in which spin-orbit coupling (SOC) is known to dramatically impact the scattering rates at the valence band edge⁷⁰ but was not included in our calculations. Additional deviation is observed for n-ZnS, where the calculated mobility is almost a factor of 4 larger than Hall measurements. We find this overestimation is largely due to the underestimation of the conduction band effective mass arising from the use of the PBE exchange-correlation functional ($$ m_c^{*,\mathrm{PBE}}=$$ 0.16 m_e) when compared to experiment ($$ m_c^{*,\exp }=$$ 0.22 m_e)⁷¹. As we detail in Supplementary Fig. 10, calculations performed using the hybrid HSE06 functional result in a larger effective mass ($$ m_c^{*,\mathrm{HSE}}=$$ 0.20 m_e) and improved agreement with the experimental mobility. Conversely, the mobility of n-type PbTe and p-type SnS are underestimated relative to the experiment. In both compounds, this originates from an overestimation of the effective mass. For example, in PbTe, HSE06 with the inclusion of spin-orbit coupling effects is required to reduce the effective mass from $$ m_c^{*,\mathrm{PBE}}=$$ 0.11 m_e to $$ m_c^{*,\mathrm{HSE+SOC}}=$$ 0.03 m_e, in line with the experimental effective mass of 0.02 m_e. In Supplementary Fig. 10, we demonstrate that more accurate treatment of the electronic structure significantly improves the agreement against the experiment across the full temperature range. Lastly, the deviation seen for CH₃NH₃PbI₃ is likely due to the use of polycrystalline thin films in experimental measurements. As highlighted by Supplementary Fig. S8, our temperature-dependent results are in excellent agreement (within 5% at all temperatures) against fully first principles calculations performed using EPW³⁷. The ability of our approach to accurately describe the electron–phonon coupling of a highly complex structure with 144 phonon-modes while remaining computationally efficient highlights its potential in high-throughput screening of transport behavior.

Accurate calculation of Seebeck coefficients is of primary interest in the prediction and analysis of thermoelectric materials. In Fig. 4b we compare calculated Seebeck coefficients against those obtained experimentally at 300 K. A comparison of the temperature dependence of the Seebeck coefficient for all materials is provided in Supplementary Fig. 11. We see reasonable agreement against experiment across the full range of materials, for both p- and n-type samples, corresponding to positive and negative Seebeck coefficients, respectively. In Supplementary Fig. 12, we demonstrate that use of the HSE06 hybrid functional can further improve the agreement against the experiment for n-type GaAs. We note that for Si and CdS we compare directly to the diffusive component of Seebeck coefficient only, ignoring the effects of phonon drag which contribute substantially even at room temperature^72–75. The Seebeck coefficient displays a weaker dependence on electron lifetimes than mobility and conductivity and so is often treated within the CRTA (in which case the specific relaxation time cancels in the equations).²⁷. Figure 4b indicates that this approximation is often justified due to the relatively small disagreements between constant relaxation time and mode-dependent relaxation time results, in-line with previous comparisons of CRTA against experimental data⁷⁶.