Access the full text.
Sign up today, get DeepDyve free for 14 days.
Material Descriptors for the Discovery of Efficient Thermoelectrics †‡ † † Patrizio Graziosi,* Chathurangi Kumarasinghe, Neophytos Neophytou School of Engineering, University of Warwick, Coventry, CV4 7AL, UK Consiglio Nazionale delle Ricerche – Istituto per lo Studio dei Materiali Nanostrutturati, CNR – ISMN, via Gobetti 101, 40129, Bologna, Italy KEYWORDS. thermoelectrics, power factor, materials design, charge transport, half-Heusler 1 ABSTRACT. The predictive performance screening of novel compounds can significantly promote the discovery of efficient, cheap, and non-toxic thermoelectric materials. Large efforts to implement machine-learning techniques coupled to materials databases are currently being undertaken, but the adopted computational methods can dramatically affect the outcome. With regards to electronic transport and power factor calculations, the most widely adopted and computationally efficient method, is the constant relaxation time approximation (CRT). This work goes beyond the CRT and adopts the proper, full energy and momentum dependencies of electron-phonon and ionized impurity scattering, to compute the electronic transport and perform power factor optimization for a group of half-Heusler alloys. Then the material parameters that determine the optimal power factor based on this more advanced treatment are identified. This enables the development of a set of significantly improved descriptors that can be used in materials screening studies, and which offer deeper insights into the underlying nature of high performance thermoelectric materials. We 𝑛 𝘀 𝑣 𝑟 have identified as the most useful and generic descriptor, a combination of the 𝐷 𝑚 𝑜 cond number of valleys, the dielectric constant, the conductivity effective mass, and the deformation potential for the dominant electron-phonon process. The proposed descriptors can accelerate the discovery of new efficient and environment friendly thermoelectric materials in a much more accurate and reliable manner, and some predictions for very high performance materials are presented. 2 1. Introduction Thermoelectric generators (TEG) convert heat flow into useful electrical power and can provide energy savings, thus reduce our dependence on fossil fuels. However, the use of toxic materials without a sufficient energy conversion efficiency is currently hampering the large-scale exploitation 1, 2, 3, 4 of the thermoelectric (TE) technology. Efficient predictive simulation of materials’ performance from first principles, by screening the myriad of possible compounds and alloys, can 1, 2, 5-10 strongly boost the discovery of novel TE materials. The accuracy and predicting ability of such efforts, however, strongly depend on the appropriateness of the screening parameters and descriptors for complex bandstructure materials, whereas the currently used descriptors are developed based on simplified parabolic bands and elastic phonon-limited scattering. The proper identification of the parameters that govern the materials’ performance, will accelerate the discovery process and drive the experimental research reliably. 𝜎 𝑆 The dimensionless figure of merit = 𝑇 quantifies the TE material conversion 𝑘 +𝑘 𝐿 𝑒 efficiency, where is the electrical conductivity, S the Seebeck coefficient, kL and ke are respectively the lattice and electronic thermal conductivities, and T is the absolute temperature. 12, Materials research to date revolved around nanostructuring to reduce the thermal conductivity. 13 2 Recently, several strategies to improve the power factor PF = S have been identified, a direction 14 - 17 which attracts significant interest as well. Indeed, an analysis of the experimental literature on 15, 18 more than 11,000 compounds. reveals that the materials with the highest ZT are generally found more often amongst the materials with high power factors. Several material families are currently under intense investigation, such as clathrates, 1, 17 - 19 skutterudites, half-Heuslers, oxides, selenides, silicides, and others. Most of these materials have complex bandstructures, with multiple valleys and bands that extend in the entire Brillouin 𝑍𝑇 Zone (BZ). In these materials, the electronic scattering processes also have complex energy, momentum and band dependencies. The most common computational approach is to extract the electronic bandstructures using ab initio methods, and afterwards use the Boltzmann Transport 1, 20 Equation (BTE) to extract the TE coefficients. As the extraction of the scattering rates used within the BTE is difficult, the constant relaxation time (CRT) or constant mean free path (CMFP) approximations are mostly used, where a CRT around 10 fs, or a CMFP of around 5 nm, are 5, 6, 10, 21, 22 routinely employed. These approximations have the advantage of being computationally efficient, but have the disadvantage of providing uncertain and rather arbitrary outcomes, both quantitatively and qualitatively, with respect to materials ranking, temperature trends, and optimal 1, 23 carrier density. Furthermore, by ignoring all energy, momentum and band dependence of the scattering processes, we have previously showed that the CRT approximation smears out the richness of the bandstructure that these complex materials possess. The importance of energy dependent scattering process treatment is gaining ground recently, as it is shown that accurately treating the energy dependencies can lead to much better, even opposite conclusions to what is predicted by simpler methods. Computational methods vary in the treatment of electron-phonon scattering, depending on the required accuracy and computational complexities. In reference (24), by considering the energy dependence of acoustic phonon scattering, the TE coefficients were computed in the effective mass approximation and a combination of parameters to identify optimized materials was suggested. Reference (25) considered elastic scattering by acoustic phonons in a full band approach and identify the important role of the energy dependences, and the role of the (acoustic) deformation potential as a performance parameter. In references (2) and (20), the authors employed a full band approach and a numerical scheme for electron-phonon scattering based on calculating the electron-phonon coupling matrix elements. Due to the computational burden, however, it was proposed that a 4 constant optical phonon frequency is used in the scattering description to sample the Brillouin zone. Other than for 2D materials, it is computationally challenging to utilize the entire phonon and electron dispersions for full transport calculations, this, it is common to employ reasonable approximations. In this work, we perform an investigation for the thermoelectric power factor of a group of 14 half-Heusler compounds using DFT bandstructures and the full energy, momentum, and band dependence of electron-phonon scattering and ionized impurity scattering. We employ an advanced simulator with such capabilities, described in Section 2. We then compare, in section 3, the appropriateness of the currently used descriptors and develop a new set of descriptors that behave in a better way in predicting materials ranking in terms of the power factor, especially at the optimal doping level where the role of the ionized impurities is dominant. Finally, we discuss the impact of these descriptors on the understanding of the material properties that lead to optimal performance. 2. Methods We use a novel simulator that considers all appropriate scattering mechanisms (acoustic phonon, optical phonon, ionized impurity scattering), and the full energy, momentum, and band dependence of the relaxation times. The simulation approach consists of three stages: 1) calculation of the electronic structure using Density Functional Theory (DFT), 2) numerical extraction of the scattering rates, 3) calculation of the transport coefficients. The electronic bands are calculated 26, 27 using the Quantum Espresso package. The projector augmented wave technique is used with the PBE-GGA functional and the spin orbit coupling (SOC) has been considered as well. A kinetic energy cut-off greater than 60 Ry is used for the wave functions and an energy convergence criterion -8 of 10 Ry for self-consistency was adopted. The computed bandgaps are very close to the 5 experimental ones where available. For instance, for the NbFeSb compound we obtain a bandgap of 0.53 eV, in excellent agreement with the experimental value of 0.51 eV. The 3D bandstructure was calculated using a 51x51x51 Monkhorst–Pack k-point mesh on the primitive unit cell of the reciprocal lattice. The coordinates of the k-points, originally described in the reciprocal unit cell coordinate system, are transformed in the orthogonal Cartesian coordinate system. We then calculate the transport coefficient tensors in the Cartesian coordinates x,y,z using 16, 20, 29 - 31 the linearized Boltzmann Transport Equation (BTE). The electrical conductivity and the Seebeck coefficient S tensors are expressed as 𝜕 𝑓 𝜎 = 𝑞 𝛯 (𝐸 )(− ) 𝑑𝐸 (1a) 𝑖𝑗 0 𝑖𝑗 𝜕𝐸 𝑞 𝑘 𝜕 𝑓 𝐸 −𝐸 0 𝐵 0 𝐹 𝑆 = ∫ 𝛯 (𝐸 )(− ) 𝑑𝐸 (1b) 𝑖𝑗 𝑖𝑗 𝜎 𝜕𝐸 𝑘 𝑇 where EF, T, q0, kB, are the Fermi level, the absolute temperature, the electronic charge, and the Boltzmann constant, respectively, f0 is the equilibrium Fermi distribution, and ij are the Cartesian ( ) directions (in all the calculations we show, i = j = x). 𝛯 𝐸 is the transport distribution function (TDF) defined as: 𝑑𝐴 2 1 𝒌 ,𝑛 ,𝐸 ( ) ∑ 𝛯 𝐸 = 𝑣 𝑣 𝜏 𝛿 = ∯ 𝑣 𝑣 𝜏 (2) 𝑖𝑗 𝒌 ,𝑛 𝑖 (𝒌 ,𝑛 ) 𝑗 (𝒌 ,𝑛 ) 𝑖 (𝒌 ,𝑛 ) (𝐸 −𝐸 ) 𝑖 (𝒌 ,𝑛 ,𝐸 ) 𝑗 (𝒌 ,𝑛 ,𝐸 ) 𝑖 (𝒌 ,𝑛 ,𝐸 ) 𝒌 ,𝑛 (2𝜋 ) ℏ |𝑣 | 𝐸 (𝒌 ,𝑛 ,𝐸 ) where the sum over all the states in the Brillouin Zone (BZ) becomes a surface integral over the constant energy surface 𝔏 of energy E defined by the delta function for the band of index n. In equation (2) 𝑣 is the i-component of the band velocity of the transport state defined by the 𝑖 (𝒌 ,𝑛 ,𝐸 ) wave vector k in the band n at energy E, 𝜏 is its momentum relaxation time (combining the 𝑖 (𝒌 ,𝑛 ,𝐸 ) relaxation times of each scattering mechanism, defined in equation (3) below, using Matthiessen’s rule). 𝑑𝐴 is the surface area element associated to the (k,n,E) transport state, 𝑣 is its band 𝒌 ,𝑛 ,𝐸 (𝒌 ,𝑛 ,𝐸 ) 𝐵𝑍 𝑑𝐴 𝒌 ,𝑛 ,𝐸 velocity and its DOS. For each transport state and each scattering mechanism ms, the ℏ |𝑣 | (𝒌 ,𝑛 ,𝐸 ) (𝑚 ) 𝑛 relaxation time 𝜏 is given by a surface integral on the final constant energy surface 𝔏 𝑖 (𝒌 ,𝑛 ,𝐸 ) 𝑑𝐴 1 1 ′ 𝑖 (𝒌 ) (𝑚 ) 𝒌 = ′ |S ′ | (1 − ) (3) (𝑚 ) 3 𝒌 ,𝒌 𝑠 𝔏 𝜏 (2𝜋 ) ℏ|𝒗 | 𝑣 ′ ′ 𝑖 (𝒌 ,𝑛 ,𝐸 ) 𝑖 (𝒌 ,𝑛 ,𝐸 ) 𝐸 (𝒌 ) where k’ represents a possible final state (the symbol k’ lumps its band index 𝑛 , that may be ′ ′ ′ different from the initial one, and its energy 𝐸 , with 𝐸 = 𝐸 for elastic scattering and 𝐸 = E ± ℏ𝜔 in the case of inelastic scattering with phonon of frequency 𝜔 ). 𝑖 (𝒌 ′) |S | is the transition rate between k and k’, and 𝑣 is the carrier velocity. The (1 − ) k,k’ 𝑖 (𝒌 ,𝑛 ,𝐸 ) 20, 23, 29 - 31 term is an approximation for the momentum relaxation. |S | is derived from Fermi’s k,k’ Golden Rule for the different scattering mechanisms, and the relevant scattering parameters are 9, 32, 33 taken from first-principle calculations. Specifically, we employ deformation potential theory for electron-phonon scattering (acoustic and optical), and ionized impurity scattering with charge screening. The corresponding transition rates for acoustic phonons, optical phonons and ionized impurity scattering are: 𝑘 𝑇 (𝐴𝐷𝑃 ) |S | = (4a) 𝒌 ,𝒌 ℏ𝜌 𝑣 𝜋 𝐷 1 1 (𝑂𝐷𝑃 ) |S | = (𝑁 + ∓ ) (4b) 𝜔 , 𝒌 ,𝒌 2 2 𝑍 𝑞 𝑁 ( ) 2𝜋 𝑆𝐼𝐼 |S | = (4c) 2 2 2 𝒌 ,𝒌 ℏ 𝘀 𝘀 𝟐 0 1 (|𝒌 −𝒌 | + ) Above, D is the acoustic deformation potential, is the mass density, v is the sound velocity ADP s 1 2 defined as 𝑣 = 𝑠 + 𝑠 , which is important when the phonon bands are not isotropic, where 𝑠 𝑙 𝑡 3 3 𝐾 + ⁄ 𝐺 𝑉 𝑉 3 𝑉 𝑠 = and 𝑠 = √ are the longitudinal and transverse sound speeds and KV and GV are 𝑙 𝑡 𝜌 𝜌 𝑖𝑚𝑝 𝜌𝜔 𝐵𝐸 𝑂𝐷𝑃 𝐴𝐷𝑃 𝜋𝐷 the bulk and shear moduli. D is the optical deformation potential, is the longitudinal optical ODP phonon frequency in the single mode approximation with constant frequency over the entire reciprocal lattice unit cell. N is its population density given by the Bose-Einstein statistics where the ‘+’ and ‘–’ signs indicate the absorption and emission processes, respectively. Nimp is the impurity density, Z the impurity charge, and and are the vacuum and the static relative 0 r 𝘀 𝘀 𝜕 𝐸 𝑟 0 𝐹 permittivities. LD the Debye screening length in 3D defined as 𝐿 = , where 𝓃 is the carrier 𝑞 𝜕 𝓃 density and 𝜕 𝓃 /𝜕 𝐸 is the variation of the carrier density with the Fermi level, which is temperature and doping dependent. The explicit use of 𝜕 𝓃 /𝜕 𝐸 enables us to apply the equation also under degenerate doping conditions. For the longitudinal optical phonon energy we adopt an averaged 2, 23 value that is assumed to be constant over the reciprocal lattice unit cell. A representative example of the compounds studied is TiCoSb, whose bandstructure is shown in Figure 1a, and its conventional Zinc Blende unit cell in Figure 1b. Multiple bands of different curvatures in different directions compose the bandstructure and participate in transport. The shaded areas in Figure 1a highlight the energy range considered in our calculations, which extends up to 0.7 eV beyond the band edge, and corresponds to about 7kBT beyond the highest Fermi level considered in the calculations at the highest temperature considered, T = 900 K. 7k T 𝜕 𝑓 is largely sufficient as the largest part of the term in equation (1) extends for about 5kBT above 𝜕𝐸 the Fermi level, and represents an optimal tradeoff between calculation accuracy and computational cost. Examples of the carrier scattering transitions we consider in the calculations are sketched in Figure 1c, where two simple parabolic bands are displayed for illustrative purposes. The selection rules and details of the strength of the electron-phonon coupling of each initial state to all other states individually are not yet well established for half-Heusler alloys in general. Thus, in this work we consider both intra-valley and inter-valley scattering events for phonon 8 scattering, and allow transitions from an initial state to all other energy conserving states of the BZ. The Acoustic Deformation Potential (ADP) scattering that we consider is elastic. The scattering with optical phonons of energy ℏ𝜔 , considered within the Optical Deformation Potential (ODP) theory, is inelastic. Phonon scattering is considered both intra and inter-band. The scattering due to ionized dopant impurities, (IIS), is elastic and considered as intra-band only, following the common 30, 35, 36 treatment in the literature. We adopt the relative dielectric constants from reference (9) where they were calculated from first principles. The full list of first principles computed input parameters is listed in the Supporting Information. Figure 1. (a) Bandstructure of TiCoSb, as a representative example. The shaded regions highlight the energy windows used in the calculations. (b) Conventional cubic zincblende unit cell of TiCoSb as a representative of all 14 compounds investigated. (c) Types of scattering transition we consider depicted using two parabolic bands. The inelastic processes are due to the absorption or emission of optical phonons of energy ℏ. (d) Reciprocal space constant energy surfaces at E = -0.125 eV for two bands (blue and orange) below the valence band edge of TiCoSb. Warped shapes with narrow wings and tubular connections between minima can be observed. (e) Experimental (filled symbols) and calculated (open symbols) power factor versus carrier concentration. Symbols with 9 blue edge are for n-type materials and with red edge are for p-type materials. The filling color is an indication of the temperature: green for T < 600 K and yellow for T ≥ 600 K. Each material is 37 38 39, 40 indicated by the symbols: HfCoSb: square, HfNiSn: circle, NbCoSn: up-pointing triangle, 28 41 37 NbFeSb: down-pointing triangle, ScNiSb: diamond, TiCoSb: left-pointing triangle, TiNiSn: 42, 43 44 41 45 right-pointing triangle, YNiBi: hexagon, YNiSb: star, ZrCoBi: pentagon, ZrCoSb: dot 37, 45, 46 47, 48 centered square, ZrNiSn: dot centered diamond. The inset compares the computed and measured PFs at similar conditions related to the carrier density and temperature conditions. The dashed line has unity slope. The TDF defined in equation (2) is a surface integral on the relevant constant energy surfaces. An example of such surfaces is displayed in Figure 1d, where two constant energy surfaces from two valence bands of TiCoSb at 0.125 eV below the band edge are displayed (blue and orange). These surfaces are warped, with elongated regions, necks and pockets, very different from the regular ellipsoids that enable analytical expressions for the relaxation times. They extend over nearly the whole reciprocal unit cell, all of which needs then to be included numerically in all scattering events. A full band numerical treatment becomes compulsory, which is done by the extraction of the k-points at constant energies, essentially converting the usual E(k) into a k(E). The transport calculations presented here can require from ~ 20 up to ~ 600 computation CPU hours for each material, depending on the bandstructure complexity, executed in parallel on ~ 20 CPUs. 3. Results and discussion 10 This section presents the main results, starting from the transport properties of the half-Heuslers considered, to the developed material descriptors for transport under: i) phonon-limited and then under phonon plus IIS transport, ii) at 300 K first and then at higher temperatures, and iii) under unipolar and bipolar transport. An initial glance at the results is shown in Figure 1e, that presents measured PF values versus carrier density, and compares them with the computed values at the same (or nearby) carrier density 28, 37 - 48 and temperature. The symbols, one for each material, denote: filled/empty for experimental/calculated, blue/red boundary for n-/p-type, and green/yellow filling for measured temperature below/above 600 K. The inset shows the correlation between computed and experimental values, with the dashed line having unity slope. In some cases, the calculation shows excellent agreement with experimental values, but in most cases the calculated PFs are higher than the measured ones by a factor of two or three. Considering that experimental samples are polycrystalline while our calculations assume perfect single crystals, the similarity in the PF values indicates the credibility of our computational method. An interesting observation is the case of 2 21 -3 p-type NbFeSb (down-pointing triangles at PF ~ 10 mW/mK and densities ~ 10 cm ). The computed PF values are very close to the measured ones, and those samples have the highest degree of crystallinity, with the largest grain sizes. 3.1 Thermoelectric charge transport properties 11 Both electron transport for n-type materials and hole transport for p-type materials are considered. We compute the conductivity , the Seebeck coefficient S, and the power factor PF = S , under three scattering scenarios: i) phonon-limited transport τph(E), ii) scattering of charge carriers with both phonons and ionized impurities τ (E), and iii) CRT with τ = 10 fs, as commonly adopted ph,IIS c 1, 5, 7, 10, 21 in the literature. The conduction and valence bands are analyzed separately assuming unipolar transport, but also combined later on, in section 3.3, to describe bipolar transport. We plot all quantities in terms of the reduced Femi level F, which denotes the distance of the Fermi level from the band edge, and has a one-to-one correspondence with the carrier density. Figure 2. n-type thermoelectric coefficients versus the relative position of the Fermi level F for the 14 half-Heusler compounds under investigation for unipolar transport. Row-wise, simulation results for three different scattering scenarios are presented: (a-c) electron phonon-limited scattering ph (E); (d-f) phonon plus ionized impurity scattering, ph,IIS (E); (g-i) constant relaxation 12 time (CRT) approximation with = 10 fs. Column-wise: electrical conductivity σ (a), (d) and (g), Seebeck coefficient S (b), (e), (h), and power factor PF (c), (f), (i). The data refer to room temperature conditions. The inset in (e) compares the computed Seebeck coefficients for TiNiSn (dark yellow) with the corresponding experimental values (magenta dots) measured at room temperature in quasi-single crystals. 3.1.1 n-type thermoelectric coefficients (unipolar transport) The dependence of the charge transport coefficients (, S and PF) of the considered half-Heuslers upon F at room temperature (T = 300 K) is shown in Figure 2. F = 0 means that the Fermi level is at the band onset, < 0 in the band gap, and > 0 is for heavily doped, degenerate conditions F F with the Fermi level into the conduction band. Each compound is depicted by a specific color as displayed in the legend in Figure 2a. The transport coefficients , S and PF are displayed column- wise, while the different rows correspond to the different scattering scenarios: scattering with phonons only (E), the adding ionized impurity scattering (IIS) (E), and the CRT ph ph,IIS approximation c. The electrical conductivity, in Figures 2a, 2d, 2g (first column), displays a large variation between the materials, while the S values, in Figures 2b, 2e, 2h (second column), are similar for all the materials for a given scattering scenario. Thus, the electrical conductivity drives the materials ranking, and it is also the quantity that is most sensitive to the specific details used to compute the relaxation times. We note that in the CRT approximation the Seebeck coefficient S does not depend on the value of the relaxation time, and it is in general considered weakly sensitive to the scattering details as well. However, when we consider the energy/wave-vector dependence of the relaxation time, we 13 observe somewhat quantitatively different results compared to the CRT approximation. The average Seebeck coefficient between materials at F = 0 in the phonon-limited case is ~ 250 V/K, in the case of phonon plus IIS ranges between 200 and 300 V/K, but in the CRT case it is ~ 300 V/K. Under IIS consideration, a larger variation between the Seebeck coefficients of the materials is observed (in the range of ~ 50 % around = 0), whereas in the (E) and CRT cases the variation F ph is only ~ 10 %. Thus, a proper scattering treatment, including both, energy dependencies and ionized impurity scattering, should be conducted also in the evaluation of the Seebeck coefficient, rather than only in the case of the conductivity. The inset of Figure 2e compares the S values computed for TiNiSn with experimental values in quasi-single crystals (magenta points). We use values of F which map the experimental carrier density. The excellent agreement supports our computational approach. The ranking of materials with respect to the optimal PF (third column, Figures 2c, 2f and 2i), varies significantly under the different scattering scenarios, driven by the variations in . The optimal Fermi level is at the band edge in the case of phonon-limited scattering (F = 0 in Figure 2c). When IIS is also considered in Figure 2f, the ideal Fermi level position moves from few tens of meV in the bandgap (negative F, extreme case for NbFeSb with optimal F around -40 meV) to few meV into the band (positive , extreme case of NbCoSn with optimal around 10 meV). On F F the contrary, the CRT approximation suggests degenerate conditions as being the optimal ones for all the materials, suggesting even up to one order of magnitude higher doping, which is also more difficult to achieve experimentally. The best candidates for n-type TE materials are NbFeSb (purple line) for its low deformation potentials (D = 1.0 eV, D = 1.6×10 eV/m), low conductivity effective mass m (0.33 m ), ADP ODP cond 0 and high dielectric constant (r = 23.0). We refer to mcond as the effective mass of an isotropic 14 parabolic band that gives the same carrier velocity as the overall material’s full bandstructure, which is a quantity experimentally accessible; see Supporting Information for more details on how we extract it. HfCoSb and HfNiSn also have quite low deformation potentials (DADP = 0.2 eV, DODP = 10 10 1.4×10 eV/m for HfCoSb and D = 0.1 eV, D = 1.8×10 eV/m for HfNiSn), but relatively ADP ODP lower r (especially HfCoSb, r = 17.5). The ranking of these materials slightly changes when IIS is added, possibly because the larger increases the charge screening and improves the conductivity (as we show below). Vice versa, the high εr (26.6, 20.8, 29.0) could favor YNiBi, HfNiSn, ScNiBi, respectively, which end up in second place after NbFeSb. It is worth mentioning that from those only HfNiSn is a known n-type half-Heusler with Sb as the donor in the lattice place 51,52 of Sn while the other compounds are usually p-type. Notice that the advantage of NbFeSb is lost in the case of the CRT approximation, which predicts that NbFeSb is the worst performer. 3.1.2 p-type thermoelectric coefficients (unipolar transport) The same analysis is now applied to hole transport for the same compounds. The dependence of the charge transport coefficients (, S and PF) for different Fermi level positions at 300 K is displayed in Figure 3, as above in Figure 2. Column-wise we show , S and PF, and row-wise the three different scattering scenarios. In the case of hole transport, > 0 describes a Fermi level residing in the bandgap and F < 0 indicates heavily doped conditions with the Fermi level into the valence band. Each compound is depicted by the same color as in Figure 2. As in the case of n-type transport, the electrical conductivity in the first column, Figures 3a, 3e, 3g, is the transport coefficient that displays the larger variation between materials and drives the materials ranking. Again, the optimal doping appears when the Fermi level is placed around the band edge when energy dependent scattering physics is considered ( = 0 in Figure 3c, and 3f). The only exceptions 15 are ScNiBi and ScNiSb, for which the PFs peak at a Fermi level slightly in the gap. Under the CRT approximation, in Figure 3i, the PFs peak at heavily doped degenerate conditions for all the materials, with Fermi levels positioned from 0.05 eV to 0.1 eV into the band, suggesting an optimal carrier density around one order of magnitude higher. The best candidates for p-type TE materials are ScNiBi for its low deformation potentials (DADP = 0.2 eV, D = 1.4×10 eV/m), low m (0.44 m ), and very high dielectric constant ( = 29.0). ODP cond 0 r TiNiSn is a top performer under phonon-limited considerations for its low deformation potentials (D = 0.2 eV, D = 0.9 × 10 eV/m). Fortunately, ScNiBi can be p-doped when substituting ADP ODP 52 39 Bi with Sn and NbCoSn and NbFeSb can be p-doped when replacing Nb with Ti or Zr, or by substitutional doping with Fe in the Co site, like in TiCoSb. When it comes to the CRT approximation, only the NbCoSn (orange line) continues to perform well, but ScNiBi and TiNiSn lose their advantage. The reason NbCoSn still performs well is its simultaneous small conductivity effective mass and high DOS effective mass, as we will highlight later on in the development of the material descriptors. Figure 3. p-type thermoelectric coefficients versus the relative position of the Fermi level F for the 14 half-Heusler compounds under investigation for unipolar transport. Row-wise, simulation results for three different scattering scenarios are presented: (a-c) electron phonon-limited scattering ph (E); (d-f) phonon plus ionized impurity scattering, ph,IIS (E); (g-i) constant relaxation time (CRT) approximation with c = 10 fs. Column-wise: Electrical conductivity σ (a), (d) and (g), Seebeck coefficient S (b), (e), (h), and power factor PF (c), (f), (i). Importantly, however, each scattering scenario can result in different predictions for the PF. In Figure 4, for example, we show the correlation between the PF values computed under the three different scattering assumptions. Blue dots are for electron transport data and red dots for hole transport data. In the legend, the re/h are the Pearson correlation coefficients for electron/hole 17 transport, that range from 1 (complete correlation) to -1 (complete anticorrelation). Figure 4a shows the correlation between the peak PF values computed under the CRT approximation (c) and under phonon-limited scattering treatment, (E). The ranking data are largely uncorrelated, even ph nearly anticorrelated in the n-type case, which reflects back to the very different materials ranking observed in the two scattering cases. A similar absence of correlation between the PF computed values still holds when comparing the CRT results with the phonon plus IIS scattering, ph,IIS(E), shown in Figure 4b. On the other hand, comparing the computed PF under the (E) and the ph ph,IIS(E) cases, in Figure 4c, a sufficient correlation appears for both the n-type and p-type materials (r = 0.82 and r = 0.79). e h 3.2 Descriptors development 3.2.1 Available Descriptors Large efforts are currently being undertaken in order to identify potential high performance TE materials and their alloys using materials screening methods, guided by machine learning 1, 54 18, 55 - 58 techniques. A tremendously large set of materials data is accumulated in databases. Such efforts define appropriate descriptors, i.e. simple parameters or a combination of them, which could signify high TE performance, for example the number of valleys, density of states, effective mass, mobility, etc., and then the material databases are screened to find materials that maximize the descriptor values. In this way, high-throughput screening can be performed, but the confidence in the descriptors needs to be high. 1, 11 A widely used descriptor is the dimensionless quality factor , or “B-factor”, defined 𝑇 𝑘 as 𝛽 = 𝑛 𝜇 , where 𝑛 is the number of valleys, 𝜇 is the mobility of the carriers in a single 𝑣 𝑣 𝑣 𝑣 𝑒 𝜅 valley, 𝜅 is the lattice thermal conductivity, T is the temperature, kB the Boltzmann constant, and 18 ∗3/2 𝜇 𝑚 11 DOS ∗ e the electron charge. A slight variation of this is , where 𝜇 is the mobility and 𝑚 is 0 DOS ∗3/2 the DOS effective mass. The numerators, 𝑛 𝜇 , or 𝜇 𝑚 , are related to the power factor. 𝑣 𝑣 0 DOS Another variation, developed from the mobility dependence on acoustic phonon-limited scattering, 𝑛 𝑐 𝑣 𝑙 is , where 𝑛 is the number of valleys, 𝑚 the density of states (DOS) effective mass for a 𝑣 𝑣 𝑚 𝐷 𝑎𝑐 single valley, 𝐷 the deformation potential for acoustic phonons, 𝑐 = 𝜌 𝑢 with 𝜌 the material 𝑎𝑐 𝑙 𝑠 density and 𝑢 the sound velocity. All these descriptors originate from a semi-empirical description of charge transport under the simple assumptions that the bands are parabolic and the 1, 11 charge transport is phonon-limited. Recent studies that use full numerical bandstructures coupled to the extraction of the electron- phonon scattering rate via the direct calculation of the interaction matrix elements, suggested the ∗ 2 inverse of the density of states effective mass, 1/𝑚 , as a single-parameter descriptor. While DOS this is a valuable outcome, it is still based on phonon-limited scattering conditions. The weak correlation in the materials ranking between each scattering scenario suggests that each case is optimized by different descriptors, which points to the necessity of accurate descriptors to obtain any meaningful results in materials screening. Under the CRT approximation, a descriptor 3⁄ 2 DOS ⁄ that offers satisfactory Pearson coefficients for both carrier types is , as shown in cond Figure 4d. The numerator directly relates to the DOS and the denominator to the squared velocity of the carriers, as it appears in the TDF, equation (2). Here, the m is the density of states effective DOS mass (the effective mass of an isotropic parabolic band that gives the same DOS as the actual bandstructure at the band onset), and mcond is the conductivity effective mass, which is related to the carrier velocity. Details on the extraction of these quantities can be found in the Supporting Information. These parameters are related to the whole bandstructure and automatically include the 19 effect of multiple, and different, valleys. Such descriptor points reasonably towards seeking materials with many and fast carriers, and indicates that the CRT approximation favors materials with higher DOS effective mass. Figure 4. Correlation between the peak PF values computed under different scattering scenarios: constant relaxation time c, electron phonon scattering ph(E), and electron-phonon plus ionized impurity scattering (E). The blue and red colors are for the n-type and p-type materials, ph,IIS respectively. (a) Correlation between CRT and ph(E) transport; (b) correlation between CRT and (E) transport; (c) correlation between (E) and (E) transport. (d) A descriptor proposed ph,IIS ph ph,IIS in this work for the CRT approximation: the ratio between DOS effective mass (m ) raised to the DOS 3/2 power over the conductivity effective mass mcond. The Pearson correlation coefficient r is displayed in the legend, r for electron transport and r for hole transport. e h 20 The situation is less straightforward when we consider energy-dependent electron-phonon scattering, because in addition to having more conducting states, a higher DOS results in a stronger 23, 60 23, 27, 60 scattering rate, and the two effects could cancel each other out. We first examine the adequacy of existing descriptors already proposed in the literature, in Figure 5. We begin with the 𝑛 𝑐 𝑣 𝑙 commonly used , in Figure 5a. 𝑚 represents an ‘average’ DOS mass of a single valley 𝑚 𝐷 𝑣 𝑎𝑐 such that if its density of states is multiplied by the number of valleys, the overall bandstructure DOS is obtained. This is given by 𝑚 = (𝑛 ) 𝑚 that gives a total DOS same as the sum DOS/1𝑣 𝑣 DOS of the DOS of the individual valleys. The unsatisfactory performance of this descriptor with respect to our calculations (r = -0.04 and r = 0.04) resides in two reasons. First, the descriptor e h relies on intra-valley scattering alone, while here we consider inter-valley/inter-band scattering as well. In the case of intra-valley scattering only and phonon-limited transport, the quality factor is likely to provide better correlation. However, in the case of inter-band/inter-valley transitions, the increased scattering into the larger DOS of the final scattering states cancels out the potential 23, 60 benefit, which removes the influence of the number of valleys and effective mass. In general, whether intra- versus inter-band scattering dominates is not known for these materials, but we have no reason to neglect the inter-band/inter-valley, especially when the dominant phonons are the optical ones. Second, acoustic phonon deformation potential in the half-Heuslers we examine is quite low compared to the optical phonon deformation potential which dominates, which is why we use it in the formation of the descriptor. Indeed, the band edge states in the conduction and valence bands of half-Heuslers, are characterized by nearly non-bonding orbitals. This is a situation under which the bonding/anti-bonding interactions of the orbitals that form the band edge states nearly cancel, making the band edges resilient to energy shifts upon structure deformations, thus providing weak electron-acoustic phonon interactions. This resilience is enhanced by the fact 21 that the interactions of those orbitals with orbitals forming other energy states in nearby or higher energies are mostly symmetry forbidden. This does not happen for the electron-optical phonon interactions, which then acquire a significantly higher deformation potential. Indeed, upon substitution of the acoustic deformation potential 𝐷 with the optical phonon one, 𝐷 , in Figure 𝑎𝑐 o 5b, much better agreement is achieved for the n-type case (r = 0.72), but still not for the p-type. The reason is that all n-type materials we consider have the same number of valleys (only X valleys). Thus, whether intra- or inter-valley scattering is considered, it affects all of them similarly. That is not the case for the p-type materials, which have several and varying number of valleys that enter the descriptor (see the Supporting Information for details), but drop out of the TDF in equation (2). Thus, the quality factor descriptor provides adequate correlation for n-type when it includes the dominant deformation potential, regardless it being 𝐷 or 𝐷 . Experimentally, the deformation 𝑎𝑐 o potentials are extracted from the temperature dependence of the mobility, and by assuming a certain effective mass value. Under the effective mass approximation, under either ADP or ODP, the same 34,61 mobility temperature dependence is achieved. Thus, one assumes that the extracted deformation potential will somehow effectively contain both. The single parameter descriptor 1/m , is displayed in Figure 5c and performs well for electron DOS transport (re = 0.72) but not for hole transport (rh = 0.25). This is consistent with reference (2), in which the relaxation time versus the m trend is clearly linear for n-type half-Heusler alloys, but DOS not for the p-type case. Figures 5d, 5e and 5f show the correlation between the peak PF computed at 300 K and the same descriptors as in Figures 5a, 5b and 5c, when we further include the IIS mechanism, ph,IIS(E). This comparison is important, because this is the default experimental situation. The quality factor, Figure 5d, and its adjustment with 𝐷 instead of 𝐷 , Figure 5e, still provides acceptable Pearson o 𝑎𝑐 22 correlation coefficients for electrons (r = 0.60), which somehow justifies its use in rationalizing 45,59 experimental measurements. As for the phonon-limited case, the 1/m descriptor works much DOS better for electrons than for holes (re = 0.76, rh = 0.19), Figure 5f. This is not surprising, as this descriptor is specifically suggested from phonon-limited scattering considerations. Figure 5. Correlation between the peak PF and common descriptors proposed in the literature for the phonon-limited transport, τ (E), top row, and the phonon plus ionized impurity scattering ph transport, τph,IIS(E), bottom row. The Pearson correlation coefficient r is displayed in the legends. (a), (d) A descriptor based on the single band mobility and the number of valleys proposed by Tang et al. (b), (e) The same as in (a) and (d) but with the acoustic deformation potential D , ac substituted by the non-polar optical deformation potential Do that is dominant in these materials. (c), (f) The DOS effective mass m descriptor. DOS 23 3.2.2 Development of new descriptors After, testing the existing descriptors, we present here new descriptors based on the energy and wave-vector dependent scattering mechanisms, as well as IIS, which would be more appropriate for materials screening. The descriptor development is depicted in Figure 6. We first examine the phonon-limited scattering case (top row) and then the phonon plus IIS case (bottom row). In the search for a single parameter descriptor, the ODP deformation potential 𝐷 represents a first intuitive guess because this inelastic scattering mechanism has much higher strength amongst these 10 2 compounds compared to 𝐷 . The correlation between the peak PF values and 1 𝐷 in Figure 𝑎𝑐 o 6a is fairly satisfactory for both n-type and p-type cases (r = 0.54, r = 0.69). The fact that the e h scattering strength is proportional to the square of the deformation potential motivates the choice of squared 𝐷 . This descriptor can be further improved by carrier velocity considerations, by including the bandstructure conductivity effective mass 𝑚 . The 1⁄ 𝐷 𝑚 descriptor in cond o cond Figure 6b offers an astonishing agreement (r = 0.96, r = 0.95). e h We deliberately avoid using parameters related to the DOS due to the reason that under energy dependent scattering the DOS cancels out with the inverse DOS dependence of the relaxation times, 23, 60 and at first order drops out of the transport coefficients. Thus, when looking for effective mass parameters, we chose the conductivity effective mass, which does not cancel out in the evaluation of the TDF in equation (2). Since the deformation potential and conductivity effective mass are two parameters related to the mobility, Figure 6c plots the peak PF values versus the computed lattice mobility (low field and low carrier density), a measurable quantity related to electron-phonon scattering processes. We do not choose the mobility at the carrier density of the PF peak, because the mobility of an experimental specimen at these carrier densities is limited by the IIS. The mobility correlation is very good for n-type materials (r = 0.96) and fairly good for p-type (r = e h 24 1,11 0.78) as also suggested in the past. The mobility can be a valid single parameter descriptor for machine learning studies in databases that collect experimental results, but in the case of computational studies it is not readily available. Figure 6. Correlations between the peak of the PF and various performance descriptors proposed in this work for phonon-limited transport, τph(E), first row, and phonon plus ionized impurity scattering transport, τ (E), second row. (a) A single parameter descriptor based on the optical ph,IIS deformation potential Do. (b) The descriptor in (a) is improved by including the conductivity effective mass mcond. (c) The correlation of the PF peak with the computed lattice low-field, low- density mobility, which is a parameter accessible in experiments. (d-f) Same as the descriptors in (a-c) but applied to the peak PF values computed when also including the IIS mechanism. The computed mobility in (f) is at the density where the peak PF is achieved for each material and carrier type, a parameter accessible in experiments. 25 Adding ionized impurity scattering (IIS) to electron-phonon scattering, (E), comprises a ph,IIS picture closer to realistic TE materials at the doping densities of interest. Testing the descriptors of Figure 6a, 6b and 6c for the peak of the PFs in the phonon plus IIS case, though, indicates reduced correlation. The descriptor in Figure 6d is inadequate for both n-type (re = 0.20) and p-type cases (rh = 0.29). The addition of the conductivity effective mass (1 𝐷 𝑚 ) in Figure 6e largely o cond improves the correlation for both n-type (re = 0.86) and p-type materials (rh = 0.78). In addition, the mobility parameter in Figure 6f (calculated at the density where the peak of the PF appears), which is useful from an experimental point of view, seems applicable for n-type half-Heusler TE materials (r = 0.83) but again not for the p-type materials (r = 0.48). Thus, in the presence of the e h dominant IIS, seems to be a satisfactory descriptor for the PF in the nearly parabolic 𝐷 𝑚 o cond conduction band of the considered half-Heuslers, but it is not adequate for the strongly warped, non-parabolic bandstructures like the valence bands of these materials, in contrast to phonon- limited transport conditions. Thus, a further improvement to the descriptor is necessary. The material dependent term that appears in the IIS scattering rate, is proportional to 1 1 ( ) , with 𝘀 being the relative permittivity, 𝐿 the Debye screening length and q = |k - 2 𝑟 D 2 1 𝜀 𝑞 + 𝑟 2 −1 𝜀 𝜀 𝜕𝓃 𝑟 0 2 23 k’| the scattering vector between initial state k and the final state k’. 𝐿 = ( ) is the 𝑒 𝜕 𝐸 screening length, with 𝘀 the vacuum permittivity, e the electronic charge, E the Fermi level and 𝓃 the carrier density. That term is then proportional to ( ) , in which case the material- 2 2 𝑞 𝐿 +1 2 2 2 dependent parameters appear through the 𝑞 𝐿 ~ 𝑞 𝘀 term, and they become important if it is 𝐷 𝑟 bigger than the unity. The screening lengths are around ~ 1 nm (L at the PF peak ranges from 0.3 nm to 1.9 nm depending on the material). With the cubic lattice constant being around ~ 0.6 nm, 26 2 the 𝑞 𝘀 term becomes important in suppressing scattering and increasing the PF for transitions for which the states are farther than 15% (from 5% to 30% depending on the material) of the Brillouin Zone. This essentially creates conditions for inter-valleys scattering suppression for IIS. We developed a method to estimate an average 〈𝒒 〉 to improve the descriptor. For each possible pair of initial and final states we compute the scattering vector (Figure 7a) and average it for all 〈 〉 the k-points of the constant energy surface E and band index n to get a 𝒒 . Then, we perform 𝐸 ,𝑛 the energy average using the density of states g and the Fermi distribution f centered at the (E) (E) ∑ 〈 〉 𝒒 𝑓 𝑔 𝐸 𝑛 ,𝐸 (𝐸 ) (𝐸 ) 〈 〉 〈 〉 band edge, where the peak PF occurs as 𝒒 = . Finally, we compute 𝒒 by 𝑓 𝑔 𝐸 (𝐸 ) (𝐸 ) 〈 〉 averaging the 𝒒 over the bands. 〈 〉 𝒒 𝘀 2 𝑟 The correlation between the peak PF and the descriptor containing 𝑞 𝘀 , as ⁄ , is 𝑟 2 𝐷 𝑚 o cond shown in Figure 7b, with re = 0.88 and rh = 0.82. Effectively, this indicates that to achieve high correlation for p-type materials, the descriptor needs to include features of the IIS screening details, which reflect on the size of the constant energy surfaces. In bands with many valleys separated by low energy maxima, the constant energy surfaces can be very large and connect points very far in the BZ, even for energies not so far from the band edge. In this situation, the scattering vector can be very large. When the scattering vector becomes very large, it can dominate screening and play an important role in conductivity, the PF, and materials’ ranking. The extraction of such average scattering vector, however, might not be straightforward, and it certainly defeats the purpose of having descriptors that are easily constructed on simple parameters. A simple way to capture the size of the energy surface is to consider the number of valleys in each band. The larger the number of valleys (from the same band), the wider the surfaces are. For this, 27 we count the number of valleys in a range of 0.1 eV from the band edge. Essentially, this says that large constant energy surfaces reduce the IIS strength by improving the charge screening, and this boosts the PF. Figure 7. (a) Schematic of a regular parabolic band, in blue, and of a warped band with two (or in general many) valleys, in red, with the ionized impurity scattering vector 𝒒 indicated. In the case of warped, irregular bands, with many valleys, 𝒒 can be much larger (at the same energy), increasing the carrier screening effectiveness and decreasing the IIS scattering rate. (b) The 〈 〉 averaged 𝒒 is added to the descriptor of Figure 6e, together with the dielectric constant to produce a generic descriptor that can be employed in different materials situations. The Pearson correlation for the hole transport cases (note the valence bands are strongly warped) is improved. (c) The average number of valleys per band counted in the 0.1 eV range from the band edge is used as an easily accessible parameter to mimic the surface width and average IIS scattering exchange vector. The Pearson correlation is improved compared to Figure 6e as well. 𝑛̃𝘀 𝑣 𝑟 The correlation between the PF and the upgraded descriptor , is shown in Figure 𝐷 𝑚 o cond 7c, where 𝑛̃ is the average number of valleys per band. The descriptor is now very satisfactory for 28 both n-type and p-type materials (r = 0.89 and r = 0.90) and can be easily used instead of the e h scattering vector. Note that including the number of valleys will make the descriptor even stronger in the general case where inter-valley scattering is suppressed, independent of the scattering mechanism. However, in the half-Heusler materials we consider the inelastic optical phonons are dominant and inter-valley/inter-band scattering could be significant. The average number of valleys 𝑛 𝘀 𝑣 𝑟 per band, 𝑛̃, however, can be changed to the total number of valleys 𝑛 as with 𝑣 𝑣 2 𝐷 𝑚 o cond almost identical correlations outcomes (which reflects to the validity of the descriptor earlier in Figure 7 as well). Importantly, however, it is understood that the number of valleys is essential in cases where inter-valley scattering is suppressed for the dominant scattering mechanisms. In materials with wide constant energy surfaces, and under IIS dominant transport (as in the optimal cases for TE materials), due to enhanced screening we are in the regime of inter-valley suppression, and the use of 𝑛 is imperative. We underline that in NbCoSn, one of the best p-type materials we computed, the valence band has the highest number of valleys amongst the materials shortlist considered. In particular, the valence bands has a valley at the W point, which is 6-fold degenerate (24 W points shared between 4 adjacent zones) and relatively far from the other high-symmetry points (except from the X point). Thus having valleys at W is very efficient in suppressing inter- valley scattering, especially when there are no valleys in X. This supports recent findings that suggest that having valleys in W is helpful in gaining high PF. 3.2.3 Validity of the new descriptors at higher temperatures The analysis up to now considered only room temperature, T = 300 K. It is necessary to evaluate the performance of the descriptors for higher temperatures, however, since the majority of TE 29 materials are more efficient at higher temperatures. Figure 8 shows the Pearson correlation coefficients for the maximum PFs of the materials under investigation for T = 300 K, 600 K, and 900 K. Several combinations of the parameters that compose the final descriptor are separately analyzed. The intent is to highlight which has the most dominant role for electrons (nearly parabolic bands) in Figure 8a and Figure 8b, and holes (complex multiple-valleys bands) in Figure 8c and Figure 8d. In Figure 8a and Figure 8c (first column) we show results for phonon-limited transport, and in Figure 8b and Figure 8d for phonon plus IIS transport conditions. The conductivity effective mass 𝑚 is noted as 𝑚 in the figures for brevity. cond c In the phonon-limited case, for n-type materials with (nearly) parabolic bands (Figure 8a), at 300 K the descriptor 1/D performs well and is improved with increasing temperature (red-dotted line). On the other hand, although the performance of 1/mc is better at 300 K compared to 1/Do , with re = 0.74 (black-squared line), its performance decreases significantly with temperature. This is possibly because at higher temperatures the non-parabolicity and warping of the bands at higher energies become more relevant. For p-type materials with large and warped constant energy surfaces, the 1/𝑚 is not a good descriptor at any temperature with rh < 0.1 (Figure 8c), possibly for the same reason as for n-type at high temperatures. The 1/Do is a strong parameter for p-type phonon-limited transport, and in both carrier-type cases the combination of 1/D m is the dominant o c one with re = 0.96 and rh = 0.95 (green-triangle lines). In the phonon plus IIS transport case for n-type nearly parabolic bands (Figure 8b), the 1/D descriptor (red dots) provides a fair correlation, which improves with temperatures. When enriched with 1/mc, the 1/Do mc descriptor performs well across temperatures with re = 0.86 at 300 K and re = 0.83 at 900 K (green-triangle line). A further improvement occurs at 300 K by incorporating ε and the number of valleys, nv, leading to nvεr/ Do mcond with re = 0.90 (blue pentagons). Its strength 30 2 weakens at medium temperatures. Still, however, n ε / D m is retained at acceptable levels r = v r o cond e 0.84 at 600 K and 0.78 at 900 K. The 1/Do mc descriptor retains its strength for the warped highly non-parabolic p-type bands (Figure 8d) with r = 0.78 at 300 K and r = 0.94 at 900 K (green triangles). This originates from h h the synergy between 1/mc and 1/Do , which individually perform very poorly. In this case, an additional benefit to the descriptor is the addition of the average exchange vector <q> , which improves the r to 0.82 at T = 300 K up to 0.94 when T = 900 K. This signifies the importance of suppressing intervalley scattering in the presence of the dominant IIS. Note that this is important in the p-type materials we investigate, as their bandstructures consist of large and warped constant energy surfaces. The exchange vector can be effectively substituted by the band-averaged number of valley (see the Supporting Information for details). However, when we use the total number of valleys (blue pentagon) the descriptor loses its performance because of the low performance of n alone (light blue diamonds). This signals that more than the number of valleys itself, it is their distance in the BZ which matters. Thus, many valleys placed far away from each other in the BZ signify a favourable situation, whereas many valleys in the same or nearby BZ points represents a non-favourable situation. Moreover, the DOS effective mass of the materials with higher number of valleys increases with temperature, signalling that when inter-valley scattering is strong, it takes away possible positive effect of high degeneracy. It turns out, however, that the n ε /D m descriptor v r o cond (blue pentagons) is the best trade-off between simplicity and accuracy also for the p-type case, with 𝑛 𝘀 𝑣 𝑟 r > 0.75 over the whole temperature range. Thus, (blue pentagon lines) is the only h ⁄ 𝐷 𝑚 o cond descriptor that provides good to very high correlation results for most cases, across carrier type, scattering mechanisms, and temperatures. Figure 8. Power factor (PF) peak value Pearson correlations coefficients with multiple descriptors, computed at temperatures of 300 K, 600 K, and 900 K in the case of phonon-limited transport (first column) and phonon plus IIS transport (second column). The first row, (a-b), shows results for n- type materials and the second row, (c-d), for p-type materials. 3.3 Bipolar transport The bandgap is historically one of the most relevant features to signal TE material performance. However, the direct consideration of the energy gap size in the ranking of TE materials does not lead to straightforward conclusions. The compounds we consider have bandgaps ranging from 0.19 eV to 1.14 eV, which makes it necessary to consider the intrinsic thermally generated carriers and the minority carriers as well, i.e. bipolar transport effects, especially at elevated temperatures. 32 In this case, the transport calculation should take into account both conduction and valence bands combined, and this requires a slight modification in the descriptors, as we show below. A simulation difference between bipolar and unipolar transport is the different shift of the Fermi level with temperature, which influences the transport calculations significantly. In the unipolar case, as the temperature increases, the Fermi level shifts down (for n-type, for example), in order for charge neutrality to be satisfied, such that the mobile carrier density equals the doping density. In the bipolar case, for charge neutrality we need to consider also the intrinsic thermally generated carriers and the minority carriers as well, whose concentration increases with temperature. The increased density of both types of carriers leads to a reduced Fermi level shift compared to the unipolar case. Importantly, in this case both electrons and holes contribute to the screening of 34, 64 ionized impurities in an additive fashion. The correlation between the peak PF computed under unipolar and bipolar transport considerations for the phonon plus IIS case is shown in Figure 9a at T = 900 K where bipolar effects are relevant. The dashed line in the figure has a slope of unity and serves as guide-for-the- eye. The peak PF values feature an extremely high correlation, r > 0.9. At the lower T = 300 K the correlation is complete (not shown). It is also almost perfect in the case of 900 K and phonon- limited scattering. The magnitude of the PF, is of course reduced in the low bandgap materials, however, this does not affect the materials ranking, except for ScNiBi, which has the smallest bandgap amongst the investigated compounds. This is because in the materials we examine, the change in the PF due to bipolar effects is not larger in general compared to the differences in the PF between materials. In Figure 9b and Figure 9c we show the Pearson correlations for several descriptors at different temperatures for n-type and p-type materials, respectively. It is clear that the same descriptors that provide optimal correlations for the unipolar case, are also the most relevant ones for the bipolar 33 1 case, i.e. the ones involving for n-type (re > 0.6) and the ones involving 𝑛 for p-type 2 𝑣 𝐷 𝑚 o cond (rh > 0.7). The temperature trends are also similar in several cases. The descriptor developed earlier, 〈 〉 𝒒 𝘀 which provided high performance for both carrier types, ⁄ still performs 𝐷 𝑚 o cond exceptionally well for the p-type case, with Pearson correlation coefficients of rh ~ 0.85 at high temperatures. For n-type, however, it drops to r ~ 0.45 at 900 K. The observation is similar for the 𝑛 𝘀 𝑣 𝑟 descriptor, which results in Pearson correlation coefficients of r = 0.73, and r = h e 𝐷 𝑚 o cond 0.82 at 900 K. Thus, for bipolar transport in the more realistic case of phonon plus IIS transport, a modification to this descriptor is not needed. However, a simple modification is generally proposed in the literature by using the square of the reduced band gap, 𝐸 = , to capture the reduction to 𝑘 𝑇 the Seebeck coefficient with the band gap (squared as the Seebeck appears squared in the PF). 𝐸 𝑛 𝘀 G v 𝑟 The descriptor then becomes , but its Pearson correlation values at 900 K drop 𝐷 𝑚 𝑜 cond to re = 0.21 for n-type and rh = -0.43 for p-type, as indicated by the center-dotted circles shown in Figure 9b,c, respectively. The inclusion of the band gap in this latter descriptor does not lead to a straightforward improvement, as already indicated in the literature, and the common consideration of it requires further studies. 34 Figure 9. (a) Correlation between power factor (PF) peak values computed under bipolar and unipolar transport considerations for phonon plus IIS transport at 900 K. (b) and (c) PF peak value Pearson correlation coefficients with multiple descriptors, computed at 300 K, 600 K, and 900 K in the case of phonon plus IIS transport for n-type and p-type materials, respectively. 3.4 Discussion We have identified promising material candidates in section 3.1, both n-type and p-type. For n-type the more promising material candidates are NbFeSb, HfNiSn, YNiBi and ScNiBi with peak PFs > 15 mW/(m·K ). Other than NbCoSn, these are usually p-type TE materials so the introduction of donor impurities could present a materials science challenge. The most promising p-type materials are ScNiBi, NbCoSn, NbFeSb and TiNiSn with PFs ≥ 12 mW/(m·K ). For the former, the use of acceptor dopants has been demonstrated. The latter is poorly studied, but it is in principle p-type dopable with Sn in the Bi lattice site. The descriptors developed indicate why these compounds are the more promising among the ones we have studied. They have a good combination of a low optical deformation potential, a low conductivity effective mass, and a high dielectric constant. In the case of charge transport in the warped, multi-valley, valence bands, their many valleys for transport also result in large constant energy surfaces with large scattering vectors, which improve the charge screening and increase the conductivity. The NbCoSn is interesting as it performs quite well for both n- and p-type. Both can be experimentally realized by substituting Sn with Sb or Bi for n-type, and by substituting Nb with Ti 39, 40, 51 or Zr and Co with Fe for p-type. This would enable a TEG with the n-type and p-type legs from the same material, which could improve the TEG’s mechanical stability and life endurance. 35 We have shown the importance of considering electron scattering with ionized dopants in power 𝑛 𝘀 𝑣 𝑟 factor studies, and proposed the descriptor. Low deformation potentials, low 𝐷 𝑚 o cond conductivity effective masses, large constant energy surfaces through many valleys placed far from each other in the reciprocal space, and high dielectric constants, can compose a material properties ‘shopping list’ (experimentally and computationally). While the dielectric constant is usually computed in ab-initio calculations and the conductivity effective mass and number of valleys can be easily extracted from the 3D bandstructure, the deformation potentials are rarely computed in first principle calculations. For example, we register only two, very recent attempts to compute 9, 24, 25, 65 from first principles the deformation potentials in half-Heusler alloys. The situation for 66, 67 other TE material families is similar. We suggest that in the predictive search for better thermoelectrics, deformation potential calculations are also performed – at least for regular nearly parabolic bandstructure material. It is worth to underline the importance of the deformation potential method. In the deformation potential theory, the electron-phonon coupling is described by a single parameter, the deformation potential. This theory can have some limitations in respect of a precise evaluation of the coupling strength for each possible pair of initial and final state and each possible phonon, especially for complex band materials. Nevertheless, the latter treatment is computationally extremely expensive and, at present, can be successfully used for a relatively fast and automated computation of the transport properties only in 2D materials and with the restriction to certain portions of the BZ, thanks to the reduced number of employed k-points. The use of deformation potentials can be a significant step forward compared to constant relaxation times. Besides, the proper and feasible computation of the deformation potential values can be extremely useful in the materials ranking 36 and pre-screening. Indeed, the predictive evaluation of the deformation potentials, without available 9, 65 – 67 experimental data for comparison, is an on-going effort in the thermoelectric community. Finally, we note that there is significant effort undertaken in order to achieve band alignment in complex bandstructure materials for high power factors. Recent works show that this is beneficial only in the case where the bands brought closer to the band edge are lighter, and when inter-valley 27, 68 scattering is weak. This is very well aligned with what our descriptors suggest, i.e. sharp bands with light conductivity effective masses and large energy surfaces with strong charge screening, which effectively leads to weakened inter-valley IIS scattering. Thus, the descriptors we developed can also be applicable in band alignment studies for improving the thermoelectric power factor. 4. Conclusions In summary, we have developed a computational method that combines energy dependent electron- phonon and electron-ionized impurity scattering mechanisms with DFT electronic structures. We then investigated the thermoelectric power factor performance of a group of 14 half-Heuslers, for both n-type and p-type transport. The comparison of the materials’ performances under the constant relaxation time approximation and energy dependent scattering assumptions shows absence of any correlation in the materials ranking. We tested the validity of a set of descriptors currently used in the literature for materials screening exercises, and have developed new, highly improved descriptors for the thermoelectric power factor, based on our detailed scattering treatment. Our 𝑛 𝘀 𝑣 𝑟 proposed descriptor is based on number of valleys, the dielectric constant, the 𝐷 𝑚 𝑜 cond conductivity effective mass, and the deformation potential for the dominant electron-phonon 37 process. In the case of strong bipolar effects, it can be modified with the inclusion of the reduced 𝐸 𝑛 𝘀 G 𝑣 𝑟 bandgap to . 𝐷 𝑚 𝑜 cond We strongly suggest that every thermoelectric materials screening study should consider in more detail: i) the energy and wave-vector dependence of the scattering processes, ii) the ionized impurity scattering, and iii) the volume that the bands occupy in the Brillouin Zone. We further show, in the case of n-type half-Heuslers, which have nearly parabolic bands, a highly relevant and simplest descriptor is 1/𝐷 𝑚 . For warped non-parabolic multi-band/valley bandstructures, as the p-type 𝑜 cond half-Heuslers, the most important transport aspect is the ability of the carriers to screen the dopant impurities, which leads to inter-valley scattering suppression. Only then, the different valleys can contribute to transport independently and improve the PF. The number of valleys, 𝑛 , can then be used as a single parameter descriptor, as it i) reflects the average size of the constant energy surface, and the larger it is, the stronger the IIS screening is, which allows the valleys to ii) independently contribute to transport without increasing scattering for carriers from other valleys. The insight gained from this work can serve as guidance to understand experimental observations. The proposed descriptors can be useful materials screening tools to advance the search for high performance thermoelectric materials. ASSOCIATED CONTENT Supporting Information. Tables with scattering input parameters and extracted effective masses, procedure to extract the effective masses. This material is available free of charge via the Internet at http://pubs.acs.org. 38 AUTHOR INFORMATION Corresponding Author *Patrizio.Graziosi@warwick.ac.uk, Patrizio.Graziosi@gmail.com Author Contributions The manuscript was written through contributions of all authors. Funding Sources European Commission under the Grant agreement ID: 788465 (GENESIS - Generic semiclassical transport simulator for new generation thermoelectric materials) and European Research Council (ERC) under the European Union's Horizon 2020 Research and Innovation Programme (Grant Agreement No. 678763). ACKNOWLEDGMENT This work has received funding from the Marie Skłodowska-Curie Actions under the Grant agreement ID: 788465 (GENESIS - Generic semiclassical transport simulator for new generation thermoelectric materials) and from the European Research Council (ERC) under the European Union's Horizon 2020 Research and Innovation Programme (Grant Agreement No. 678763). The authors thank Laura De Sousa Oliveira, Dhritiman Chakraborty, Alberto Riminucci and Alek Dediu for reading the manuscript. 39 REFERENCES (1) Gorai P., Stevanović V., Toberer E. S., Computationally guided discovery of thermoelectric materials, Nat. Rev. Mater. 2017, 2, 17053. (2) Samsonidze G., Kozinsky B., Accelerated Screening of Thermoelectric Materials by First‐ Principles Computations of Electron–Phonon Scattering, Adv. Energy Mater. 2018, 8, (3) Thermoelectric Roadmap, 2018, editors: Robert Freer and Anthony Powell. (4) Future Market Insights, Report ID REP-GB-3652. (5) Madsen G. K. H., Automated Search for New Thermoelectric Materials: The Case of LiZnSb, J. Am. Chem. Soc. 2006, 128, 12140. (6) Wang S., Wang Z., Setyawan W., Mingo N., Curtarolo S., Assessing the Thermoelectric Properties of Sintered Compounds via High-Throughput Ab-Initio Calculations, Phys. Rev. X 2011, 1, 021012. (7) Yan J., Gorai P., Ortiz B., Miller S., Barnett S. A., Mason T., Stevanović V., Toberer E. S., Material descriptors for predicting thermoelectric performance, Energy Environ. Sci. 2015, 8, (8) Pizzi G., Cepellotti A., Sabatini R., Marzari N., Kozinsky B., AiiDA: automated interactive infrastructure and database for computational science, Comput. Mater. Sci. 2016, 111, 218. (9) Zhou J., Zhu H., Liu T. H., Song Q., He R., Mao J., Liu Z., Ren W., Liao B., Singh D. J., Ren Z., Chen G. , Large thermoelectric power factor from crystal symmetry-protected non-bonding orbital in half-Heuslers, Nat. Commun. 2018, 9, 1721. (10) Berland K., Shulumba N., Hellman O., Persson C., Løvvik O. M., Thermoelectric transport trends in group 4 half-Heusler alloys, J. Appl. Phys. 2019, 126, 145102. 40 (11) Mahan G. D., in Solid State Physics, edited by H. Ehrenreich and F. Spaepen (Academic Press, New York, 1998), Vol. 51, p. 81. (12) Biswas K., He J., Blum I. D., Wu C. I., Hogan T. P., Seidman D. N., Dravid V. P., Kanatzidis M., High-performance bulk thermoelectrics with all-scale hierarchical architectures, Nature 2012, 489, 414. (13) Zhao H., Cao B., Li S., Liu N., Shen J., Li S., Jian J., Gu L., Pei Y., Snyder G. J., Ren Z., Chen X., Engineering the Thermoelectric Transport in Half‐Heusler Materials through a Bottom‐Up Nanostructure Synthesis, Adv. Energy Mater. 2017, 7, 1700446. (14) Narducci D., Do we really need high thermoelectric figures of merit? A critical appraisal to the power conversion efficiency of thermoelectric materials, Appl. Phys. Lett. 2011, 99, 102104 (15) Gaultois M. W., Sparks T. D., Borg C. K. H., Seshadri R., Bonificio W. D., Clarke D. R., Data-Driven Review of Thermoelectric Materials: Performance and Resource Considerations, Chem. Mater. 2013, 25, 2911. (16) Neophytou N., Karamitaheri H., Kosina H., Atomistic calculations of the electronic, thermal, and thermoelectric properties of ultra-thin Si layers, J. Comput. Electr. 2013, 12, 611. (17) Beretta D., Neophytou N., Hodges J. M., Kanatzidis M. G., Narducci D., Martin-Gonzalez M., Beekman M., Balke B., Cerretti G., Tremel W., Zevalkink A., Hofmann A. I., Müller C., Dörling B., Campoy-Quiles M., Caironi M., Thermoelectrics: From history, a window to the future, Mater. Sc. Eng. R 2019, 138, 100501. (18) 18. http://www.mrl.ucsb.edu:8080/datamine/thermoelectric.jsp, Access March 2020. (19) 19. Zeier W. G., Zevalkink A., Gibbs Z. M., Hautier G., Kanatzidis M. G., Snyder G. J., Thinking Like a Chemist: Intuition in Thermoelectric Materials, Angew. Chem., Int. Ed. 2016, 55, 6826. 41 (20) Sohier T., Campi D., Marzari N., Gibertini M., Mobility of two-dimensional materials from first principles in an accurate and automated framework, Phys. Rev. Materials 2018, 2, 114010. (21) Pizzi G., Volja D., Kozinsky B., Fornari M., Marzari N., BoltzWann: A code for the evaluation of thermoelectric and electronic transport properties with a maximally-localized Wannier functions basis, Comput. Phys. Commun. 2014, 185, 422. (22) Markov M., Rezaei E., Sadeghi S. N., Esfarjani K., Zebarjadi M., Thermoelectric properties of semimetals, Phys. Rev. Materials 2019, 3, 095401. (23) Graziosi P., Kumarasinghe C., Neophytou N., Impact of the scattering physics on the power factor of complex thermoelectric materials, J. Appl. Phys. 2019, 126, 155701. (24) Jia, T., Feng, Z., Guo, S., Zhang, X., Zhang, Y., Screening Promising Thermoelectric Materials in Binary Chalcogenides through High-Throughput Computations. ACS Appl. Mater. Interfaces 2020, 12, 11852. (25) Xi, L., Pan, S., Li, X., Xu, Y., Ni, J., Sun, X., Yang, J., Luo, J., Xi, J., Zhu, W., Li, X., Jiang, D., Dronskowski, R., Shi, X., Snyder, G. J., Zhang, W., Discovery of High- Performance Thermoelectric Chalcogenides through Reliable High-Throughput Material Screening. J.Am. Chem. Soc. 2018, 140, 10785. (26) Giannozzi P., Baroni S., Bonini N., Calandra M., Car R., Cavazzoni C., Ceresoli D., Chiarotti G. L., Cococcioni M., Dabo I., Dal Corso A., de Gironcoli S., Fabris S., Fratesi G., Gebauer R., Gerstmann U., Gougoussis C., Kokalj A., Lazzeri M., Martin-Samos L., Marzari N., Mauri F., Mazzarello R., Paolini S., Pasquarello A., Paulatto L., Sbraccia C., Scandolo S., Sclauzero G., Seitsonen A. P., Smogunov A., Umari P., Wentzcovitch R. M., QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials, J. Phys.: Condens. Matter 2009, 21, 395502; Giannozzi P., Andreussi O., Brumme T., Bunau O., Buongiorno Nardelli M., Calandra M., Car R., Cavazzoni C., Ceresoli D., Cococcioni M., 42 Colonna N., Carnimeo I., Dal Corso A., de Gironcoli S., Delugas P., DiStasio R. A., Ferretti A., Floris A., Fratesi G., Fugallo G., Gebauer R., Gerstmann U., Giustino F., Gorni T., Jia J., Kawamura M., Ko H. Y., Kokalj A., Küçükbenli E., Lazzeri M., Marsili M., Marzari N., Mauri F., Nguyen N. L., Nguyen H. V., Otero-de-la-Roza A., Paulatto L., Poncé S., Rocca D., Sabatini R., Santra B., Schlipf M., Seitsonen A. P., Smogunov A., Timrov I., Thonhauser T., Umari P., Vast N., Wu X., Baroni S., Advanced capabilities for materials modelling with Quantum ESPRESSO, J. Phys.: Condens. Matter 2017, 29, 465901. (27) Kumarasinghe C., Neophytou N., Band alignment and scattering considerations for enhancing the thermoelectric power factor of complex materials: The case of Co-based half-Heusler alloys, Phys. Rev. B 2019, 99, 195202. (28) He R., Kraemer D., Mao J., Zeng L., Jie Q., Lan Y., Li C., Shuai J., Kim H. S., Liu Y., Broido D., Chu C. W., Chen G., Ren Z., Achieving high power factor and output power density in p- type half-Heuslers Nb1-xTixFeSb, Proc. Natl. Acad. Sci. 2016, 113, 13576. (29) Fischetti M. V., Laux S. E., Band structure, deformation potentials, and carrier mobility in strained Si, Ge, and SiGe alloys, J. Appl. Phys. 1996, 80, 2234. (30) Neophytou N. and Kosina H., Atomistic simulations of low-field mobility in Si nanowires: Influence of confinement and orientation, Phys. Rev. B 2011, 84, 085313. (31) Neophytou N. and Kosina H., Effects of confinement and orientation on the thermoelectric power factor of silicon nanowires, Physical Review B 2011, 83, 245305. (32) Jain A., Ong S. P., Hautier G., Chen W., Richards W. D., Dacek S., Cholia S., Gunter D., Skinner D., Ceder G., Persson K. A., Commentary: The Materials Project: A materials genome approach to accelerating materials innovation, APL Materials 2013, 1, 011002. 43 (33) Wang D., Wang G., First-principles study the elastic constant, electronic structure and thermoelectric properties of Zr Hf NiPb (x = 0, 0.25, 0.5, 0.75, 1) , Phys. Lett. A 2017, 381, 1−x x (34) Nag B. R., Electron Transport in Compound Semiconductors Vol. 11 Springer Series in Solid- State Sciences (Springer-Verlag Berlin Heidelberg New York, New York, 1980). (35) Brown D. M. and Bray R., Analysis of Lattice and Ionized Impurity Scattering in p-Type Germanium, Phys. Rev. 1962, 127, 1593. (36) Fischetti M. V., Vandenberghe W. G. (2016) Scattering with Ionized Impurities, in: Advanced Physics of Electron Transport in Semiconductors and Nanostructures, pages 315-325. Graduate Texts in Physics book series, Springer, Cham. (37) Rausch E., Balke B., Ouardi S., Felser C., Enhanced thermoelectric performance in the p-type half-Heusler (Ti/Zr/Hf)CoSb Sn system via phase separation, Phys. Chem. Chem. Phys. 0.8 0.2 2014, 16, 25258. (38) Kimura Y., Kuji T., Zama A., Lee T., Mishima Y., Thermoelectric properties of half-Heusler compounds N-type MNiSn and P-type MPtSn (M= Hf, Zr) , MRS Online Proc. Libr. Arch. 2006, 980. (39) Ferluccio D. A., Smith R., Buckman J., Bos J.-W., Impact of Nb vacancies and p-type doping of the NbCoSn–NbCoSb half-Heusler thermoelectrics, Phys. Chem. Chem. Phys. 2018, 20, (40) He R., Huang L., Wang Y., Samsonidze G., Kozinsky B., Zhang Q., Ren Z., Enhanced thermoelectric properties of n-type NbCoSn half-Heusler by improving phase purity, APL Mater. 2016, 4, 104804. (41) Oestreich J., Probst U., Richardt F., Bucher E., Thermoelectrical properties of the compounds VIII VIII VIII ScM Sb and YM Sb (M = Ni, Pd, Pt), J. Phys.: Condens. Matter 2003, 15, 635. 44 (42) Muta H., Kanemitsu T., Kurosaki K., Yamanaka S., High-temperature thermoelectric properties of Nb-doped MNiSn (M= Ti, Zr) half-Heusler compound, J. Alloys Compd. 2009, 469, 50. (43) Berry T., Fu C., Auffermann O. G., Fecher G. H., Schnelle W., Serrano-Sanchez F., Yue Y., Liang H., Felser C., Enhancing thermoelectric performance of TiNiSn half-Heusler compounds via modulation doping, Chem. Mater. 2017, 29, 7042. (44) Li S., Zhao H., Li D., Jin S., Gu L., Synthesis and thermoelectric properties of half-Heusler alloy YNiBi, J. Appl. Phys. 2015, 117, 205101. (45) Zhu H., He R., Mao J., Zhu Q., Li C., Sun J., Ren W., Wang Y., Liu Z., Tang Z., Sotnikov A., Wang Z., Broido D., Singh D. J., Chen G., Nielsch K., Ren Z., Discovery of ZrCoBi based half Heuslers with high thermoelectric conversion efficiency, Nat. Commun. 2018, 9, 2497. (46) Sekimoto T., Kurosaki K., Muta H., Yamanaka S., High-Thermoelectric Figure of Merit Realized in p-Type Half-Heusler Compounds: ZrCoSnxSb1-x, Jpn. J. Appl. Phys. 2007, 46, L673. (47) Xie H., Wang H., Fu C., Liu Y., Snyder G. J., Zhao X., Zhu T., The intrinsic disorder related alloy scattering in ZrNiSn half-Heusler thermoelectric materials, Sci. Rep. 2014, 4, 6888. (48) Muta H., Yamaguchi T., Kurosaki K., Yamanaka S., Thermoelectric properties of ZrNiSn based half Heusler compounds, 24th International Conference on Thermoelectrics, ICT 2005, (49) Zilber T., Cohen S., Fuks D., Gelbstein Y., TiNiSn half-Heusler crystals grown from metallic flux for thermoelectric applications, J. Alloys Compd. 2019, 781, 1132. (50) Armakavicius N., Stanishev V., Knight S., Kuhne P., Schubert M, Darakchieva V., Electron effective mass in In0.33Ga0.67N determined by mid- infrared optical Hall effect, Appl.Phys. Lett. 2018, 112, 082103. 45 (51) Kimura Y., Tamura Y., Kita T., Thermoelectric properties of directionally solidified half- Heusler compound NbCoSn alloys, Appl. Phys. Lett. 2008, 92, 012105. (52) Huang L., Zhang Q., Yuan B., Lai X., Yan X., Ren Z., Recent progress in half-Heusler thermoelectric materials, Materials Research Bulletin 2016, 76, 107. (53) Wu T., Jiang W., Li X. Y., Zhou Y. F., Chen L. D., Thermoelectric properties of p-type Fe- doped TiCoSb half-Heusler compounds, J. Appl. Phys. 2007, 102, 103705. (54) Gaultois M. W., Oliynyk A. O., Mar A., Sparks T. D., Mulholland G. J., Meredig B., Perspective: Web-based machine learning models for real-time screening of thermoelectric materials properties, APL Mater. 2016, 4, 053213. (55) https://www.materialscloud.org/ (56) https://materialsproject.org/ (57) http://oqmd.org/ (58) https://nomad-coe.eu/ (59) Tang Y., Gibbs Z. M., Agapito L. A., Li G., Kim H.-S., Buongiorno Nardelli M., Curtarolo S., Snyder G. J., Convergence of multi-valley bands as the electronic origin of high thermoelectric performance in CoSb skutterudites, Nature Materials 2015, 14, 1223. (60) Rudderham C., Maassen J., Analysis of simple scattering models on the thermoelectric performance of analytical electron dispersions, J. Appl. Phys. 2020, 127, 065105. (61) Ravich Yu. I., Efimova B. A., Tamarchenko V. I., Scattering of Current Carriers and Transport Phenomena in Lead Chalcogenides, physica statu solidi (b) 1971, 43, 11. (62) Dylla M. T. , Dunn A., Anand S., Jain A., Snyder G. J., Machine Learning Chemical Guidelines for Engineering Electronic Structures in Half-Heusler Thermoelectric Materials, Research 2020, 8, 6375171. 46 (63) Foster S., Neophytou N., Doping optimization for the power factor of bipolar thermoelectric materials, J. Electron. Mater. 2019, 48, 1889. (64) Combescot M., Combescot R., Bok J., Screening of the Electron-Phonon Interaction in Semiconductors, EPL (Europhysics Letters) 1986, 2, 31. (65) Naydenov G. A., Hasnip P. J., Lazarov V. K., Probert M. I. J., Huge power factor in p-type half-Heusler alloys NbFeSb and TaFeSb, J. Phys.: Materials 2019, 2, 035002. (66) Cao J., Querales-Flores J. D., Murphy A. R., Fahy S., Savić I., Dominant electron-phonon scattering mechanisms in n-type PbTe from first principles, Phys. Rev. B 2018, 98, 205202 (67) Murphy A. R., Murphy-Armando F., Fahy S., Savić I., Acoustic deformation potentials of n- type PbTe from first principles, Phys. Rev. B 2018, 98, 085201. (68) Kumarasinghe C., Neophytou N., “Thermoelectric Power Factor Under Strain-Induced Band- Alignment in the Half-Heuslers NbCoSn and TiCoSb”, in: E. Levchenko, Y. Dappe, G. Ori (eds) Theory and Simulation in Physics for Materials Applications. Springer Series in Materials Science, 2020, vol 296. Springer, Cham (69) Neophytou N., Kosina H., Large Enhancement in Hole Velocity and Mobility in p-Type  and  Silicon Nanowires by Cross Section Scaling: An Atomistic Analysis, Nano Lett. 2010, 10, 4913. (70) Rahman A., Guo J., Datta S., Lundstrom M. S., Theory of ballistic nanotransistors, IEEE Trans. Electron Devices 2003, 50, 1853. Supporting Information Patrizio Graziosi*, Chathurangi Kumarasinghe, Neophytos Neophytou * firstname.lastname@example.org Table S1. The scattering parameters as used in the simulations presented in this work, extracted 9, 32, 33 from the literature. a) a) Material us ADP ODP ℏ r 10 kg/m × [m/s × ] [eV] [eV/m ×10 ] [eV] HfCoSb 9.54 5.64 0.2 / 0.4 1.4 / 1.85 0.028 17.51 HfNiSn 10.3 3.39 0.1 / 0.7 1.8 / 2 0.028 20.88 NbCoSn 8.43 5.36 0.2 / 0.5 2.6 / 2.16 0.034 22.6 NbFeSb 8.45 3.99 1 / 0.8 1.6 / 2.1 0.036 22.99 ScNiBi 8.48 3.1 0.8 / 0.2 2.2 / 1.4 0.028 29 ScNiSb 6.6 3.9 0.7 / 0.4 2.3 / 1.8 0.032 19.66 TiCoSb 7.4 4.04 1 / 0.5 1.75 / 2.2 0.036 19.09 TiNiSn 7.2 3.9 0.7 / 0.2 2.3 / 0.9 0.034 22.92 YNiBi 8.6 2.9 0.5 / 0.9 2.3 / 2.7 0.022 26.55 YNiSb 6.9 3.5 0.4 / 0.9 2.35 / 3.14 0.026 19.74 ZrCoBi 9.83 3.2 0.2 / 0.8 2.1 / 1.8 0.028 20.37 ZrCoSb 7.14 5.55 0.2 / 1 2.05 / 1.75 0.028 17.87 ZrNiPb 9.6 3.39 0.36/0.60 2.8/1.6 0.024 23.6 ZrNiSn 7.62 3.97 0.35 / 0.3 2.8/ 1.8 0.029 20.9 a) electrons / holes Table S2. Extracted bandstructure parameters for n-type materials. The effective masses are in units of the electron rest mass and the total number of valleys considers the valleys within 0.1 eV from the band edge. nb is the number of bands and 𝑛̃ the average number of valleys per band. Material mDOS mDOS mcond mcond mcond nv nb 𝑛 ̃ mDOS T = 300 K T = 900 K T = 300 K T = 600 K T = 900 K T = 600 K HfCoSb 4.31 3.99 3.8 1.42 1.33 1.26 3 1 3 HfNiSn 4.1 3.5 3.29 0.93 0.96 1.0 3 1 3 NbCoSn 3.64 3.68 3.8 1.33 1.29 1.32 6 2 3 NbFeSb 0.92 0.88 0.9 0.33 0.37 0.61 3 1 3 ScNiBi 3.12 2.83 2.69 0.75 0.73 0.72 3 1 3 ScNiSb 3.02 2.78 2.65 0.79 0.78 0.77 3 1 3 TiCoSb 4.97 4.92 5.05 2.12 2.03 2.03 3 1 3 TiNiSn 2.43 2.4 2.46 0.84 0.89 0.94 3 1 3 YNiBi 3.26 2.79 2.57 0.63 0.61 0.6 3 1 3 YNiSb 3.45 3.03 2.82 0.76 0.73 0.71 3 1 3 ZrCoBi 4.49 4.13 3.94 1.44 1.34 1.28 3 1 3 ZrCoSb 4.36 4.07 3.9 1.59 1.49 1.42 3 1 3 ZrNiPb 1.93 1.86 1.84 0.55 0.57 0.58 3 1 3 ZrNiSn 1.94 1.87 1.85 0.58 0.6 0.61 3 1 3 Table S3. Extracted bandstructure parameter for p-type materials. The effective masses are in units of the electron rest mass and the total number of valleys considers the valleys within 0.1 eV from the band edge. nb is the number of bands and 𝑛̃ the average number of valleys per band. Material mDOS mDOS mcond mcond mcond nv nb 𝑛 ̃ mDOS T = 300 K T = 900 K T = 300 K T = 600 K T = 900 K T = 600 K HfCoSb 7.8 8.31 8.58 1.44 1.55 1.73 3.67 11 3 HfNiSn 1.88 1.8 1.78 0.55 0.56 0.57 3 3 1 NbCoSn 5.79 6.67 7.2 1.08 1.2 1.34 16 3 5.33 NbFeSb 5.07 4.84 4.85 1.05 1.08 1.13 4 8 2 ScNiBi 1.92 1.71 1.64 0.44 0.49 0.52 3 3 1 ScNiSb 2.2 1.96 1.87 0.49 0.55 0.58 3 3 1 TiCoSb 7.35 8.6 9.16 1.58 1.88 2.14 3.67 11 3 TiNiSn 6.11 5.33 5.11 1.64 1.64 1.69 3 3 1 YNiBi 1.43 1.28 1.25 0.3 0.36 0.41 3 3 1 YNiSb 1.79 1.62 1.59 0.39 0.46 0.51 3 3 1 ZrCoBi 6.59 6.74 7.09 1.4 1.45 1.58 3.67 11 3 ZrCoSb 6.63 6.63 6.96 1.46 1.54 1.7 8 3 2.67 ZrNiSn 3.33 2.97 2.88 0.93 0.95 1.01 3 3 1 ZrNiPb 4.49 3.98 3.94 1.15 1.2 1.36 1 3 3 Extraction of the density of states and conductivity effective masses from bandstructure The data analysis and proposed descriptors make use of the DOS and conductivity effective masses m and m , respectively. These are effective quantities that include information from the entire DOS cond bandstructure, rather just specific valleys. To extract these, we consider non-degenerate conditions, where the Fermi level is placed in the bandgap, several kBT from the bands. Then the mDOS is the effective mass of an isotropic parabolic band that provides the same carrier density as the actual bandstructure. The m is the effective mass of an isotropic parabolic band that has the same cond average velocity as all that of the entire bandstructure. To compute mDOS we utilize the common expression for the carrier density in a non-degenerate semiconductor: 𝐸 −𝐸 𝑐 𝐹 𝑘 𝑇 𝑁 = 𝑁 𝑒 (S1) where 𝐸 is the band edge, 𝐸 is the Fermi level, 𝑘 is the Boltzmann constant, T is the temperature 𝑐 𝐹 B 3/2 𝑚 𝑘 𝑇 DOS B and 𝑁 = 2 ( ) is the effective density of states in the conduction band. This expression 2𝜋 ℏ is derived by approximating the Fermi Dirac statistics with Boltzmann statistics. The carrier density in the numerical DFT bandstructure is given by: 𝑁 = ∑ 𝑓 𝑑𝑉 (S2) 𝒌 ,𝑛 𝐸 𝒌 3 𝒌 ,𝑛 (2𝜋 ) where the sum runs over all the k points in the reciprocal unit cell (mesh used for the DFT calculations) and all the bands, 𝑓 is the Fermi-Dirac distribution evaluated at the energy of each 𝒌 ,𝑛 state (approximated by Boltzmann statistics in the non-degenerate limit), and 𝑑𝑉 is the volume element in the k space that depends only on the mesh. By equating the carrier density from equations (S1) and (S2) for the same Fermi level and band edge, we can extract the effective density of states in the conduction band, 𝑚 . The process is the same for the valence band. DOS The conductivity effective mass of the bandstructure, 𝑚 , is extracted from the so-called cond injection velocity 𝑣 of the carriers in non-degenerate conditions. As this method is derived for inj 51 70 MOSFET devices, it considers the ballistic current that the positive velocity states of the bandstructure will allow for a specific Fermi level in the non-degenerate conditions as: 1 2 𝐼 = 𝑒 𝑓 |𝑣 |𝑑𝑉 (S3) +,𝑛 𝒌 (𝐸 −𝐸 ) 𝒌 𝒌 𝑛 𝒌 𝐹 ,S 𝑛 𝑛 2 (2𝜋 ) where e is the electron charge, the ½ factor is to take into account only half of the states, and, |𝑣 | is the band carrier velocity in absolute terms in order to only account for positive velocities. Then, the injection velocity is extracted by dividing the positive going ballistic current with the carrier density (to get the average velocity per carrier): +,𝑛 𝑣 = (S4) inj 1 2 𝑒 ∑ 𝑓 𝑑𝑉 𝒌 (𝐸 −𝐸 ) 𝒌 3 𝑛 𝑛 2 𝒌 𝐹 ,S (2𝜋 ) where the denominator in equation (S4) is the charge occupying the positive velocity states alone. Once the injection velocity is extracted, the conductivity effective mass can be computed as: 2𝑘 𝑇 𝑚 = (S5) cond 2 𝑣 𝜋 inj 69, 70 which is the corresponding quantity for a simple parabolic band. This calculation is performed along the three Cartesian directions x, y, and z. Thus, three conductivity affective masses are extracted, 𝑚 , 𝑚 , and 𝑚 and a final average cond,𝑥 cond,𝑦 cond,𝑧 conductivity affective mass is then formed by 𝑚 = . The process cond −1 −1 −1 𝑚 +𝑚 +𝑚 cond,𝑥 cond,𝑦 cond,𝑧 is performed twice, once for the conduction and once for the valence bands. As an example, we show in Figure S1a the numerically extracted DOS for HfCoSb (blue lines) and the DOS provided by the extracted mDOS, assuming that it forms an isotropic, parabolic band (dashed-dot green line), indicating adequate match at low energies. In Figure S1b we show the HfCoSb valence band injection velocity using the numerical bandstructure (blue line), and the injection velocity for an isotropic parabolic band with the extracted 𝑚 value (dash-dot green cond line), again indicating adequate match in the low-density region. The relevance of 𝑚 becomes clear when comparing the TDF computed for the actual cond numerical bandstructure and for an isotropic parabolic band of an effective mass equal to the extracted 𝑚 , Figure S1c. The blue solid line (numerical bandstructure) and the green dashed- cond 52 dot line (parabolic band with 𝑚 ) are for the ADP-limited scattering mechanism. The lines cond overlap for a large energy region, indicating that 𝑚 captures the core transport information for cond ADP scattering in complex bandstructures. In the ODP case, the TDFs for the solid orange line (for the numerical bands) and the purple dash-dot line (for the parabolic band of 𝑚 ) also clearly cond overlap. The TDFs deviate far from the band onset, when the parabolicity is lost. Figure S1: Comparisons between quantities extracted from the full numerical bandstructure and the corresponding quantities computed from the extracted m and m for HfCoSb. (a) The DOS of HfCoSb, blue line, and the DOS DOS cond of parabolic bands with the isotropic m values, green dashed lines. (b) Injection velocity versus carrier concentration DOS for the positive going carriers of the HfCoSb valence band under ballistic conditions using the full bandstructure, blue solid line, and the parabolic isotropic band with the extracted effective mass value m , green dash-dotted line. (c) cond Transport distribution function versus energy for the HfCoSb valence band using the full bandstructure and the band formed with the extracted m . The case of acoustic phonons (ADP – Acoustic Deformation Potential) transport and cond optical phonons (ODP – Optical Deformation Potential), transport are plotted. The blue and yellow solid lines indicate the full bandstructure TDFs. The green and purple dash-dotted lines indicate the TDFs for the parabolic isotropic band having an effective mass of m . cond
Physics – arXiv (Cornell University)
Published: Jun 4, 2020
Access the full text.
Sign up today, get DeepDyve free for 14 days.