Carrier scattering considerations and thermoelectric power factors of half-Heuslers

Abstract

The electronic and thermoelectric (TE) transport properties of 13 n-type and p-type half-Heusler alloys are computationally examined using Boltzmann transport. The electronic scattering times resulting from all relevant phonon interactions and ionized impurity scattering (IIS) are fully accounted for using ab initio extracted parameters. We find that at room temperature the average peak TE power factors (PFs) of all materials we examine reside between 5 and 10 mW mK⁻². We also find that the combination of IIS and the long range polar optical phonon (POP) scattering is more influential in determining the electronic transport and PF over all other non-polar phonon interactions (acoustic and optical phonon transport). In fact, the combination of POP and IIS determines the thermoelectric power factor of the half-Heuslers examined on average by about 65%. The results highlight the crucial impact of Coulombic scattering processes (POP and IIS) on the TE properties of half-Heusler alloys and provide profound insight for understanding transport, which can be applied widely in other complex bandstructure materials. In terms of computation expense, the computationally cheaper POP and IIS provide an acceptable first-order estimate of the power factor of these materials, while the non-polar contributions, which require more expensive ab initio calculations, could be of secondary importance.

---

1. Introduction

Thermoelectric materials, capable of converting heat directly into electricity and vice versa, have emerged as promising candidates for sustainable energy solutions, including power generation, refrigeration, and waste heat recovery systems.^1–4 The performance of thermoelectric materials is quantified using the dimensionless figure of merit ZT, defined as ZT = σS²T/κ, where σ is the electrical conductivity, S is the Seebeck coefficient, T is the absolute temperature, and κ is the total thermal conductivity (κ = κ_e + κ_l), i.e., the combined electronic and lattice contributions. The product of σS² is called the power factor (PF) and it is a measure of the power generation of the conversion process. Achieving a high ZT requires a material to exhibit a high electrical conductivity and Seebeck coefficient while retaining a low thermal conductivity. This is a challenging endeavour due to the well-documented adverse interdependence of these parameters.³

In the last few decades, thermoelectric (TE) research has focused on reducing the thermal conductivity of materials, using methods such as nanostructuring,^5–9 exploitation of anharmonic bonding within crystal structures,^10–14 boundary scattering,¹⁵ etc. These methods have proved successful in improving ZT from ZT = 1 for only a handful of materials to ZT > 2 for many materials over many different temperatures.^16,17 The thermal conductivities of many prominent thermoelectric materials have now reached very low values, approaching the amorphous limit and even below in many cases. This, however, makes it difficult to achieve further ZT improvements by reducing the thermal conductivity. In recent years, efforts have been equally directed towards improving the thermoelectric power factor, which turns out to be a rather complex task due to the well-known adverse interrelation between σ and S with carrier density.

The majority of strategies to improve the power factor revolve around band engineering, such as band convergence,^16,18,19 alloying,²⁰ introduction of specific orbital interactions to introduce band-gap in metallic systems and increase the valley-degeneracy of the conduction band,^21,22 increasing the entropy of the crystal structure to increase the electronic band convergence,²³ introduction of topological bands to increase σ,^24,25 enhanced inter-band scattering that increases S,^26–28 etc. Among these, band convergence is one of the most widely adopted strategies, as it increases the number of transport channels.^{16,19,22,29,30} This improves the electrical conductivity and carrier density at the band edges, which can also enhance the Seebeck coefficient, ultimately leading to higher PFs.^31–34

For optimizing the PF of complex bandstructure materials, an understanding of electronic transport, and more specifically of the electronic scattering processes, is required. Electrons in general scatter with acoustic phonons, optical phonons, polar optical phonons, ionized impurities, etc. For this, accurate computational methods are pivotal in understanding, designing, and optimizing thermoelectric materials. A hierarchy of methods with increasing accuracy (and computational cost) is available for the electronic scattering times calculation. The most widely used method considers the constant relaxation approximation (CRTA) as implemented in BoltzTrap,^35,36 which, however, introduces an arbitrary quantitative error in the electrical conductivity and PF.³⁷ On the other side of the spectrum, first-principles calculations, particularly those based on density functional perturbation theory (DFPT) together with Wannierization to describe the electron-phonon interactions, provide high accuracy and predictability.³⁸ However, these come at the expense of enormous and difficult-to-scale computational costs, as they involve the calculation of millions of matrix elements.³⁸ Intermediate methods have also recently been developed that require far fewer matrix elements compared to full ab initio methods, which they later post-process to extract effective scattering rates.^39,40 Other methods use deformation potential theory with the use of matrix elements,⁴¹ or without the use of matrix elements.⁴² While accuracy in the quantification of electron-phonon interactions is considered essential for understanding electronic and thermoelectric transport properties, any method that extracts matrix elements is computationally expensive and hard to scale, thus less practical in material optimization studies. Quantifying the strength of the different mechanisms and their influence on TE performance can help prioritize computational tasks for fast evaluation of materials at the level of being able to drive machine learning studies as well.^43–45

Among the emerging materials for thermoelectric applications, half-Heusler (HH) alloys are promising candidates, due to their stability, abundance, and high intrinsic power factors (PFs), which compete with or even surpass those of some of the best traditional TE materials.^3,46–49 The reason for their high performance lies in their complex, multi-band, multi-valley electronic structures, which contain multiple carrier pockets, leading to large conductivity and Seebeck coefficients. Half-Heusler alloys have a general formula of XYZ, where X and Y are typically transition or rare earth metals and Z is a main group element. They crystallize in a cubic MgAgAs (C_1b) structure, characterized by a face-centered cubic lattice composed of four interpenetrating sub-lattices. The versatility of HH alloys allows for tunable electronic^43,50 and thermal properties^49,51,52 through elemental substitution and doping. This provides excellent opportunities for further optimization of their pristine thermoelectric performance through bandstructure engineering for PF improvements and alloy scattering for thermal conductivity reduction. Band alignment to increase σ is an important direction of research in HHs (and full-Heuslers)^46,53–55 as well as the exploration of topological bands,⁵⁶ flat bands⁵⁷ and carrier filtering from flat bands that increases S,²⁸ broken bands that also increase S,⁵⁸ etc.

However, in order to understand and further optimize the PF in these materials, we need knowledge of the full details of the processes that determine electronic transport and how the specific bandstructure features determine those transport properties. For this we typically employ the Boltzmann transport equation (BTE) using relevant solvers, but accurate descriptions of the scattering rates including their full energy/momentum/band-dependences are required. In a typical pristine TE material scenario, i.e., before alloying, the relevant processes that electrons scatter off are acoustic phonon (ADP), non-polar optical phonon (ODP), polar optical phonon (POP), and ionized impurity scattering (IIS). While we know that IIS is a strong mechanism in most TEs,^37,59–64 not much is known about the relative strength of all these scattering mechanisms. POP is also suspected to be a strong mechanism because of its Coulombic nature, but its strength has not yet been quantified. For example, we know that in materials such as Mg₃Sb₂, POP is as strong as, if not stronger than, IIS and overshadows all other non-polar mechanisms.⁶⁵ Accurately quantifying these contributions is vital for understanding material behavior and guiding the development of alloys with enhanced ZT values. For example, materials in which IIS and POP dominate have the advantage of PF benefits upon band alignment when valleys are far from each other in the Brillouin zone, since these mechanisms are strongly anisotropic and favor small momentum exchange vectors.^60,66–68 In addition, knowledge about the strength of each mechanism can prioritize the considerations of computational studies; i.e. non-polar scattering times require ab initio treatment, and if they are significantly weaker, their extraction could be avoided with large computational savings at a small expense in accuracy.

Here, we quantitatively investigate electronic and thermoelectric transport in 13 n-type and p-type half-Heusler materials at room temperature. The work is organized with a focus on (i) providing transport-related details for HHs such as deformation potentials, the nature of overlap integrals and intra-/inter-valley transitions. (ii) Establishing the influence of the Coulombic IIS + POP scattering mechanisms versus that of the non-polar ADP + ODP phonon mechanisms. We show that, put together, IIS and POP determine the PF performance more strongly on average by about 65%, with a stronger influence on electrons compared to holes. (iii) Providing the PFs of these materials as accurate and as computationally efficient as possible. In section II, we discuss the computational methodology we follow. In section III we discuss the results with a focus on the strength of POP compared to that of non-polar phonon processes and their effect on the PF. Section IV compares the electron-phonon scattering strength with that of IIS and provides full PF results. Finally, Section V summarizes and concludes the work.

2. Methodology

We consider a group of 13 HHs and perform calculations for both their n-type and p-type polarities. We use density functional theory as implemented in the Quantum Espresso package to calculate the electronic structures. We use optimized norm-conserving Vanderbilt (ONCV) based pseudopotentials. For exchange–correlation, we use the Perdew–Burke–Ernzerhof (PBE) functional based generalized gradient approximation (GGA) over the local density functional. This approach is validated for these half-Heusler compounds by the good agreement of our calculated lattice parameters with the literature^69–72 (see Table 1). A cut-off of 120 Ry for wave-function expansion and an energy convergence criterion of 10⁻⁸ for self-consistency are employed. A dense mesh of 81 × 81 × 81 is used for non-self-consistent calculations of the band structure to increase the accuracy of the transport calculations. We have used our BTE simulator ElecTra for the transport properties.⁴¹

Table 1 Scattering parameters used in the calculation of the transport coefficients for each of the materials considered. Equilibrium values of lattice constant (a₀), values for acoustic deformation potentials (D_ADP), optical deformation potentials (D_ODP), and inter-valley deformation potentials (D_IVS) are shown. ℏω_ODP and ℏω_POP are the highest energies of the transverse optical (TO) and longitudinal optical (LO) phonon modes at the Γ point. k₀ and k_∞ are the static and high frequency dielectric constants. Their ratio in the last column provides the splitting of the LO-TO modes and an indication of the strength of the polar dipole interaction (1 means zero interaction, while 2 signifies a very polar material)

Material	a0	DADP (eV)		DODP (eV Å^−1)		DIVS (eV Å^−1)		ℏωODP	ℏωPOP	Dielectric constants
	(Å)	n-type	p-type	n-type	p-type	n-type	p-type	(eV)	(eV)	k0	k∞	k0/k∞
ZrNiSn	6.15	2.67	2.70	2.12	1.34	0.86	—	0.029	0.032	27.19	21.95	1.24
HfNiSn	6.15	1.87	3.32	2.09	1.40	1.02	—	0.028	0.030	25.94	20.94	1.24
TiCoSb	5.96	3.19	2.29	3.23	1.34	0.80	—	0.036	0.037	31.84	20.92	1.52
ZrCoSb	6.09	2.43	2.37	2.07	1.35	1.02	1.24	0.028	0.033	27.56	18.87	1.46
HfCoSb	6.07	2.35	2.24	3.64	2.07	1.17	1.1	0.028	0.033	25.72	18.12	1.42
ScNiSb	6.12	3.29	2.42	2.56	1.93	1.14	—	0.027	0.030	22.86	18.79	1.22
YNiSb	6.37	3.72	2.72	1.44	1.89	1.22	—	0.023	0.025	21.72	18.82	1.15
ScNiBi	6.27	3.20	2.33	2.28	1.67	0.91	—	0.025	0.027	26.95	22.86	1.14
YNiBi	6.49	2.63	2.68	1.19	1.62	1.17	—	0.021	0.023	26.23	23.53	1.11
ZrCoBi	6.22	3.34	2.98	2.37	2.03	0.86	1.22	0.025	0.031	29.60	20.71	1.43
NbFeSb	6.00	3.72	2.50	2.03	2.80	0.77	0.6	0.036	0.042	42.84	24.95	1.72
NbCoSn	5.97	3.44	2.04	2.10	1.85	1.26	2.07	0.030	0.037	36.27	24.76	1.46
ZrNiPb	6.25	2.32	1.90	1.65	1.49	0.86	—	0.024	0.027	29.46	23.91	1.23
Average (not used)	—	2.94	2.50	2.21	1.67	1.0	1.25	0.026	0.031	28.79	21.47	1.33

---

Standard expressions for the transition rates of the different scattering mechanisms are considered. For acoustic deformation potential (ADP) we use

where ρ is the mass density of the material and g_k′(E) is the final density of states for the scattering event (k → k′).

For optical deformation potential (ODP) transition rates we use

where ω is the frequency of the optical phonon, which is considered constant in the entire Brillouin zone. N_ω stands for the Bose-Einstein statistical distribution function for phonons. The symbols ‘+’ and ‘−’ indicate the emission and absorption processes. Inter-valley scattering (IVS) needs to be accounted for, as in the case of n-type carriers for scattering between X-valleys (the CBM of most HHs is at the X high symmetry point) and some p-type materials where the VBM is at a high symmetry point other than Γ, such as for NbFeSb at the L-point, NbCoSn at the L- and W-points, etc. The strength of the inter-valley scattering contribution is evaluated using all the phonon modes as discussed in our previous studies.^73–75 Scattering due to IVS is evaluated as

In this work, we combine both ODP and IVS, so hereafter, ODP scattering refers to the combined contribution from both ODP and IVS.

We computed the values used for D_ADP and D_ODP using matrix elements from density functional perturbation theory, followed by Wannierization as in ref. 38 and 73–75. While computing the transitions for all phonon branches and directions separately, for easiness we combine all values to a ‘global’ one for D_ADP, one for D_ODP and one for D_IVS for each material, i.e. one single value applied to all transitions between all bands/valleys. These values are given in Table 1 for all materials. An example of how we perform the matrix element calculations is shown in the SI in Fig. S1–S3.

For the calculation of the polar optical phonon (POP) scattering, we use the Fröhlich formalism as⁷⁶

where e is the electronic charge and ω is the dominant frequency of polar optical phonons. ε₀ is the free space permittivity and k_∞ and k₀ are the static and high-frequency dielectric constants, respectively. is the square of the wavefunction overlaps, which is, in principle, band, energy and momentum dependent. We have computed multiple combinations of these from a large number of initial states to a large number of final states in the Brillouin zone and show these values in Fig. 1 versus their separation in the BZ for both intra-valley (left column) and inter-valley (right column) transitions. These values vary from 0 to 1. The first row shows data for the three valleys in the VB of HfNiSn, where the averaged values for the squared overlaps for both intra- and inter-valley transitions are around 0.33. The second row shows data for the VB of NbFeSb where two bands are present, and the averaged values are around 0.5. These values are well understood and explained in ref. 75. For simplicity, we used these averaged values in the scattering rates across all states for each relevant material (0.33 for three-fold, 0.5 for two-fold, and 1 for single degenerate bands).

Fig. 1 (a and b) Square of the wavefunction overlaps between the states in the three-fold degenerate valence bands of HfNiSn at the Γ-point for intra-band and inter-band transitions, respectively. (c and d) Square of the wavefunction overlaps between the two-fold degenerate valence bands of NbFeSb at the L high symmetry point for intra-band and inter-band transitions. The insets depict the bands and the transitions with the same color dots.

To account for the effect of screening of the polar dipole interaction by the accumulation of charge carriers, the additional term is also included (we multiply the Fröhlich expression above by this term). Here, is the generalized screening length, where E_F is the Fermi level and n is the carrier density.⁶⁸ This additional term makes the calculation more costly since the screening length depends on density and the Fermi level, but it is necessary since the PF peaks at large densities and screening will be important.⁶⁸ After multiplication of the screening term, the final scattering due to POP becomes

To account for ionized impurity scattering (IIS) we use the Brooks-Herring model as⁷⁷

where Z stands for the electric charge of the ionized impurity (Z = 1 is used in the present work) and N_imp is the density of the ionized impurities. To calculate the total scattering transition rate, we use Matthiessen's rule as

After calculating the transition rates, the thermoelectric coefficients, namely the electrical conductivity, σ, and Seebeck coefficient, S, are calculated as

Here, f₀ is the equilibrium Fermi–Dirac distribution function and Ξ_ij(E) is the transport distribution function (TDF) defined as

where v(E) is the bandstructure velocity, g(E) is the density of states at energy E, and i, j are the Cartesian coordinate indices, for which we set i = j = x), and τ_k,k′ is the scattering relaxation time (inversely proportional to the transition rates), which is calculated for each scattering mechanism (m′s) as

We start our analysis by examining the main features of a typical HH material by considering HfNiSn as an example. First, we examine the effect of electron-phonon scattering processes on transport, as in this part of the paper we focus on investigating the effect of POP compared to all other non-polar mechanisms. We consider IIS later on in the manuscript. The band structure for HfNiSn is shown in Fig. 2(a). The conduction band consists of a single degenerate band with its valley minimum (CBM) placed at the X high symmetry point. A three-fold degenerate valence band is observed, with its valence band maxima (VBM) located at the Γ high symmetry point. These features agree with other reports in the literature.^78,79 Fig. 2(b)(i) and (ii) illustrate, respectively, the Fermi surface of HfNiSn at energy E = 60 meV above the CBM, indicating the six-fold degenerate X-valley, and at 60 meV below the VBM, indicating the three degenerate bands at the Γ point. After calculating the band structure, we use the Boltzmann transport code ElecTra for the calculation of the TDF and the thermoelectric properties.

Fig. 2 (a) The band structure of HfNiSn. (b) (i) and (ii) show the Fermi surface of the conduction and valence bands, respectively, at energy E_C + 60 meV and E_V − 60 meV. (c and d) The transport distribution functions (TDFs) for n-type and p-type for different phonon scattering-limited transport considerations as noted in the legend. (e and f) The scattering times (τ) for n-type and p-type HfNiSn, respectively for the scattering cases in (c and d).

3. Phonon scattering and the power factor

Fig. 2(c) and (d) show the calculated TDFs versus energy for both n- and p-type carriers, respectively. We compute these separately as unipolar systems, i.e. we don't consider bipolar effects, because we perform all simulations at room temperature T = 300 K and because we focus more on examining the isolated band structure and quantifying the scattering processes. Here, we plot separately and compare the TDF from the combined non-polar scattering mechanisms, i.e. ADP and ODP (blue lines), the polar contributions alone (green lines), and finally the total TDF from all phonon mechanisms (black line). Among the ADP + ODP, for n-type both ADP and ODP have similar contributions of relaxation time, whereas for p-type the major contribution comes from the ADP scattering process, as shown in the SI, Fig. S4. The kinks in the TDFs at around 30 meV indicate the starting point of the phonon emission of the inelastic ODP and POP scattering process (the kink in the blue line is overshadowed by the elastic ADP process). In the n-type case, the contribution to the TDF from the POP (green lines) is very close to the total TDF (black line) that considers all phonon scattering mechanisms, indicating that POP has more influence in determining electronic transport. The TDF from the combined contribution of ADP + ODP is much higher, indicating significantly higher scattering times compared to POP. The overall scattering rate is computed using Matthiessen's rule, and thus the process with the largest rate (smallest time) dominates. In the case of p-type, on the other hand, the trend is reversed, with the POP being the weaker mechanism in determining transport. This is further corroborated by the comparison of the relaxation times (τ) for both n-type and p-type polarities, as shown in Fig. 2(e) and (f), respectively, which follow the same trend. Here, we would like to mention that POP is inclusive of screening, unless otherwise mentioned, which, as we will show, has important implications (even though it increases computational cost).

The HfNiSn calculated TE coefficients versus reduced Fermi level η_F, defined as E_F – E_C/V for the conduction/valence bands, are shown in Fig. 3 (left and right columns, respectively) at T = 300 K (the temperature considered throughout this work). In each case, we assume that the energy at the band edge (either valence or conduction) is set to η_F = 0 eV, while η_F > 0 eV values indicate that the Fermi level is placed in the bands, as implemented in ElecTra.⁴¹ For both n-type and p-type polarities, it is evident from Fig. 3(a) and (d) that the electrical conductivity (σ) follows the TDFs. POP is the strongest of the phonon mechanisms (green line closer to the overall black line) for n-type, especially at lower η_F values, and around η_F = 0 eV, where the PF peaks. The non-polar phonon scattering contributions are stronger for p-type in the entire Fermi level range (blue line closer to the black line). This trend is typical for most HHs examined.

Fig. 3 Thermoelectric coefficients for HfNiSn for n-type (left column) and p-type (right column) carriers. Three phonon scattering limited transport situations are presented: ADP + ODP, POP, and all phonon scattering. (a and d) Electrical conductivity. (b and e) Seebeck coefficient. (c and f) Power factor. Insets of (c) and (f): Comparison between the power factor computation results with (solid) and without (dashed) screening in the POP calculation.

In the case of the Seebeck coefficient (S) (Fig. 3(b) and (e)), for n-type the values of polar and non-polar contributions are essentially identical, while for p-type the POP-limited result is slightly lower (as is expected under higher conduction). Overall, the n-type and p-type Seebeck coefficients are also similar in absolute terms. This is also a general trend.

The PFs resulting from the combinations of the electron-phonon scattering mechanisms considered are shown in Fig. 3(c) and (f) for the two polarities. In the case of considering only the non-polar (ADP + ODP) scattering mechanisms, the peak PF values are 29.16 mW mK⁻² and 26.73 mW m⁻¹ K⁻² for n-type and p-type carriers, respectively. In the case of POP-limited transport, the trend is different, with the PF being 28.49 mW mK⁻² and 40.85 mW m⁻¹ K⁻² for n-type and p-type carriers, respectively. The overall PF peaks are further reduced when combining both polar and non-polar scattering mechanisms to 12.60 mW mK⁻² and 15.97 mW mK⁻² for n-type and p-type, respectively. The PF values comprising all phonon scattering mechanisms are found to be closer to the POP-limited values in n-type and to the ADP + ODP values in p-type, indicating what has more influence on performance in each case.

The insets of Fig. 3(c) and (f) demonstrate the importance of screening in POP. They compare the PFs by considering all phonon scattering mechanisms for the cases where the POP scattering rate includes (solid line) and excludes (dashed line) the screening term. When screening is not included in the scattering rate calculation, scattering is stronger, leading to a 41% and 38% PF decrease in n-type and p-type materials, respectively, in the region of maximum PF. However, for n-type carriers, around η_F = 0.1 eV and higher, the PF ordering of the two computation treatments changes, with the influence of the unscreened treatment increasing compared to the screened one. This is mainly due to the higher value of the Seebeck coefficient under the unscreened POP treatment, as compared to the S of the screened POP (stronger scattering in general results in larger S). Thus, at higher values of η_F, the PF is governed by the higher values of the Seebeck coefficient, as shown in the SI in Fig. S5.

After presenting an example of a HH material and demonstrating the stronger role of POP for n-type and the stronger role of the non-polar ADP + ODP for p-type, we will proceed to conduct a similar investigation on the list of 12 additional half-Heusler alloys to derive broader insight. Fig. 4 illustrates the PF for both n-type (first row) and p-type (second row) carriers in these HHs, which are noted in the legend of Fig. 4(d). We consider the same three scattering cases as above for HfNiSn. Essentially, the first column presents simulations where scattering is limited by ADP + ODP, the second column presents simulations where scattering is limited by POP, and the third column presents simulations where the scattering includes all non-polar and polar contributions (ADP + ODP + POP).

Fig. 4 The power factor versus density for all 13 HHs considered, for n-type (first row, a–c) and p-type (second row, d–f) polarities. The simulation results for three phonon scattering limited transport considerations are shown column-wise: (a and e) non-polar ADP + ODP scattering, (b and d) POP scattering, and (c and f) scattering due to all phonon processes combined (ADP + ODP + POP).

The average PF values in each case at η_F = 0 eV, near the Fermi level position for maximum PF, are all indicated in each sub-figure (individual PF values are shown in Fig. S6). In the case of n-type, as shown for the ADP + ODP scattering limited transport in Fig. 4(a), the average PF value is 29.16 mW mK⁻². In the case of n-type POP-limited transport, the average PF value as shown in Fig. 4(b) is 28.49 mW mK⁻². This indicates that for electrons, overall, the POP has a similar influence on the PF as the ADP + ODP (HfNiSn earlier was an outlier from this general trend). Fig. 4(c), for the PF of the materials under full phonon considerations, shows an average PF value of 12.60 mW mK⁻². In the n-type case, since the POP and ADP + ODP have a similar influence on the PF, the overall phonon-limited PF is around half from the POP-limited one (and ADP + ODP).

For the p-type materials, under ADP + ODP as shown in Fig. 4(d), the average PF value is 26.73 mW mK⁻². This is a similar value to n-type materials under the same ADP + ODP. Despite slightly lower deformation potential values for p-type (as shown in Table 1), the overall averaged PF value is determined by many complex factors. In the case of POP-limited transport, the average p-type PF value as shown in Fig. 4(e) is 40.85 mW mK⁻². This is a higher value compared to both p- and n-type ADP + ODP and n-type POP-limited transport as well, indicating the weaker overall POP influence on holes in these 13 HHs, following the earlier discussion for HfNiSn. The difference between n-type and p-type is a consequence of the different bandstructure shapes between the CB and VB (since the other parameters that control POP, the dielectric constant and phonon energies, are the same for both carriers). We will be discussing this further below. Fig. 4(f), for the PF of the p-type materials under full phonon considerations, shows an average PF value of 15.97 mW mK⁻². This is very similar (only slightly higher) to the corresponding n-type value, reflecting the fact that POP is weaker for holes.

Finally, notice the exceptionally high value for p-type NbCoSn (purple squared line) under ADP + ODP due to its high degeneracy from L- and W-valleys and lower D_ADP value. This larger advantage is lost under POP, however, due to the smaller distances between the various W-valleys in the BZ, which reduce the POP exchange vectors (and the fact that POP is independent of D_ADP). In contrast, under POP, HfCoSb with L- and Γ-valleys performs higher, as transitions between these valleys involve larger exchange vectors. This advantage is also lost when all phonons are considered, but still NbCoSb remains the best PF performing p-type material (phonon-limited) and HfCoSb also performs well compared to others.

To better understand the PF trends, and especially the comparison between the strengths of POP and the non-polar ADP and ODP mechanisms, in Fig. 5(a and b), we show the momentum relaxation scattering times for all 13 HHs under consideration for these three mechanisms separately. Fig. 5(a) shows the scattering times for n-type and Fig. 5(b) for p-type carriers, respectively. ADP relaxation times are shown by the red lines, ODP relaxation times by the blue lines, and POP relaxation times by the green lines. Note that here we computed the POP rates for the case where the Fermi level is aligned to the band edge, i.e., η_F = 0 eV, where the PF typically peaks (the POP screening depends on density). Also note that our purpose here is to reach general conclusions, rather than perform analyses for individual materials; thus for simplicity we don't label each material, but we group the different mechanisms together with different colors.

Fig. 5 Scattering relaxation times for different mechanisms and their effect on transport for all 13 HHs considered. (a and d) The scattering times under ADP- (red lines), ODP- (blue lines), and POP-limited electron-phonon scattering conditions. (a) Results for n-type carriers and (d) for p-type carriers. (b and c) The averaged value of the relaxation times for each category in (a and d). (c) Scattering times for the POP-limited case for n-type (open circles) and p-type (full circles). (f) The averaged values of the TDFs due to ADP- (red), ODP- (blue), and POP (green)-limited electron-phonon scattering for n-type (open circles) and p-type (full circles) carriers.

To better illustrate the strength of the different mechanisms, in Fig. 5(b) and (e), we calculate the average relaxation times for the three phonon mechanisms for n-type and p-type carriers, respectively (i.e. we simply averaged the results of Fig. 5(a) and (b)). Despite the large spread in the individual materials' data, the averaged values clearly indicate an overall trend, with the POP being the strongest mechanism for n-type compared to ADP and ODP (lower green line), although not so drastically, so the case of ADP + ODP has a similar influence on the PF. In the p-type case (Fig. 5(e)), the POP is similar to the ADP and ODP individually, such that when the latter two combine, they become stronger.

We will now compare the averaged POP rates for n-type and p-type in Fig. 5(c), to quantify how much n-type materials are affected more by POP compared to p-type materials (empty-symbol line for electrons vs. filled-symbol line for holes). Overall, the relaxation time for n-type carriers is about 3× lower than that of p-type carriers, indicating stronger POP influence for n-type overall, which is also reflected in their average TDF values in Fig. 5(f) and the PF values in Fig. 4. The ratio between the average PF of p-type carriers to that of n-type carriers when only POP scattering is considered is 1.43 (40.85/28.49), indicating through that comparison that it is more complicated than just comparing scattering times – we don't refer to the Seebeck coefficient much here because it is very similar for n-type and p-type under POP scattering conditions. Also note that on average, other than for the case of POP, the influence of the ADP and ODP on the TDF is very similar for n-type and p-type materials (red and blue lines in Fig. 5(f)).

While we showed that the POP scattering rates for electrons are more substantial compared to those of holes by direct calculation, it is useful to provide illustrative explanations based on bandstructure features. The CB of HHs typically consists of a single band whose valleys reside at degenerate high symmetry points (i.e., other than Γ), at least for most of the systems we consider (as seen in the band structures in the SI file, Fig. S7). It is also typical for the valence bands of HHs to consist of degenerate bands of different effective mass, but at the same high symmetry point or high symmetry points. By examining the expression for the POP scattering rate, the differences that different band structures bring are related to (i) the screening length L_D and (ii) the exchange vector q. The smaller the screening length is, the more reduced the POP dipole interaction's strength is. This is determined by the density of the material, which in turn is determined by the DOS. The DOS in general is higher for p-type channels (see the extracted m_DOS in the SI in Fig. S8). The L_D values for all n-type and p-type materials we examine are shown in Fig. S9 in the SI file. In Fig. 6(a) below, however, we show the averaged L_D values for n-type materials (empty symbols) and p-type materials (filled symbols) versus the reduced Fermi level. As η_F increases, the carrier density increases, which reduces L_D. The p-type L_D, however, is lower than that of the n-type in the entire range by more than 1.4× (ratio is shown in the inset of Fig. 6(a), which indicates that POP will be weaker for holes compared to electrons. Note that this also applies to IIS further below.

Fig. 6 (a) The averaged screening length for all 13 HHs considered versus Fermi energy for n-type (open circles) and p-type (full circles) carriers. Inset: the ratio of the two quantities. (b) Zoom-in for the band structure of HfNiSn, also showing the energy surfaces of the different bands (labelled V_i) at 60 meV in the band extrema. (c) Same as in (c) for NbFeSb.

In addition, the larger the q-vector is, the weaker the scattering, since POP scattering is an anisotropic mechanism that favors small-angle scattering, i.e., it is inversely proportional to the exchange q-vector between the initial and final scattering states. The most detrimental to transport are the vectors that can cause backscattering (or the larger the scattering angles), but those in the case of POP (and IIS) have the weakest scattering rates. Thus, the size of the energy surfaces can provide a measure of the strength of POP, with larger energy surfaces having larger q-vectors, and reduced backscattering. Fig. 6(b and c) show the band structures of HfNiSn and NbFeSb as generic examples, together with the energy surfaces for the CBs and VBs at 60 meV in the bands. The VB energy surfaces are larger, indicating that larger q-vectors are associated with the scattering events, leading to reduced POP scattering. Note that the three bands at the VB edge are degenerate at the Γ-point for HfNiSn (twofold for NbFeSb at the L-point), but the degeneracy is lifted as the energy increases (in the negative direction) and this increases the energy surface size.

4. IIS strength and the overall power factor

In addition to electron–phonon scattering, ionized impurity scattering (IIS) is also an important scattering mechanism, which is typically very strong due to its Coulombic nature. We now investigate the role of IIS in the thermoelectric properties of the half-Heusler alloys under consideration, with the goal not only to provide information about their ultimate PF values, but also to provide a comparison of IIS to the strength of POP and the non-polar phonon scattering mechanisms ADP + ODP.

Fig. 7(a) and (d) show a comparison between the relaxation scattering times for POP (green lines) and IIS (purple lines) for n-type and p-type carriers, respectively, for all 13 HHs under consideration. Since both mechanisms involve screening, which depends on the charge density, we are performing this comparison here specifically for η_F = 0 eV. In the low-energy region, IIS is stronger compared to POP scattering for the materials considered, leading to lower relaxation scattering times, due to the well-known divergence of IIS rates at low energies. The trend changes with the POP rate becoming increasingly stronger at higher energies with lower scattering times overall, while IIS times increase with energy. When comparing the screened POP scattering rate eqn (5) with the IIS rate eqn (6), POP has an additional q-vector squared (k − k′)² term in the numerator. The exchange vector, in general, increases at higher energies and leads to higher POP scattering rates compared to IIS (lower scattering times for POP). Comparing the n-type to the p-type cases, IIS is a stronger mechanism for n-type carriers across energies, again due to larger screening lengths L_D,n and narrower energy surfaces and smaller scattering exchange vectors, q, as in the case of POP as well. Thus, both POP and IIS are weaker in the VB compared to the CB.

Fig. 7 (a and d) The momentum relaxation scattering times for POP (purple lines) and IIS (green lines) for n-type and p-type carriers, respectively. All 13 HHs under consideration are shown indiscriminately. (b and e) The PF of n-type and p-type carriers, respectively, under IIS-limited scattering considerations for the group of all 13 HHs. (c and f) The overall PF of all n-type and p-type materials under all phonon + IIS scattering considerations.

Fig. 7(b) and (e) plot the PF due to IIS limited transport conditions. The averaged values of the PF due to IIS at η_F = 0 eV are 20.12 mW mK⁻² and 21.62 mW mK⁻², very similar for both n-type and p-type carriers, respectively. Despite the IIS rate being slightly higher for n-type, which lowers the conductivity compared to p-type, the n-type IIS-limited Seebeck coefficient is slightly higher, and the averaged PFs are similar. These data are shown in the SI in Fig. S10. The averaged PF values here are lower compared to the POP-limited transport PF values in Fig. 4, also indicating that for the Fermi level position where the PF peaks, IIS is a more determining mechanism compared to POP. For p-type this difference between POP and IIS is particularly high, since the IIS rate depends on the impurity density, which is equal to the carrier density here, and is higher in p-type materials due to their higher DOS.

Finally, in Fig. 7(c) and (f), we show the overall PFs of these materials at room temperature for n-type and p-type polarities, respectively, including all phonon processes and IIS. The average PF values at η_F = 0 eV are more than halved now with the inclusion of all phonon processes and IIS to PF_n = 7.04 mW mK⁻² for n-type and PF_p = 8.58 mW mK⁻² for p-type. These differences between n-type and p-type and the differences in PF between materials themselves, at this point, originate from the conductivity. As shown in the SI in Fig. S11, the Seebeck coefficients of all materials (n-type and p-type) are quite similar, with p-type being only 5% lower.

These values are close to those observed in experiments,^80–82 but on the upper level as expected. In general, experimentally measured PFs are lower compared to computed ones due to the multiple defects that exist in the real material and not accounted for in simulations. Experimental samples are typically polycrystalline and include grain boundaries, defects, dislocations, etc. Grain boundaries are common in polycrystalline materials, which scatter carriers and reduce mobility and conductivity. Vacancies can act as scattering centres, affecting the carrier concentration and mobility. Extra atoms occupying interstitial sites can disrupt the crystal lattice and scatter carriers. There are various other factors like antisite defects or secondary phases formed during synthesis, which can lead to additional scattering. These factors almost certainly reduce the electrical conductivity and mobility significantly. In our previous work on NbFeSb,⁷⁵ we found that, overall, the electrical conductivity of the experimental data is around a factor of three lower compared to our computed results. On the other hand, the measured Seebeck coefficient is larger compared to the computed Seebeck coefficient, following the well-expected adverse trend compared to the electrical conductivity. The electrical conductivity suffers much more in the experiment compared to the improvements in the Seebeck coefficient (compared to the simulated data), so the PF is overall lower in the experiment (especially at higher temperatures). Here we provide a comparison for each of the materials we consider in S12 in the SI, using a large number (over 1400) of measured PF data points, which also demonstrate reasonable agreement (on the upper level) between our computed results and measured data.

All the band structures of the materials we consider are plotted in the SI, Fig. S7. On examining the CB of these HHs, we observe similar features across materials with the presence of the conduction band minima (CBM) at the X high-symmetry point. The presence of multiple valleys, far away from each other in the Brillouin zone, makes most n-type HHs favorable for high thermoelectric PFs, since inter-valley scattering is suppressed.⁶⁴ In addition to multiple valleys and deformation potentials, other parameters that influence transport are the conductivity effective masses of these materials m_cond (shown in Fig. S13(a) in the SI, as extracted using the EMAF code⁸³) and scattering parameters, such as the energy of the phonons involved in inelastic scattering, the static dielectric constant (for ionized impurity scattering), the difference between the ionic and electronic dielectric constant contributions (which determines the strength of the POP scattering), etc. (see all corresponding values in Table 1). It appears that the ordering of the PF has a strong correlation with the inverse of m_cond, which determines the carrier velocities, as discussed in ref. 64, rather than the m_DOS (Fig. S13(b)). With regard to the scattering parameters, most of the materials considered have similar values for these properties. However, a few systems, namely YNiBi, ScNiBi, ScNiSb, YNiSb, and HfNiSn, not only exhibit a smaller m_cond, but also exhibit a much smaller difference between their ionic and electronic dielectric constants, resulting in smaller POP scattering (arising from a reduced value of [1/k_∞- 1/k₀] – also see the last column of Table 1 for a measure of the polar strength). Thus, these perform better in terms of PF as n-type thermoelectric materials. On the other hand, materials such as n-type TiCoSb and NbCoSn have higher m_cond, and while they also have two bands near the conduction band edge (as shown in the SI in Fig. S7) they still have lower power factors. Complex interactions between overlapping bands as in these materials are a disadvantage to transport.⁶⁷

In the case of p-type materials, the valence band minima (VBM) also consist of multiple bands with different valley degeneracies and other band shape complexities as well. The materials with X-, L- or W-valleys have multiple degeneracies with electronic states farther from each other in the Brillouin zone and have higher PFs. Note that this is a feature of IIS, which is an elastic and intra-valley mechanism. This might not always be the case, for example, in the case of inelastic-limited processes like POP, which couples higher/lower energy states and can lower S, or IVS-limited processes, which significantly smoothen the advantage of degeneracy. These, together with different deformation potentials and dielectric constants, will make the PF ranking more material-specific. However, IIS-limited transport (with the strongest contribution) is favoured significantly by degeneracy. In general, p-type PF ranking seems to be more degeneracy controlled rather than determined by the lower m_cond.

One characteristic is that for many of the materials we examine, the VBM resides at the Γ-point, which consists of three degenerate bands (of different effective masses). One would think that holes in those bands experience around three times the scattering as well, through inter-valley/band scattering, which would make them behave as if they are single degenerate and have lower performance. This, however, is not the case. As we showed above, the square of the wavefunction overlap integrals that appear in the scattering rate equations (and indirectly in the electron-phonon matrix elements) turns out to be on average 0.33 for the case of three-fold degenerate bands, both for intra- and inter-valley scattering, and on average 0.5 for double degenerate bands.⁷⁵ These valleys at first order behave in terms of transport as if they are independent. Instead of having only intra-band scattering and no inter-band, scattering is intra- and inter-band but each with one-third of the scattering strength. Thus, at first order, what matters for transport is the valley degeneracy, whereas the positioning of the valleys in the Brillouin zone can provide a secondary effect to transport by broadening the energy surfaces, increasing the exchange vectors and affecting POP and IIS as discussed earlier. Nevertheless, still, the material with multiple degenerate bands at the L-point (NbFeSb: 4-fold degenerate) or at L- and W-points (NbCoSn: 10-fold degenerate, 4 from L and 6 from W) tops the PF, while the ones with three-fold degenerate valleys at Γ have weaker PF performance.

In Fig. 8 we show this better, by grouping and averaging the PFs from Fig. 7(f) into different categories. One group of materials is the one which has the VBM only at the Γ-point (blue line), in which case three bands of different masses exist (YNiSb, ScNiSb, ScNiBi, YNiBi, ZrNiPb, TiCoSb, and HfNiSn). These have lower PFs as only three bands participate in transport. The group of materials with the immediately next higher PF have their VBM at the L-point (red line), which has a degeneracy of 4 (ZrCoSb, ZrCoBi, and NbFeSb). Materials with higher valley degeneracy generally exhibit higher power factors. Thus, the next higher PF material is the one with its VBM at both the L and Γ points (HfCoSb). Finally, the highest PF is observed for NbCoSn, with its VBM at both the L and W high-symmetry points, with a combined degeneracy of 10 (4 from L and 6 from W). Although there is some spread in the PFs within each group that can extend into other groups (as shown by the error bars in Fig. S9), what we illustrate here is a general trend that highly degenerate bands provide higher PFs. Note, however, that the variation in the PF is not fully one to one with degeneracy; thus, any increment in degeneracy leads to a smaller PF increment on average. One reason is that the presence of multiple bands/valleys could enhance inter-valley scattering, thus limiting the overall increase in its PF, and this is more pronounced for materials such as NbCoSn for which the W-valleys reside close by in BZ and are affected by smaller exchange vectors. As shown in the SI in Fig. S14, we also perform the grouping for the cases of POP and IIS separately and POP and IIS combined as well, using data from Fig. 4(e) and 7(e),

Fig. 8 The average values of the PF across p-type materials with their VBs having different high symmetry points as noted in the legend, for all phonon + IIS scattering considerations with the error bars showing the spread in the peak PF for each group.

Since the Coulombic IIS and POP are both the strong scattering mechanisms in HHs, it is interesting to examine how much these two alone determine the overall PF performance of the material. This is important for computational studies, since such calculations do not require the heavy computational matrix element extraction that non-polar phonon scattering involves. For this, we calculate the average PF values for all 13 HHs under consideration for transport limited by many combinations of scattering mechanisms. These are shown in Fig. 9(a) and (b) for n-type and p-type materials, respectively. We show the non-polar ADP + ODP limited transport (blue lines), the IIS (pink lines), the POP (green lines), the POP + IIS (orange lines), the all phonon case (black lines), and the all phonon + IIS (magenta lines) limited transport. Indeed, the PF of the IIS-limited transport case is the lowest compared to the ADP + ODP and POP, showing the importance of IIS for both n-type and p-type carriers, but it is weaker compared to the combination of all phonons. Our interest here is also to compare the PF extracted from the Coulombic IIS + POP combined to the total PF where all phonons + IIS are included. The former overestimates the PF by 30% and 37% in the cases of n-type and p-type carriers, respectively (as noted in Fig. 9). The information that POP + IIS overestimates the PF by an average of 35% is important, since large scale computational studies can be performed using such mechanisms alone and still have a reasonable estimation of the PF. On the other hand, non-polar phonon scattering requires the extraction of matrix elements, which is computationally extremely expensive and non-scalable, especially for complex materials.⁶⁵ Alternatively, POP and IIS mechanisms combined determine the thermoelectric PF by around 70% for n-type materials and 63% for p-type materials (approximately 65% overall). It is also interesting to see that in the case of both n-type and p-type materials, most scattering combinations (IIS, POP + IIS, all phonon, and all phonon + IIS) result in very similar PF values.

Fig. 9 The averaged power factor of all 13 HHs considered for multiple scattering case considerations as labelled in the legend. (a) n-type materials. (b) p-type materials. The overestimation of the Coulombic POP + IIS from the total all phonon + IIS cases is indicated.

Finally, we would like to mention the approximations that we employ. We employed the standard PBE-GGA functional without the addition of a Hubbard U correction (DFT+U). This approach is validated for these half-Heusler compounds by the good agreement of our calculated lattice parameters with the literature^69–72 (see Table 1) We mainly focussed on the role of the shape of the electronic structure for n-type and p-type separately and the role of scattering mechanisms, which makes our conclusions independent of the exact bandgap values. In this way, we do not consider minority carrier contributions (i.e. bipolar effects). We would expect some contribution of bipolar effects at high temperatures for the small bandgap HHs, but in this work we only consider room temperature operation, T = 300 K. Furthermore, the effect of spin–orbit coupling on the band structure of the half-Heuslers we consider is negligible near the band edges, which determine transport, as reported in other studies by us and others.^53,69 Thus, it is not included in the calculations. Transport calculations were performed using the Boltzmann transport equation within the relaxation time approximation, which is common practice, providing adequate accuracy compared to methods beyond the relaxation time approximation. Our overall PF values are similar to what were experimentally obtained in the literature; thus, we believe that the approximations we employ do not affect our results quantitatively in a noticeable manner.

5. Summary and conclusion

In summary, we have studied the thermoelectric power factor (PF) properties of a group of 13 half-Heusler (HH) alloys, with the goal of identifying the scattering mechanisms that mainly determine transport and provide as computationally accurate predictions as possible for the PF. DFT-based ab initio simulations were carried out to calculate the electronic structures and static value of the dielectric constants, whereas the deformation potentials used were obtained from density functional perturbation theory (DFPT) with Wannierization. Boltzmann transport simulations were performed using the code ElecTra. We find that the PFs of HHs at 300 K reside mostly between 5 and 10 mW/mK2 for both n-type and p-type materials, with the higher performers being the ones with high band degeneracies (p-type) and/or weaker polar nature (n-type). This could provide a roughly generic range for the power factor of HHs. We find that the positioning of the valleys in the Brillouin zone is of less importance, since inter-valley scattering is effectively suppressed to a large degree across materials. We compared in detail the strength and influence that different scattering mechanisms have on the thermoelectric properties. We found that the Coulombic scattering mechanisms, ionized impurity scattering (IIS) and polar optical phonon (POP) scattering combined determine transport more significantly compared to their non-polar counterparts, namely, ADP and ODP, with IIS being the stronger of all mechanisms. In fact, we showed that the POP and IIS mechanisms combined determine the thermoelectric PF on average by around 65%. We believe that the findings of this work are not only applicable for the group of HHs we have studied, but they can also be applied to the understanding, design and optimization studies for other HHs and polar materials in general. For computational purposes, knowledge of the importance of POP and IIS compared to ADP and ODP (which are computationally much more expensive ab initio) would suggest that transport examination of such materials can be performed at first order based on POP and IIS, which are computationally much cheaper.

Conflicts of interest

There are no conflicts to declare.

Data availability

The code for ElecTra can be found at https://zenodo.org/records/6951369 with https://doi.org/10.5281/zenodo.6951369. The version of the code employed for this study is version ElecTra v1.4.3.2.

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information: complementary simulation results regarding: (i) extraction of matrix elements, (ii) effects of scattering times, (iii) effects of screening and screening lengths, (iv) clarifying the power factor peak values for each material separately, (v) the behavior of the Seebeck coefficient for all materials under consideration for different scattering cases, (vi) the bandstructures of the materials considered, (vii) the density of states (DOS) and values for conductivity and DOS effective masses for the materials considered, and finally (viii) a grouping of the averaged values of the power factors in a way to understand the bandstructure features of the p-type materials that provide the highest power factors. See DOI: https://doi.org/10.1039/d5ta02998h.

Acknowledgements

This work received funding from the UK Research and Innovation fund (project reference EP/X02346X/1).

References

1. D. Beretta, N. Neophytou, J. M. Hodges, M. G. Kanatzidis, D. Narducci, M. Martin- Gonzalez, M. Beekman, B. Balke, G. Cerretti, W. Tremel, A. Zevalkink, A. I. Hofmann, C. Müller, B. Dörling, M. Campoy-Quiles and M. Caironi, Mater. Sci. Eng., R, 2019, 138, 100501.
2. C. Artini, et al., Nanotechnology, 2023, 34, 292001.
3. G. J. Snyder and E. S. Toberer, Nat. Mater., 2008, 7, 105.
4. J.-W. G. Bos, T. Mohanty, T. D. Sparks, W. Xie, A. Weidenkaff, S. Grasso, R. Zhang, M. J. Reece, T. Wang, J. S. Son, S. Akbar, I. S. Nandhakumar, R. Tuley, C. Koz, R. He, P. Ying, A. Bahrami, V. Pacheco, K. Nielsch, R. Grau-Crespo, L. M. Antunes, K. T. Butler, N. Neophytou, R. Dutt, B. Sahni, N. S. Chauhan, T. Mori, M. Parzer, F. Garmroudi, A. Riss, E. Bauer, C. Zeng, E. Bilotti, C. You, O. Fenwick, P. Vaqueiro, E. Guilmeau, A. Das, K. Biswas, Y. Liu, C. Fu, T. Zhu, G. Rogl, P. F. Rogl, P. Mangelis, T. Kyratsi and R. Funahashi, J. Phys. Energy, 2025 DOI:10.1088/2515-7655/ae2d98.
5. S. Chen, K. C. Lukas, W. Liu, C. P. Opeil, G. Chen and Z. Ren, Adv. Energy Mater., 2013, 3, 1210.
6. X. Yan, W. Liu, S. Chen, H. Wang, Q. Zhang, G. Chen and Z. Ren, Adv. Energy Mater., 2013, 3, 1195.
7. T. Ghosh, M. Dutta, D. Sarkar and K. Biswas, J. Am. Chem. Soc., 2022, 144, 10099.
8. A. Bhui, P. V. Matukumilli, S. Biswas, A. Ahad, A. Ghata, D. Rawat, M. Dutta, A. Soni, U. V. Waghmare and K. Biswas, J. Am. Chem. Soc., 2025, 147, 758–3768.
9. K. Biswas, J. He, I. D. Blum, C.-I. Wu, T. P. Hogan, D. N. Seidman, V. P. Dravid and M. G. Kanatzidis, Nature, 2012, 489, 414.
10. S. Kukreti, S. Ramawat and A. Dixit, Phys. Rev. B, 2024, 110, 245202.
11. J. He, M. Amsler, Y. Xia, S. S. Naghavi, V. I. Hegde, S. Hao, S. Goedecker, V. Ozoliņš and C. Wolverton, Phys. Rev. Lett., 2016, 117, 046602.
12. Y. Wang, Y. Zhao, J. Ni and Z. Dai, Adv. Funct. Mater., 2024, 34, 2410983.
13. A. Suwardi, J. Cao, L. Hu, F. Wei, J. Wu, Y. Zhao, S. H. Lim, L. Yang, X. Y. Tan and S. W. Chien, et al., J. Mater. Chem. A, 2020, 8, 18880.
14. S. Ramawat, S. Kukreti and A. Dixit, Phys. Rev. Mater., 2023, 7, 085403.
15. S. Bhattacharya, M. J. Skove, M. Russell, T. M. Tritt, Y. Xia, V. Ponnambalam, S. J. Poon and N. Thadhani, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 77, 184203.
16. Y. Pei, J. Lensch-Falk, E. S. Toberer, D. L. Medlin and G. J. Snyder, Adv. Funct. Mater., 2011, 21, 241.
17. Y. Pei, A. D. LaLonde, N. A. Heinz and G. J. Snyder, Adv. Energy Mater., 2012, 2, 670.
18. Y. Tang, Z. M. Gibbs, L. A. Agapito, G. Li, H.-S. Kim, M. B. Nardelli, S. Curtarolo and G. J. Snyder, Nat. Mater., 2015, 14, 1223.
19. T. Koga, S. B. Cronin and M. S. Dresselhaus, MRS Online Proc. Libr., 2000, 626, Z4.
20. K. Imasato, S. D. Kang, S. Ohno and G. J. Snyder, Mater. Horiz., 2018, 5, 59.
21. A. Li, M. K. Brod, Y. Wang, K. Hu, P. Nan, S. Han, Z. Gao, X. Zhao, B. Ge and C. Fu, et al., Advanced Science, 2023, 10, 2302086.
22. S. Guo, S. Anand, M. K. Brod, Y. Zhang and G. J. Snyder, J. Mater. Chem. A, 2022, 10, 3051.
23. K. Li, L. Sun, W. Bai, N. Ma, C. Zhao, J. Zhao, C. Xiao and Y. Xie, J. Am. Chem. Soc., 2024, 146, 14318.
24. M. Y. Toriyama and G. J. Snyder, Mater. Horiz., 2024, 11, 1188.
25. Z. Li, K. Pal, H. Lee, C. Wolverton and Y. Xia, Nano Lett., 2024, 24, 5816.
26. F. Garmroudi, M. Parzer, A. Riss, C. Bourgès, S. Khmelevskyi, T. Mori, E. Bauer and A. Pustogow, Sci. Adv., 2023, 9, eadj1611.
27. A. Riss, F. Garmroudi, M. Parzer, C. Eisenmenger-Sittner, A. Pustogow, T. Mori and E. Bauer, Phys. Rev. Mater., 2024, 8, 095403.
28. P. Graziosi, K.-I. Mehnert, R. Dutt, J.-W. G. Bos and N. Neophytou, PRX Energy, 2024, 3, 043009.
29. H.-S. Kim, N. A. Heinz, Z. M. Gibbs, Y. Tang, S. D. Kang and G. J. Snyder, Mater. Today, 2017, 20, 452.
30. Y. Feng, J. Li, Y. Li, T. Ding, C. Zhang, L. Hu, F. Liu, W. Ao and C. Zhang, J. Mater. Chem. A, 2020, 8, 11370.
31. Z.-Z. Luo, S. Cai, S. Hao, T. P. Bailey, H. Xie, T. J. Slade, Y. Liu, Y. Luo, Z. Chen and J. Xu, et al., J. Am. Chem. Soc., 2022, 144, 7402.
32. J. Zhu, Q. Ren, C. Chen, C. Wang, M. Shu, M. He, C. Zhang, M. D. Le, S. Torri and C.-W. Wang, et al., Nat. Commun., 2024, 15, 2618.
33. C. Kumarasinghe and N. Neophytou, Phys. Rev. B: Condens. Matter Mater. Phys., 2019, 99, 195202.
34. Y. Pei, X. Shi, A. LaLonde, H. Wang, L. Chen and G. J. Snyder, Nature, 2011, 473, 66.
35. G. K. Madsen and D. J. Singh, Comput. Phys. Commun., 2006, 175, 67.
36. G. K. Madsen, J. Carrete and M. J. Verstraete, Comput. Phys. Commun., 2018, 231, 140.
37. P. Graziosi, C. Kumarasinghe and N. Neophytou, J. Appl. Phys., 2019, 126, 155701.
38. S. Poncé, E. R. Margine, C. Verdi and F. Giustino, Comput. Phys. Commun., 2016, 209, 116.
39. T. Deng, G. Wu, M. B. Sullivan, Z. M. Wong, K. Hippalgaonkar, J.-S. Wang and S.-W. Yang, npj Comput. Mater., 2020, 6, 46.
40. G. Samsonidze and B. Kozinsky, Adv. Energy Mater., 2018, 8, 1800246.
41. P. Graziosi, Z. Li and N. Neophytou, Comput. Phys. Commun., 2023, 287, 108670.
42. A. M. Ganose, J. Park, A. Faghaninia, R. Woods-Robinson, K. A. Persson and A. Jain, Nat. Commun., 2021, 12, 2222.
43. S. Bhattacharya and G. K. Madsen, J. Mater. Chem. C, 2016, 4, 11261.
44. R. Gautier, X. Zhang, L. Hu, L. Yu, Y. Lin, T. O. Sunde, D. Chon, K. R. Poeppelmeier and A. Zunger, Nat. Chem., 2015, 7, 308.
45. X. Jia, Y. Deng, X. Bao, H. Yao, S. Li, Z. Li, C. Chen, X. Wang, J. Mao and F. Cao, et al., npj Comput. Mater., 2022, 8, 34.
46. R. J. Quinn and J.-W. G. Bos, Mater. Adv., 2021, 2, 6246.
47. R. Chen, H. Kang, R. Min, Z. Chen, E. Guo, X. Yang and T. Wang, Int. Mater. Rev., 2024, 69, 83.
48. C. Fu, T. Zhu, Y. Liu, H. Xie and X. Zhao, Energy Environ. Sci., 2015, 8, 216.
49. C. Fu, T. Zhu, Y. Pei, H. Xie, H. Wang, G. J. Snyder, Y. Liu, Y. Liu and X. Zhao, Adv. Energy Mater., 2014, 4, 1400600.
50. R. Dutt, J. Bhattacharya and A. Chakrabarti, J. Phys. Chem. Solids, 2022, 167, 110704.
51. S. Han, S. Dai, J. Ma, Q. Ren, C. Hu, Z. Gao, M. Duc Le, D. Sheptyakov, P. Miao and S. Torii, et al., Nat. Phys., 2023, 19, 1649.
52. Y. Liu, H. Xie, C. Fu, G. J. Snyder, X. Zhao and T. Zhu, J. Mater. Chem. A, 2015, 3, 22716.
53. C. Kumarasinghe and N. Neophytou, Phys. Rev. B, 2019, 99, 195202.
54. Q. Xiong, G. Han, G. Wang, X. Lu and X. Zhou, Adv. Funct. Mater., 2024, 34, 2411304.
55. T. Zhu, C. Fu, H. Xie, Y. Liu and X. Zhao, Adv. Energy Mater., 2015, 5, 1500588.
56. W. Feng, D. Xiao, Y. Zhang and Y. Yao, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 235121.
57. D. I. Bilc, G. Hautier, D. Waroquiers, G.-M. Rignanese and P. Ghosez, Phs. Rev. Lett., 2015, 114, 136601.
58. R. J. Quinn, Y. Go, A. B. Naden, A. Bojtor, G. Paráda, A. K. Shawon, K. Domosud, K. Refson, A. Zevalkink and N. Neophytou, et al., Adv. Phys. Res., 2025, 2400179.
59. M. Fischetti, Phys. Rev. B: Condens. Matter Mater. Phys., 1991, 44, 5527.
60. M. Lundstrom, Meas. Sci. Technol., 2002, 13, 230.
61. L. Pan, S. Mitra, L.-D. Zhao, Y. Shen, Y. Wang, C. Felser and D. Berardan, Adv. Funct. Mater., 2016, 26, 5149.
62. J. Mao, J. Shuai, S. Song, Y. Wu, R. Dally, J. Zhou, Z. Liu, J. Sun, Q. Zhang and C. Dela Cruz, et al., Proc. Natl. Acad. Sci. U. S. A., 2017, 114, 10548.
63. F. Garmroudi, M. Parzer, A. Riss, A. Pustogow, T. Mori and E. Bauer, Phys. Rev. B, 2023, 107, L081108.
64. P. Graziosi, C. Kumarasinghe and N. Neophytou, ACS Appl. Energy Mater., 2020, 3, 5913.
65. Z. Li, P. Graziosi and N. Neophytou, npj Comput. Mater., 2024, 10, 8.
66. I. Welland and D. K. Ferry, Semicond. Sci. Technol., 2019, 34, 064003.
67. S. E. A. Akhtar and N. Neophytou, ACS Appl. Energy Mater., 2025, 8, 1609–1619.
68. Y. Go, R. Dutt and N. Neophytou, Phys. Rev. B, 2025, 111, 195211.
69. J. Zhou, H. Zhu, T.-H. Liu, Q. Song, R. He, J. Mao, Z. Liu, W. Ren, B. Liao and D. J. Singh, et al., Nat. Commun., 2018, 9, 1721.
70. H. Zhu, R. He, J. Mao, Q. Zhu, C. Li, J. Sun, W. Ren, Y. Wang, Z. Liu and Z. Tang, et al., Nat. Commun., 2018, 9, 2497.
71. K. Ciesielski, I. Wolańska, K. Synoradzki, D. Szymański and D. Kaczorowski, Phys. Rev. Appl., 2021, 15, 044047.
72. A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner and G. Ceder, et al., APL Mater., 2013, 1, 011002.
73. Z. Li, P. Graziosi and N. Neophytou, Phys. Rev. B: Condens. Matter Mater. Phys., 2021, 104, 195201.
74. Z. Li, P. Graziosi and N. Neophytou, npj Comput. Mater., 2024, 10, 8.
75. B. Sahni, Y. Zhao, Z. Li, R. Dutt, P. Graziosi and N. Neophytou, Mater. Horiz., 2025, 12, 10255.
76. H. Fröhlich, Adv. Phys., 1954, 3, 325.
77. C. Herring and H. Brooks, Phys. Rev. B: Condens. Matter Mater. Phys., 1951, 83, 879–887.
78. X. Ai, W. Xue, L. Giebeler, N. Pérez, B. Lei, Y. Zhang, Q. Zhang, K. Nielsch, Y. Wang and R. He, Adv. Energy Mater., 2024, 14, 2401345.
79. D. Zou, S. Xie, Y. Liu, J. Lin and J. Li, J. Appl. Phys., 2013, 113, 193705.
80. R. He, D. Kraemer, J. Mao, L. Zeng, Q. Jie, Y. Lan, C. Li, J. Shuai, H. S. Kim and Y. Liu, et al., Proc. Natl. Acad. Sci. U. S. A., 2016, 113, 13576.
81. C. Fu, S. Bai, Y. Liu, Y. Tang, L. Chen, X. Zhao and T. Zhu, Nat. Commun., 2015, 6, 8144.
82. H. Zhu, J. Mao, Y. Li, J. Sun, Y. Wang, Q. Zhu, G. Li, Q. Song, J. Zhou and Y. Fu, et al., Nat. Commun., 2019, 10, 270.
83. P. Graziosi and N. Neophytou, arXiv, 2019, preprint, arXiv:1912.10924v3,  DOI:10.48550/arXiv.1912.10924.

---