Efficient first-principles electronic transport approach to complex band structure materials: the case of n-type Mg₃Sb₂

Abstract

We present an efficient method for accurately computing electronic scattering rates and transport properties in materials with complex band structures. Using ab initio simulations, we calculate a limited number of electron–phonon matrix elements, and extract scattering rates for acoustic and optical processes based on deformation potential theory. Polar optical phonon scattering rates are determined using the Fröhlich model, and ionized impurity scattering rates are derived from the Brooks-Herring theory. Subsequently, electronic transport coefficients are computed within the Boltzmann transport theory. We exemplify our approach with n-type Mg₃Sb₂, a promising thermoelectric material with a challenging large unit cell and low symmetry. Notably, our method attains competitive accuracy, requiring less than 10% of the computational cost compared to state-of-the-art ab initio methods, dropping to 1% for simpler materials. Additionally, our approach provides explicit information on individual scattering processes, offering an alternative that combines efficiency, robustness, and flexibility beyond the commonly employed constant relaxation time approximation with the accuracy of fully first-principles calculations.

Introduction

The discovery of novel materials has evolved into the paradigm of big data-driven science^1,2,3, which can benefit from theory (physical and chemical laws), experiment (empirical trial and error method), and computer simulation such as density functional theory (DFT) calculations. With the continuous development of artificial intelligence, new approaches, such as machine learning^4,5, have broad applications, ranging from image recognition to predictive analytics in materials science. Machine learning, fueled by vast datasets, requires high-quality training data to avoid erroneous, missing, or redundant information. To obtain large quantities of high-quality data, it is necessary to balance the accuracy of computing fundamental physical properties with the efficiency of identifying candidate performance. This is particularly the case for the electronic transport properties of solid-state materials, where the central quantities are the carrier scattering rates associated with intrinsic (electron–electron, electron–phonon) and extrinsic (impurities, grain boundaries, alloy disorder) scattering mechanisms. Even though electron–phonon (e–ph) interaction dominates the device performance at elevated temperatures, its detailed description is still complicated. However, the key transport parameters are crucial for novel material deployment in a variety of technological applications, including solar cells^6,7,8, solid-state batteries^9,10,11, light-emitting diodes (LED)^12,13, photocatalysis^14,15,16, thermoelectrics^17,18,19, and many more.

In recent years, there has been an increase of interest in the computational prediction of the electronic transport behavior in the myriad complex band materials that have emerged^20,21,22,23. One of the earliest and most common approaches is to solve the Boltzmann transport equation (BTE) in the constant relaxation time (CRT) approximation. This is due to the absence of knowledge about the actual relaxation times, and it assumes energy, momentum, band, and temperature independence of the relaxation time. Traditionally, the most commonly used code is BoltzTraP^24,25, based on a smoothed Fourier interpolation of the DFT electronic band structures. This code, for example, has computationally driven the field of thermoelectric (TE) materials for years^26,27,28,29. Its computational cost can be further reduced using a maximally localized Wannier function basis set as in the BoltzWann code³⁰. The CRT has been successfully applied to high-throughput searches for computing transport behavior³¹. However, the relaxation time is generally not a constant. Comparisons of electronic transport properties indicate significant quantitative and qualitative differences between those using energy/momentum/band-dependent relaxation times and the CRT³², which implies that optimizing and screening materials need to be accomplished beyond the CRT level of approximation.

DFT and density functional perturbation theory (DFPT) have enabled calculations of e–ph interactions from first principles. This procedure can be accelerated within the EPW code^33,34. By using Fourier interpolation techniques to capture the short-range contributions of the e–ph matrix elements and by computing the dipole Fröhlich term to capture their long-range e–ph contributions³⁵, this DFPT + Wannier method can accurately predict the relaxation time for all the acoustic and optical phonon modes, including the intra-valley and inter-valley scattering, and in many cases shows a remarkable agreement with experimental measurements^36,37,38. Most recently, a way to include the quadrupole e–ph interaction beyond the Fröhlich interaction has been developed as well^39,40. The quadrupole e–ph interaction is found to be sizable in certain non-polar materials. Taking it into account corrects the dipole term in polar materials and better describes the long-range interactions for obtaining more accurate transport results. However, to capture the delicate features of the scattering near the band edges, ultra-fine grids for the sampling of electron and phonon wavevectors in the Brillouin zone are needed, which results in a vast number of computed matrix elements. For example, for cubic Si, which contains two atoms in the unit cell, hundreds of thousands of k-points are needed to obtain convergence in the intrinsic mobility³⁷, and typically the number of matrix elements is of the order of ${N}_{{{{\bf{k}}}}}^{2}$ (N_k is the number of k-points). Thus, currently, the computation of e–ph relaxation times and charge transport within the BTE using the DFPT + Wannier method is highly resource-intensive for materials with larger unit cells (containing more atoms and basis functions) and lower symmetry (featuring larger non-equivalent k-space regions)^41,42. To reduce computation cost while maintaining good accuracy based on the DFPT + Wannier method, several methods are available that involve averaging or integrating matrix elements over the entire Brillouin zone, such as EPA²⁵ and EPIC STAR⁴³. Additionally, the dual interpolation method has been found to be very effective, being several orders of magnitude faster compared to the DFPT + Wannier method⁴⁴.

Instead of using the DFPT + Wannier method, one can also address this challenge by using deformation potential theory to compute the relaxation times. The initial concept of the deformation potential for the e–ph interaction in the long wavelength limit for acoustic phonons was introduced by Bardeen and Shockley in 1950⁴⁵, and it has been a driving force in simulations of semiconductor devices ever since. This method has also recently been included within first-principles calculations by computing the band shift of the valence band maximum (VBM) or conduction band minimum (CBM) due to volume change and using that to calculate transport properties, for example, as performed in the AMSET code⁴⁶. In this case, a single acoustic deformation potential is extracted from the shift of electron energies under external strain. No matrix elements are computed and these common materials parameters used can be obtained through relatively cheap ab initio calculations or found in materials databases⁴⁷. Thus, the method is computationally inexpensive, enabling studies of a large number of materials with much better accuracy compared to the CRT.

In this work, we introduce an efficient formalism to combine the DFPT + Wannier method with the deformation potential theory, offering an alternative direction to calculate transport properties which provides efficiency, robustness, and flexibility. We calculate the acoustic, optical, and inter-valley deformation potentials from e–ph matrix elements using first-principles calculations. Scattering rates are then calculated using deformation potential theory and Fermi’s Golden Rule by considering their full energy, momentum, and band dependence. We complete the calculation of the overall scattering rates by computing polar optical-phonon and ionized impurity scattering rates. Using our own-developed open-source Boltzmann transport code ElecTra⁴⁸, we validate our approach by performing an in-depth investigation of scattering rates, mobility, and Seebeck coefficient for the promising TE material n-type Mg₃Sb₂ and compare to experiments, the CRT, and DFPT + Wannier calculations. The material we use as an example is chosen for its band structure complexity, unit cell size, and degree of symmetry, which make fully ab initio calculations computationally prohibitive. We achieve excellent agreement with the DFPT + Wannier method while utilizing no more than 10% of its computational cost. Applying the same approach to Si, a simpler material, we demonstrate once again that we can attain ab initio accuracy, this time at less than 1% of the corresponding ab initio computational cost. Our method belongs to the category of methods that compute and process matrix elements. However, it distinguishes itself through advancements in accuracy and flexibility. Firstly, we ensure accuracy by selectively computing crucial matrix elements at specific energies and wavevectors, focusing on regions responsible for electronic transitions. This allows us to afford dense grids around these significant areas. Secondly, our approach provides explicit information on individual scattering processes (acoustic, optical, intra- and inter-valley), offering valuable insights and capabilities that are particularly advantageous for designing materials with optimal multi-valley electronic structures⁴⁹.

Results

Carrier scattering approach

Figure 1a displays the trigonal structure of Mg₃Sb₂ (space group $P\bar{3}m1$), which is the one of the CaAl₂Si₂-type Zintl compounds with complicated chemical bonds⁵⁰. Even though almost all Zintl materials are found to be p-type, Mg₃Sb₂ has recently been synthesized as n-type using excess Mg atoms⁵¹. Figure 1b shows the band structure of Mg₃Sb₂. Note that the band gap is undervalued within the PBE functional used but can be increased by using hybrid density functional HSE06 while the shape of the band structure is basically unaltered⁵². Since the experimental band gap is relatively large at 0.8 eV⁵³, bipolar effects are not considered. Consequently, the band gap is not factored into our calculations for unipolar transport in the conduction band. The CBM is located along the L^*–M^* line, slightly off the L–M high symmetry line (see Supplementary Fig. 1). In addition, the K-point valley is only 0.18 eV above the CBM. Thus, the vicinity of its conduction band has a high-degeneracy Fermi surface, with six-fold carrier pockets from the CBM and six more pockets from the K point. Only one-third of each K-valley is located within the first Brillouin zone, giving an overall K-valley degeneracy of two. This conduction band complexity makes it possible to achieve outstanding n-type TE performance^54,55, as a large valley degeneracy is one of the factors for high electronic transport^32,56.

Fig. 1: Carrier scattering.

a Atomic structure for Mg₃Sb₂. b Electronic band structure for Mg₃Sb₂, where the orange zone shows valleys contributing to electron transport. c A schematic depicting intra-valley scattering, e.g., at the conduction band minimum (CBM). d Illustration for inter-valley scattering from the CBM to other equivalent CBM valleys and to the K valley. The Fermi surfaces depicted for the CBM and K valleys are positioned on different energies, i.e., separately at energies 0.05 eV higher than their corresponding minima for better illustration purposes. See Supplementary Fig. 3 for the original figures. e A schematic showing the long-range and short-range matrix element contributions to electron–phonon scattering. The blue and green zones show the small and large ∣q∣ related to the intra-valley and inter-valley deformation potential scattering, respectively. The inset shows the phonon spectrum for Mg₃Sb₂, where the purple and yellow zones indicate the optical and acoustic phonon modes considered in extracting the ODP and ADP, respectively.

The complexity of the conduction bands of Mg₃Sb₂ is explicitly chosen to demonstrate our method’s flexibility in considering both intra-valley and inter-valley scattering. Figure 1c shows the various intra-valley scattering mechanisms that take part, i.e., acoustic deformation potential (ADP) scattering, optical deformation potential (ODP) scattering, polar optical phonon (POP) scattering, and ionized impurity scattering (IIS). Figure 1d shows the various inter-valley scattering (IVS) processes (see also Supplementary Fig. 2), where the 2nd Brillouin zone is also considered and allows the selection of the nearest scattering distance⁵⁷. Both the CBM and K valley are considered as the final states for IVS.

For non-polar materials, we can derive all the deformation potentials from the e–ph matrix elements M_mnν(k, q)⁵⁷. The M_mnν(k, q) is defined as^58,59

$$\begin{array}{r}{M}_{mn\nu }({{{\bf{k}}}},{{{\bf{q}}}})=\left\langle {u}_{m{{{\bf{k}}}}+{{{\bf{q}}}}}\left\vert {\delta }_{{{{\bf{q}}}}\nu }V\right\vert {u}_{n{{{\bf{k}}}}}\right\rangle \end{array}$$

(1)

where u_mk+q and u_nk are the lattice-periodic part of the electronic wavefunctions for the final state on band m with wavevector k + q and initial state on band n with wavevector k, respectively. δ_qνV is the the derivative of the self-consistent potential associated with a phonon of wavevector q and branch index ν that can be computed by DFPT calculations. Bardeen and Shockley defined the ADP as the shift of the band edges due to local strain⁴⁵. Acoustic phonons under the long wavelength limit generate a corresponding dilation of the atomic spacing, which gives rise to a perturbing potential that tends to shift the electronic band energy. In this case, the ADP can be defined by ${D}_{{{{\rm{ADP}}}}}^{{{{\rm{b}}}}}=\delta H| {{{\bf{S}}}}{| }^{-1}$, where b denotes ‘band shift’, δH is the perturbing potential, and S is the strain. The ADP for non-polar materials can also be calculated from DFPT calculations by ${D}_{{{{\rm{ADP}}}}}^{{{{\rm{e}}}}-{{{\rm{ph}}}}}={M}_{mn\nu }({{{\bf{k}}}},{{{\bf{q}}}})| {{{\bf{q}}}}{| }^{-1}$. Khan and Allen proved that the two deformation potential quantities from the band shift and DFPT are the same at the semiconductor band extrema⁶⁰. On the other hand, Harrison extended the ADP theory by considering the non-polar optical phonon interaction with carriers⁶¹. In this case, the perturbing potential is proportional to the atomic displacement u. The ODP can be defined as D_ODP = δH∣u∣^{−1 62}, obtained from DFPT by ${D}_{{{{\rm{ODP}}}}}^{{{{\rm{e}}}}-{{{\rm{ph}}}}}={M}_{mn\nu }({{{\bf{k}}}},{{{\bf{q}}}})$.

For polar materials as the Mg₃Sb₂ at hand, the relevant e–ph matrix elements M_mnν(k, q) for longitudinal phonons consist of two parts. The first is the short-range part, which leads to the deformation potential scattering described above. The second is the long-range interaction, i.e., the Fröhlich interaction, which is related to the coupling between electrons and polar optical phonons⁶³. In general, at small ∣q∣ vectors, the Fröhlich interaction due to long-wavelength longitudinal optical phonons obeys ∣q∣⁻¹, while interactions from polar acoustic phonons are of piezoelectric nature and obey ∣q∣^−1/2, as illustrated in Fig. 1e. Note that piezoelectric scattering is important at very low temperatures in very pure semiconductors, so we will not discuss much piezoelectric scattering here since we consider highly doped Mg₃Sb₂ at room temperature, plus our acoustic mode matrix element calculations did not show any ∣q∣^−1/2 behavior for this material (Supplementary Fig. 4). The recently investigated quadrupole e–ph interaction could provide corrections beyond the dipole Fröhlich term^39,40, and it is present in both polar and non-polar materials. In comparison to the Fröhlich interaction, which follows an order of ∣q∣⁻¹ at small ∣q∣, the quadrupole matrix elements converge to a constant value of order ∣q∣^{0 64}. Consequently, they do not exert a dominant influence over the long-range part in this context. The Fröhlich interaction is defined as

$$\begin{array}{ll}{M}_{mn\nu }^{L}({{{\bf{k}}}},{{{\bf{q}}}})\,=\,i\frac{{{m}_{0}}^{1/2}{e}^{2}}{\Omega {\varepsilon }_{0}}\mathop{\sum}\limits_{k}{{M}_{k}}^{-1/2}\\\qquad\qquad\qquad\times\,\mathop{\sum}\limits_{{{{\bf{G}}}}\ne -{{{\bf{q}}}}}\frac{({{{\bf{q}}}}\,+\,{{{\bf{G}}}})\cdot {{{{\bf{Z}}}}}_{k}^{* }\cdot {{{{\bf{e}}}}}_{k\nu }({{{\bf{q}}}})}{({{{\bf{q}}}}\,+\,{{{\bf{G}}}})\cdot {k}_{\infty }\cdot ({{{\bf{q}}}}\,+\,{{{\bf{G}}}})}\\\qquad\qquad\qquad\times\,\langle {u}_{m{{{\bf{k}}}}+{{{\bf{q}}}}}| {e}^{i({{{\bf{q}}}}+{{{\bf{G}}}}){{{\bf{r}}}}}| {u}_{n{{{\bf{q}}}}}\rangle \end{array}$$

(2)

where m₀ is the sum of the masses of all the atoms in the unit cell, Ω is the unit cell volume, ε₀ is the vacuum permittivity, k_∞ is the high-frequency dielectric constant, M_k is the mass for atom k, G indicates the reciprocal lattice vector, and ${{{{\bf{Z}}}}}_{k}^{* }$ is the Born effective charge tensor.

By subtracting the long-range part from the total matrix elements³⁵, we can get the short-range matrix elements for the relevant phonons, as shown in Fig. 1e, which are then used to extract deformation potentials. The deformation potentials excluding the polar singularity are defined as ${D}_{{{{\rm{ADP}}}}}^{{{{\rm{e}}}}-{{{\rm{ph}}}}}=[{M}_{mn\nu }({{{\bf{k}}}},{{{\bf{q}}}})-{M}_{mn\nu }^{L}({{{\bf{k}}}},{{{\bf{q}}}})]| {{{\bf{q}}}}{| }^{-1}$ and ${D}_{{{{\rm{ODP}}}}}^{{{{\rm{e}}}}-{{{\rm{ph}}}}}={M}_{mn\nu }({{{\bf{k}}}},{{{\bf{q}}}})-{M}_{mn\nu }^{L}({{{\bf{k}}}},{{{\bf{q}}}})$.

In the case of IVS, which scatters a carrier from a given valley into another one with very large changes in momentum, only phonons with wave vectors near the zone boundary participate. These transitions can be acoustically or optically mediated and can happen within the first or even into the second equivalent Brillouin zone. Away from the Γ point, both the acoustic and optical phonon modes are flat and of high energy. Thus, for both, the matrix elements for IVS resemble ODP scattering, and the inter-valley deformation potentials can be defined as ${D}_{{{{\rm{IVS}}}}}^{{{{\rm{e}}}}-{{{\rm{ph}}}}}={M}_{mn\nu }({{{\bf{k}}}},{{{\bf{q}}}})-{M}_{mn\nu }^{L}({{{\bf{k}}}},{{{\bf{q}}}})$. Note that the long-range contributions are reduced significantly at these large phonon wavevectors required for IVS, as shown in Fig. 1e.

Matrix elements for Mg₃Sb₂

To obtain the deformation potentials and scattering rates for Mg₃Sb₂, we first identify all possible electronic transitions, the participating phonon modes, and their wavevectors. Then, we compute the total matrix elements for all those. Next, we calculate the long-range part of the e–ph matrix elements as in Eq. 2. For generality and computational consistency, we perform this operation for all modes (acoustic and optical, longitudinal and transverse)³⁵, where in the case of non-polar modes, the long-range coupling strength is suppressed. After subtracting the long-range part from the total matrix elements, what remains is the short-range part, which is used for deriving the deformation potentials. Figure 2 shows examples of the long-range, short-range, and inter-valley matrix elements, where the first two are relevant for intra-valley transitions. In Fig. 2a, the long-range matrix elements for the optical phonon modes along the Γ–M direction show that only two out of the four longitudinal optical phonon modes are found to follow the ∣M∣ ~ ∣q∣⁻¹ relation at small ∣q∣ vectors, thus being polar (see Supplementary Fig. 5 for these along other high-symmetry directions as well). The other two longitudinal modes (and all transverse modes) have long-range values of nearly zero. For the polar modes identified, we then compute their phonon energies. The inset of Fig. 2a shows the energy of those two polar optical phonon modes along the Γ–M direction. By calculating the average value of this polar optical phonon energy along different high-symmetry directions and then by averaging those values from all polar modes, we extract one dominant frequency for the overall polar optical phonon scattering process, which in this case is identified to be 26 meV.

Fig. 2: Electron–phonon coupling matrix elements for Mg₃Sb₂.

a Long-range part of the matrix elements along the Γ–M line for optical phonons. The dashed-doted colored lines show the ∣M∣ ~ ∣q∣⁻¹ dependence for polar optical phonons near the Γ point. Only two modes have polar behavior. The energy of these two modes is shown in the inset. b Short-range part of the matrix elements along the Γ–M line for the LA and (c) optical phonon modes, with the initial electron state located at the CBM. d Overall inter-valley deformation potentials for the four inter-valley scattering processes (green squares): three from the CBM to the other equivalent CBM and one from the CBM to the K valley. The gray dots show the individual contributions of every phonon mode to each process. The inset illustrates the four processes from the CBM to the final scattered states within the first Brillouin zone.

Figures 2b, c show examples of the short-range part of e–ph coupling matrix elements for acoustic and optical phonon modes. In the case of acoustic modes, the short-range longitudinal acoustic matrix elements increase almost linearly. Then, we consider all high-symmetry equivalent crystallographic orientations⁵⁷ (see Supplementary Figs. 6, 7, and Supplementary Table 1). Afterwards, the overall ADP is calculated from both longitudinal acoustic and transverse acoustic phonon modes by ${D}_{{{{\rm{ADP}}}}}={v}_{s}{({D}_{{{{\rm{ADP}}}},{{{\rm{LA}}}}}^{2}{v}_{l}^{-2}+{D}_{{{{\rm{ADP}}}},{{{\rm{TA}}}}}^{2}{v}_{t}^{-2})}^{1/2}$, where v_s, v_l, and v_t are average sound velocity, longitudinal sound velocity, and transverse sound velocity, respectively. In the case of optical modes, the short-range optical matrix elements fluctuate slightly around a constant value. We compute the overall ODP by merging the short-range part of all the optical modes ν near the Γ point as

$$\begin{array}{r}{D}_{{{{\rm{ODP}}}}}^{{{{\rm{e}}}}-{{{\rm{ph}}}}}=\sqrt{{\omega }_{{{{\rm{ODP}}}}}\mathop{\sum}\limits_{\nu }\frac{{\left[{M}_{mn\nu }({{{\bf{k}}}},{{{\bf{q}}}})\,-\,{M}_{mn\nu }^{L}({{{\bf{k}}}},{{{\bf{q}}}})\right]}^{2}}{{\omega }_{\nu }}}\end{array}$$

(3)

where ω_ODP is the average value of the optical phonon energy involved in ODP scattering. For electrons in Mg₃Sb₂, the overall ADP and ODP are calculated to be D_ADP = 3.37 eV and D_ODP = 1.36 eV Å⁻¹, respectively, with ω_ODP= 23.00 meV.

Regarding IVS, the Fermi surface for electrons in Mg₃Sb₂ consists of six full CBM valleys and six one-third K valleys within the first Brillouin zone. To capture the nearest transition path, we consider the final points in both the first and second Brillouin zones (Fig. 1d). Four types of inter-valley transitions can be identified. From those, the first three are transitions from the initial CBM valley to its first, second, and third nearest neighbor CBM valleys. The fourth is a transition from the CBM to the K point. Starting from the general consideration that both acoustic and optical phonons contribute to inter-valley scattering, we calculate the matrix elements versus the phonon energy involved for all phonon modes that participate (Fig. 2d), where the initial electronic state k is located at one CBM valley and the final electronic state k + q is located at another CBM or the K valley. For each of the four scattering processes, the overall inter-valley deformation potential D_IVS can be calculated by combining the short-range parts of all the relevant acoustic and optical phonons as

$$\begin{array}{r}{D}_{{{{\rm{IVS}}}}}^{{{{\rm{e}}}}-{{{\rm{ph}}}}}=\sqrt{{\omega }_{{{{\rm{IVS}}}}}\mathop{\sum}\limits_{\nu }\frac{{\left[{M}_{mn\nu }({{{\bf{k}}}},{{{\bf{q}}}})\,-\,{M}_{mn\nu }^{L}({{{\bf{k}}}},{{{\bf{q}}}})\right]}^{2}}{{\omega }_{\nu }}}\end{array}$$

(4)

where ω_IVS is the weighted average value for phonons in a given IVS process. Transitions involved in these four processes have different strengths, as seen from the calculated inter-valley deformation potentials. For scattering from the initial CBM valley to the other three equivalent CBM valleys, D_IVS = 0.83, 0.69, and 1.15 eV Å⁻¹, with ω_IVS = 23.5, 21.3, and 17.6 meV, respectively. For transition processes from the CBM to the K valley, D_IVS = 1.08 eV Å⁻¹ with ω_IVS = 18.9 meV.

Scattering rates and transport calculation

The deformation potentials and the dominant frequency for polar optical phonons will be used as inputs to the BTE simulator for transport. Other than those, additional input parameters needed (dielectric constants, mass densities, etc.) are listed in Supplementary Table 2, which can be obtained through relatively cheap ab initio calculations or are already tabulated in databases such as the Materials Project⁴⁷. The scattering rate expressions are described in the Methods below, and the transport calculations are carried out using our open-access ElecTra code⁴⁸. In Fig. 3a, we compare our calculated total energy-dependent e–ph scattering rates for n-type Mg₃Sb₂ at room temperature (300 K) against fully first-principles calculations (DFPT + Wannier) obtained using the EPW code. The shape and magnitude of the phonon-limited scattering rates are well reproduced, particularly at low energies. For carriers with higher energy (>0.1 eV), which generally do not affect transport significantly, our scattering rates are somewhat lower compared to the DFPT + Wannier calculations, as we use only a few relevant matrix elements around the valley extrema to derive the deformation potentials, and don’t extend to higher energies. The overall e–ph scattering rates are found to be much higher than the ionized impurity scattering (IIS) rates (see Fig. 3a), even at large impurity concentrations. This agrees well with previous experiments⁶⁵, where the mobility of n-type Mg₃Sb₂ near 300 K is found to be dominated by e–ph scattering instead of IIS. This is different from that observed in common semiconductors such as Si⁵⁷.

Fig. 3: Calculated scattering rates and transport properties for Mg₃Sb₂.

a Comparison of the calculated scattering rates against those obtained using density functional perturbation theory combined with Wannier interpolation (DFPT + Wannier) at 300 K from ref. ⁶⁶. Calculated scattering rates for ionized impurity scattering at different impurity concentrations are also shown. b Scattering rates for acoustic deformation potential (ADP), optical deformation potential (ODP), inter-valley deformation potential (IVS), and polar optical phonon (POP) scattering. The vertical dotted line indicates the energy of the effective polar phonon frequency. c Comparison of the electron mobility limited by only electron–phonon scattering (e–ph) and by both electron–phonon and ionized impurity scattering (e–ph + IIS) against experiments from refs. ^{54,65,67,68,69}. d Comparison of the electronic conductivity against those obtained by DFPT + Wannier⁶⁶ and experiments^{54,65,67,68,69}. e Comparison of the Seebeck coefficient at 300 K using the constant relaxation time (CRT) approximation, e–ph, and e–ph + IIS conditions, and against those obtained by DFPT + Wannier⁶⁶ and experiments^{51,54,65,67,68,69,70}. f Comparison of the thermoelectric power factor at 300 K between calculations and experiments from refs. ^{54,65,67,68,69}.

Furthermore, our method considers different e–ph processes separately, so that we can evaluate their strength independently (Fig. 3b). It is clear that polar optical phonon (POP) scattering dominates the e–ph scattering at 300 K. The ADP scattering is the second strongest mechanism. It is somewhat higher but comparable to the IVS processes. The ODP scattering rates are lower. Although POP dominates the phonon scattering processes in this material, in general, it is essential to consider all processes when assessing electronic transport, which provides an in-depth understanding of the internal scattering processes. The POP-dominant scattering in n-Mg₃Sb₂ leads to the relatively weak influence of IIS on the electron mobility, as shown in Fig. 3c. This supports previous experimental data that the IIS is not dominant in Mg₃Sb₂, even at the high charge carrier concentrations near those needed for TE applications⁶⁵. Since we are able to separate all processes, it is useful to note the strength of intra- versus inter-valley transitions, which is an essential piece of information in studies of band alignment in TE materials. As an indication for the relative strength of these processes, by averaging the scattering rates in the energy range of 0–0.1 eV, without POP we find that the ratio of the intra- to inter-valley transitions is r_intra/inter = 3.0, whereas if POP is included in the intra-valley transitions, the ratio jumps to r_intra/inter = 24.8, clearly indicating that intra-valley processes dominate.

Since Mg₃Sb₂ is a very promising TE material for room temperature applications⁵⁴, we now compute the TE transport properties, i.e., electrical conductivity, Seebeck coefficient, and power factor (PF), as shown in Figs. 3d–f, respectively. We compare our calculations with experimental data and other calculations^{51,54,65,66,67,68,69,70}. Figure 3d shows the phonon-limited conductivity, which reproduces very well the results from DFPT + Wannier⁶⁶. Note that in both our calculations and DFPT + Wannier, the conductivity is larger compared to experiments^{54,65,67,68,69}, even when the IIS is included in our calculations. This difference is expected, as these materials typically have a large number of defects (e.g., grain boundaries, neutral impurities, and crystallographic defects). The influence of such defects is known to be less for higher carrier densities⁷¹. Also, note that even in single crystals, there are uncertainties arising from the difference between the Hall mobility and drift mobility, leading to a corresponding error in the indicated carrier concentrations, and additional differences in computational versus measured data⁷².

On the other hand, the Seebeck coefficient is not sensitive to the details of the relaxation time. All simulation methods⁶⁶, even the CRT approximation, are adequate to provide comparable results to experiments^{51,54,65,67,68,69,70}, as shown in Fig. 3e. Thus, the calculated PF (Fig. 3f) is higher compared to experiments^{54,65,67,68,69}, a reflection of the higher conductivity in both our calculations and DFPT + Wannier methods. Importantly, for thermoelectrics, which operate at relatively high densities, we also find that the IIS has a negligible effect on the PF for this material. In contrast, for materials such as Si, the PF can be decreased by 64% (see Supplementary Fig. 8) at high densities. Therefore, using materials like n-type Mg₃Sb₂, which appear to be immune to IIS, could represent a promising direction for thermoelectrics. Indeed, this is one of the reasons that make Mg₃Sb₂ a competitive alternative to the state-of-the-art n-type TE alloys near room temperature.

It is interesting to note that the transport properties of n-type Mg₃Sb₂ are nearly isotropic (see Supplementary Fig. 9), despite the fact that the crystal structure is anisotropic and every single carrier pocket of the CBM is anisotropic along the three axial directions in k-space. This is because the average effective mass tensor is nearly isotropic after considering the six equivalent carrier pockets. The conductivity effective masses along the different directions are calculated to be m_cond,x = 0.1646 m_e, m_cond,y = 0.1845 m_e, m_cond,z = 0.1656 m_e (m_e is electron mass). The nearly isotropic effective mass has also been reported in previous works using the parabolic band model, where the effective masses are estimated by fitting the near-edge band dispersions: m_in-plane = 0.22–0.24 m_e, m_out-of-plane = 0.16–0.17 m_e^52,66,73; m_in-plane = 0.32 m_e, m_out-of-plane = 0.28 m_e^54,55,74 It should be noted that the effective masses obtained using the parabolic band model ignore the effect of non-parabolicity, warping, as well as the presence of other bands and valleys at slightly higher energies and in different orientations. The mass values we report above are calculated by averaging over all contributing electronic bands weighted by their occupancy using our home-developed Effective Mass Finder (EMAF) code, which accounts for all band non-uniformities^32,75,76. In the same way, we also compute the overall density of states effective mass for n-type Mg₃Sb₂ to be m_DOS = 0.6788 m_e.

Computational performance

We now compare the performance of our method with that of the DFPT + Wannier for mobility calcuations in Si and Mg₃Sb₂, and demonstrate the significant computational savings our method offers. The total computation can be divided into two parts: the first is the computation of e–ph matrix elements using DFT and DFPT calculations, which is common in both methods; the second is the actual transport computation, where one needs a relatively dense mesh after interpolation to achieve convergent results in the DFPT + Wannier method. For example, for Si, around 500,000 k-points are found to be necessary for mobility convergence³⁷. For GaAs, using a 320 × 320 × 320 grid increases mobility by over 20% compared to the 150 × 150 × 150 mesh³⁶, indicating that a denser mesh for k-points could be necessary. This demonstrates the difficulty in achieving convergent results and, in general, the actual feasibility of fully ab initio electronic transport studies in complex materials because such meshes present enormous computational challenges. In the case of Mg₃Sb₂, a mesh of 75 × 75 × 50 was used for k-points in DFPT + Wannier (Fig. 3)⁶⁶. No convergence information is available due to the expensive cost of complex materials. However, the excellent match to our results in Fig. 3 could signal that convergence is reached.

Figure 4 shows the time required using the EPW code for transport simulations with respect to the number of electron wavevectors in the Brillouin zone. To be consistent with the definition of numbers of k-points N_k in the convergence test as in ref. ³⁷, we use the product of the meshes to be interpolated k_xk_yk_z as the N_k. Uniform grids are used with the k-points being twice as dense as the q-points in each direction. Reliable converged transport results using this DFPT + Wannier method require computational times at the right side of Fig. 4a, where nearly 500,000 k-points are used. This is a computationally expensive process because the mesh of q-points typically increases with that of k-points. Thus, the number of matrix elements to be computed is proportional to their product and increases with the mesh size. Although these result from interpolation, their number is very large, and thus the computation is quite heavy.

Fig. 4: Comparison of computation time and accuracy in transport calculations.

a The CPU time needed for the mobility calculation of Mg₃Sb₂ and Si by using density functional perturbation theory combined with Wannier interpolation (DFPT + Wannier) within the EPW code with different k-points. The green-shaded region indicates the mesh size required for convergent mobility calculations. The inset shows the initial DFPT time requirements. Our method requires the time in the inset plus the left-most values in the main figure. b Mobility computed by the ElecTra code considering only ADP and POP scattering, then all relevant scattering mechanisms, i.e., ADP, ODP, IVS, and POP, and compared against those obtained using DFPT + Wannier^37,66,77,81.

On the other hand, in our deformation potential method, the transport simulation is limited by the extraction of the deformation potentials, for which only a few dozen matrix elements are needed. In this case, the computational time required for transport calculations coincides with that on the left side of Fig. 4a. For Si, our transport time is only 10 core hours, significantly outperforming the DFPT + Wannier method, which is expected to require 50,000 core hours. This remarkable efficiency arises from a substantial reduction in the number of the matrix elements involved. For Mg₃Sb₂, the computational cost is considerably higher than for Si. The transport simulation is projected to demand 100,000 core hours in the DFPT + Wannier method, while our method completes the task in just 2,000 core hours.

For the first part of the overall computation, we show the computational time required for the dominant DFPT part in the inset of Fig. 4a. For Si, the initial 12 × 12 × 12 k-points and 6 × 6 × 6 q-points are used, while for Mg₃Sb₂, the initial 12 × 12 × 8 k-points and 6 × 6 × 4 q-points are used. These are the initial meshes used in the transport calculation for computation time in Fig. 4a. We also show the time with the initial 5 × 5 × 3 q-points for the more computationally demanding Mg₃Sb₂, a size used typically in literature⁶⁶, whose data we have shown in Fig. 3.

The computational time required at this first step is significant, but it is still at least an order of magnitude lower compared to the converged transport simulations. Since the transport part in our method is computationally much less expensive, the DFPT + Wannier part dominates the overall computation for both Si and Mg₃Sb₂, irrespective of the q-grid used. Thus, overall, for simple materials (small unit cells, high symmetry) such as Si, our deformation potential method requires less than 1,000 core hours while DFPT + Wannier is limited by the transport part, which makes our method 100 times more efficient. For more complex materials such as Mg₃Sb₂, both methods require 10,000 core hours for the DFPT calculation. Our method is overall 10 times faster than DFPT + Wannier. This difference can be more significant if a coarser q-grid is used, e.g., the 5 × 5 × 3. We also note that in cases involving larger unit cells, the absolute cost can be significantly higher compared to simpler materials, particularly in terms of the extra effort required to establish confidence in transport convergence. Thus, the high computational costs of the DFPT + Wannier method would make it daunting for its widespread use and quantitative confidence in its results.

Concerning accuracy, in Fig. 4b, we compare the mobility calculated at 300 K using our deformation potential method in the ElecTra code versus the fully ab initio DFPT + Wannier method in the EPW code. We compare calculations for three materials (full study for Si presented in ref. ⁵⁷), i.e., n-type Si, p-type Si, and n-type Mg₃Sb₂. Note that different values are reported from DFPT + Wannier works in various publications, as shown by the error bars in Fig. 4b, which may arise from details in the first-principles calculation of band structures, different grids for Brillouin-zone sampling, and various pseudopotentials for exchange and correlation³⁷. The mobilities calculated using our method for these three examples agree very well with DFPT + Wannier^{37,66,77,78,79,80,81}, which comes from the consideration of all scattering physics elaborately, by including contributions from all phonon processes (ADP, ODP, IVS, and POP) and distinguishing between intra- versus inter-valley processes. Notably, despite extracting only a few matrix elements and processing them into deformation potentials, we strategically focus on the essential transport regions, resulting in a locally dense matrix element grid. This is potentially advantageous compared to selecting matrix elements on a sparse grid across the entire Brillouin zone. In Fig. 4b, we also present results excluding ODP and IVS, considering only ADP for Si and ADP + POP for Mg₃Sb₂. In the case of Si, the mobility is highly overestimated, while for the polar Mg₃Sb₂, a better match is found, which signals the variation in behavior observed when ignoring ODP and IVS.

Discussion

With regard to the calculation of transport properties, there are many approaches available with varying degrees of complexity, reflecting the accuracy of the results versus the computational cost. Codes such as EPW and PERTURBO use DFPT + Wannier^34,82, accounting for all the scattering mechanisms but at a substantial computational cost, especially for materials with large unit cells and complex band structures. To reduce computational costs, the EPA method follows a simplified approach by averaging the matrix elements g across the Brillouin zone into an effective lifetime⁸³. However, such treatment may not be accurate enough in cases where g² varies rapidly. To overcome this problem, EPIC STAR uses the generalized Eliashberg function from g²ω^{−1 43}, which varies slower than g² at a small ∣q∣, although g²ω⁻¹ for acoustic phonons could also vary linearly for long wavelength and introduce some uncertainties. The dual interpolation method recently reported allows for a couple of orders of magnitude computational cost reduction compared to state-of-the-art EPW calculations⁴⁴. Such methods are appropriate when focusing on overall mobility and transport properties. However, their accuracy depends on the grid used, and the detailed information for the individual processes can also be further refined.

On the other hand, codes such as AMSET focus more on computational efficiency and scalability without the calculation of e–ph matrix elements. AMSET extracts the ADP from band shifts after applying stress along specific directions and uses that to form the scattering rate expressions⁴⁶, which is a very effective method to obtain transport properties at a very low computational cost. However, it pays the price of larger potential inaccuracies in capturing the details of the underlying physical mechanisms correctly, e.g., ODP and IVS, which may be critical for materials with strong ODP scattering or multiple valleys. The latter (multi-valley character) is one of the most prominent directions for obtaining high TE performance³².

Our method belongs to the category of methods that compute matrix elements, but does not require computations across the entire Brillouin one. For the extraction of deformation potentials, we use a few matrix elements near the specific q vectors for intra/inter-valley scattering, as shown in Fig. 5a, and these can be selected on a dense mesh. More importantly, by treating each q and phonon branch separately, we have access to the information about the underlying physical processes, such as intra-valley transitions versus inter-valley transitions. For example, Fig. 5b shows the scattering rates from each of these transitions, a capability that other codes and methods that integrate/average all matrix elements in the entire Brillouin one do not provide. Clearly, we can identify the strength of intra-valley versus inter-valley deformation potential scattering, as well as the relevant strength between the different inter-valley transitions themselves (see Fig. 5c). This capability is essential to evaluate the potential benefits of the band alignment strategy for novel TE materials to improve PF, one of the main directions that the community takes to enhance the materials’ performance.

Fig. 5: Comparison of intra-valley and inter-valley scattering in Mg₃Sb₂.

a 2D view of the first Brillouin zone for electrons (white) and phonons (blue), with the q points involved in the calculations depicted by colored dots. The initial state is at the center of the phonon Brillouin zone (purple dots). b Scattering rates for intra-valley deformation potential scattering (ADP + ODP) and three types of inter-valley scattering. c Comparison between the scattering strength of intra-valley and inter-valley deformation potential scattering. In (a), the colored dots refer to the q points used in our method to derive deformation potentials. In (c), the scattering strength is the sum of the scattering rates within the energy range 0–0.1 eV.

It should also be mentioned that our method is based on the deformation potential theory and can be applied to semiconductors and insulators but not to metals. For the electronic conductivity of metals, the usual method is to perform the computation using the Eliashberg function. A formalism that replaces the matrix elements with effective deformation potentials could reduce computation costs to the levels described in this work, as it might not be necessary to compute all matrix elements throughout the Brillouin zone, but we are not aware of developments in that direction.

In conclusion, we have introduced a first-principles approach for calculating deformation potentials and scattering rates, which are subsequently employed in a Boltzmann transport simulator to compute transport properties of semiconductors and insulators. By encompassing all common scattering processes, i.e., ADP, ODP, IVS, POP, and IIS, our method achieves comparable accuracy to state-of-the-art fully ab initio DFPT + Wannier methods with a drastic computational reduction, at least 100 times for simple materials with small unit cells and at least 10 times for more complex materials with larger unit cells and lower symmetry, owing to the limited number of matrix elements required for deriving deformation potentials. Furthermore, our method provides explicit information about all scattering processes individually, even distinguishing between different inter-valley transitions. We anticipate that our approach will offer a versatile alternative, combining efficiency, robustness, and flexibility beyond the commonly employed CRT approximation, while maintaining the accuracy of fully first-principles calculations. This method can prove invaluable for gaining insights into individual transport properties and guiding machine learning studies by facilitating the development of accurate test sets and descriptors.

Methods

DFT calculations

The electronic band structure, phonon dispersion, and e–ph coupling matrix elements are calculated from DFT and DFPT using the Quantum ESPRESSO package⁸⁴. The optimized norm-conserving Vanderbilt (ONCV)⁸⁵ pseudopotentials are used under the generalized gradient approximation (GGA) with the Perdew-Burke-Ernzerhof (PBE)⁸⁶ functional. The EPW package³⁴ is used to perform Wannier function interpolation for the e–ph coupling matrix elements.

Scattering rates

The e–ph scattering rates due to ADP are extracted using Fermi’s golden rule as^57,87,88

$$\begin{array}{r}\left\vert {S}_{{{{\bf{k}}}},{{{{\bf{k}}}}}^{{\prime} }}^{{{{\rm{ADP}}}}}\right\vert =\frac{\pi }{\hslash }{D}_{{{{\rm{ADP}}}}}^{2}\frac{{k}_{B}T}{\rho {v}_{{{{\rm{s}}}}}^{2}}g(E)\end{array}$$

(5)

Here, ρ is the mass density, and g(E) is the density of states for the final scattering state.

The corresponding ODP scattering rates are computed as^87,88

$$\begin{array}{r}\left\vert {S}_{{{{\bf{k}}}},{{{{\bf{k}}}}}^{{\prime} }}^{{{{\rm{ODP}}}}}\right\vert =\frac{\pi {D}_{{{{\rm{ODP}}}}}^{2}}{2\rho \omega }\left({N}_{\omega }+\frac{1}{2}\mp \frac{1}{2}\right)g(E\pm \hslash \omega )\end{array}$$

(6)

where ω is the dominant frequency of the optical phonons, which is considered to be constant over the entire reciprocal unit cell. N_ω is the phonon Bose-Einstein statistical distribution and the + and − signs indicate the emission and absorption processes, respectively.

The IVS rates follow the optical phonon description and are again computed as

$$\begin{array}{r}\left\vert {S}_{{{{\bf{k}}}},{{{{\bf{k}}}}}^{{\prime} }}^{{{{\rm{IVS}}}}}\right\vert =\frac{\pi {D}_{{{{\rm{IVS}}}}}^{2}}{2\rho \omega }\left({N}_{\omega }+\frac{1}{2}\mp \frac{1}{2}\right)g(E\pm \hslash \omega )\end{array}$$

(7)

Note that ω here is the phonon frequency associated with the corresponding inter-valley scattering processes.

The POP scattering rates due to the Fröhlich interaction can be computed from the dielectric constants, which capture the matrix element in a polarizable continuum as^87,88

$$\begin{array}{r}\left\vert {S}_{{{{\bf{k}}}},{{{{\bf{k}}}}}^{{\prime} }}^{{{{\rm{POP}}}}}\right\vert =\frac{\pi {e}^{2}\omega }{| {{{\bf{k}}}}-{{{{\bf{k}}}}}^{{\prime} }{| }^{2}{\varepsilon }_{0}}\left(\frac{1}{{k}_{\infty }}-\frac{1}{{k}_{0}}\right)({N}_{\omega }+\frac{1}{2}\mp \frac{1}{2})g(E\pm \hslash \omega )\end{array}$$

(8)

where e is the electronic charge, and ω is the dominant frequency of polar optical phonons over the whole Brillouin zone, which has been validated to be a satisfactory approximation⁸³. k₀ is the static dielectric constant.

The Brooks-Herring model⁸⁹ is used to describe the elastic scattering rate due to ionized dopants. The screened IIS is then given by

$$\begin{array}{r}\left\vert {S}_{{{{\bf{k}}}},{{{{\bf{k}}}}}^{{\prime} }}^{{{{\rm{IIS}}}}}\right\vert =\frac{2\pi }{\hslash }\frac{{Z}^{2}{e}^{4}}{{k}_{0}^{2}{\varepsilon }_{0}^{2}}\frac{{N}_{{{{\rm{imp}}}}}}{{(| {{{\bf{k}}}}-{{{{\bf{k}}}}}^{{\prime} }{| }^{2}\,+\,\frac{1}{{L}_{D}^{2}})}^{2}}g(E)\end{array}$$

(9)

where Z is the electric charge of the ionized impurity, N_imp is the density of the ionized impurities, and L_D is the generalized screening length defined as

$$\begin{array}{r}{L}_{D}=\sqrt{\frac{{k}_{0}{\varepsilon }_{0}}{e}\left(\frac{\partial {E}_{{{{\rm{F}}}}}}{\partial n}\right)}\end{array}$$

(10)

where E_F is the Fermi level and n is the carrier density.

Transport

Transport calculations with the DFPT + Wannier method were conducted using the EPW code³⁴, while our method employed the ElecTra code⁴⁸. The k-point mesh for the band structure of Mg₃Sb₂ was initially set to 60 × 60 × 40 and later interpolated to 180 × 180 × 120 (see Supplementary Figs. 10, 11). The kernel quantity, known as the transport distribution function (TDF), is defined as

$$\begin{array}{r}{\Xi }_{ij}(E)=\frac{s}{{(2\pi )}^{3}}\mathop{\sum}\limits_{{{{\bf{k}}}},n}{v}_{i({{{\bf{k}}}},n)}{v}_{j({{{\bf{k}}}},n)}{\tau }_{i({{{\bf{k}}}},n)}g(E)\end{array}$$

(11)

where s is the spin degeneracy and s = 2 is used as the two spin sub-bands are degenerate in the materials we consider here, v_i(k, n) is the i-component of the band velocity, and τ_i(k, n) is the overall relaxation time which is derived from the transition rate $| {S}_{{{{\bf{k}}}},{{{{\bf{k}}}}}^{{\prime} }}|$ between the considered initial state and all the possible final states⁹⁰. The overall $| {S}_{{{{\bf{k}}}},{{{{\bf{k}}}}}^{{\prime} }}|$ is calculated by combining the strength of all scattering mechanisms using Matthiessen’s rule as⁹¹

$$\begin{array}{r}| {S}_{{{{\bf{k}}}},{{{{\bf{k}}}}}^{{\prime} }}| =\left\vert {S}_{{{{\bf{k}}}},{{{{\bf{k}}}}}^{{\prime} }}^{{{{\rm{ADP}}}}}\right\vert +\left\vert {S}_{{{{\bf{k}}}},{{{{\bf{k}}}}}^{{\prime} }}^{{{{\rm{ODP}}}}}\right\vert +\left\vert {S}_{{{{\bf{k}}}},{{{{\bf{k}}}}}^{{\prime} }}^{{{{\rm{IVS}}}}}\right\vert +\left\vert {S}_{{{{\bf{k}}}},{{{{\bf{k}}}}}^{{\prime} }}^{{{{\rm{POP}}}}}\right\vert +\left\vert {S}_{{{{\bf{k}}}},{{{{\bf{k}}}}}^{{\prime} }}^{{{{\rm{IIS}}}}}\right\vert \end{array}$$

(12)

The TDF is crucial in understanding electronic transport in advanced materials. All quantities defining the TDF, including the scattering relaxation times, band velocities, and density of states, are in principle dependent on energy, momentum, and band. Then, within the electronic BTE, one can define the moments of the TDF as

$$\begin{array}{r}{R}^{(\alpha )}={e}^{2}\int\nolimits_{{E}_{0}}^{\infty }dE\left(-\frac{\delta {f}_{0}}{\delta E}\right)\Xi (E){\left(\frac{E-{E}_{{{{\rm{F}}}}}}{{k}_{{{{\rm{B}}}}}T}\right)}^{\alpha }\end{array}$$

(13)

where f₀ is the equilibrium Fermi distribution. From all relevant transport coefficients, the mobility μ, conductivity σ, Seebeck coefficient S, and electronic thermal conductivity κ_e (see Supplementary Fig. 12), are obtained as

$$\begin{array}{r}\mu =\frac{{R}^{(0)}}{ne}\end{array}$$

(14)

$$\begin{array}{r}\sigma ={R}^{(0)}\end{array}$$

(15)

$$\begin{array}{r}S=\frac{{k}_{{{{\rm{B}}}}}}{e}\frac{{R}^{(1)}}{{R}^{(0)}}\end{array}$$

(16)

$$\begin{array}{r}{\kappa }_{{{{\rm{e}}}}}=\frac{{k}_{{{{\rm{B}}}}}^{2}T}{{e}^{2}}\left[{R}^{(2)}-\frac{{[{R}^{(1)}]}^{2}}{{R}^{(0)}}\right]\end{array}$$

(17)