Computation-driven materials search for thermoelectric applications

The advancement of computational tools for material property predictions enables broad search of novel materials for various energy-related applications. However, challenges still exist in accurately predicting the mean free paths of electrons and phonons in a high throughput frame for thermoelectric property predictions, which largely hinders the computation-driven material search for novel materials. In this work, this need is eliminated under the small-grain-size limit, in which these mean free paths are restricted by the grain sizes within a bulk material. A new criterion for Z T evaluation is proposed for general nanograined bulk materials and is demonstrated with representative oxides. © The Author(s) 2017. Published by ECS. This is an open access article distributed under the terms of the Creative Commons Attribution 4.0 License (CC BY, http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse of the work in any medium, provided the original work is properly cited. [DOI: 10.1149/2.0141703jss] All rights reserved.

Solid-state thermoelectric (TE) devices have the ability to directly convert heat into electricity for power generation. 1 Despite its many energy-harvesting applications, the potential impact of TE technology is largely hindered by the relatively low performance of commercial materials. In physics, the effectiveness of TE materials is evaluated by their TE figure of merit (Z T ), defined as Z T = S 2 σT /k, where S, σ, k, and T represent Seebeck coefficient, electrical conductivity, thermal conductivity, and absolute temperature, respectively. Here k can be split into lattice contribution k L and electronic contribution k E . Within the same material, challenge exists in obtaining a low k but a high power factor (PF) S 2 σ. As the result, Z T s of commercial materials are still around unity after decades of research though Z T > 2 is preferred for TE to compete with other technologies. Along this line, nanostructured bulk materials, which are synthesized by hot pressing nanopowder into a bulk material, have been widely studied as an effective approach to improve Z T s of existing or novel materials. [2][3][4][5][6][7][8][9][10] By introducing nanostructured interfaces within a bulk material, the interdependent electron and phonon transport can be decoupled to dramatically reduce k L but still maintain S 2 σ, resulting in enhanced Z T . Unrestricted to conventional materials, this approach may also lead to high Z T s in unconventional TE materials with both a high S 2 σ and a high k L , such as bulk Si. 5,6 This will reach beyond conventional materials that heavily depend on toxic, rare, and expensive elements, e.g. Te used in Bi 2 Te 3 and PbTe. In addition, the state-of-the-art TE materials are mostly restricted to sub-1000 K due to their low melting points, poor thermal stability, and/or serious oxidation over 1000 K. This restricts high-temperature (>1000 K) TE power generation with much higher Carnot efficiency and thus more effective energy conversion though such heat sources are available in industrial applications. With particular attention to high-temperature applications, broad material search is in urgent need for novel TE materials using nontoxic, cheap, and abundant elements. Such material discovery can be largely accelerated by the Materials Genome approach that uses first-principles computation to predict the TE properties of a new material. [11][12][13][14][15][16] For computation-driven TE material search, challenge still exists in accurately predicting the energy-dependent electron bulk mean free path (MFP) l E,bulk and phonon bulk MFP l P,bulk for TE property predictions. For phonon transport, there exist first-principles-based studies for standard materials such as Si, [17][18][19] Ge, 18,19 GaN, 20 BAs, 21 and PbTe. 22 However, such calculations are intrinsically very timeconsuming and are limited to relatively simple material structures. In heavily doped TE materials, the strong point-defect scattering z E-mail: qinghao@email.arizona.edu; zhaohb@email.arizona.edu and electron scattering of phonons will add more uncertainties to the modeling. On the other hand, charge carriers are also affected by various scattering mechanisms, including acoustic/optical deformation potential scattering, ionized/neutral impurity scattering, nonpolar/polar optical phonon scattering, and piezoelectric scattering. 23 In this situation, accurate electron modeling requires identifying scattering mechanisms for a given material and fitting the experimental data to determine parameters used in the model of each scattering mechanism. Without fitting parameters, first-principles methods have been developed to predict the scattering rates of charge carriers in Si, 24 Si 1−x Ge x alloys, 25 GaAs and GaP. 26 These scattering rates can be further combined with the Boltzmann transport equation for the electrical property predictions. 27 Similar to phonon studies, such first-principles electrical property predictions are still restricted to well-studied materials and are thus unsuitable for high-throughput material search, in which thousands of complicated materials will be screened based on the computed Z T s. As the result, high-throughput TE material search is mostly based on S 2 σ/τ or the Z T s further estimated with guessed k L /τ values, [13][14][15][16] here τ is the averaged relaxation time of all charge carriers. Similar simplification is also used for first-principles Z T estimation of standard TE materials. [28][29][30] For materials screening, experimental inputs of σ and k L are often required and the guidance from computations becomes very limited.
For nanograined bulk materials, however, the need of computing electron and phonon MFPs can be eliminated when these MFPs are largely restricted by the grain size. In such materials, the optimum grain size a should satisfy l E,bulk < a < l P,bulk to reduce k L without deteriorating σ. When a < l P,bulk the Casimir limit is reached, 31 assuming 0.5 as grain-boundary phonon transmissivity based on the diffuse mismatch model. The intrinsic multi-phonon scattering within grains becomes negligible, the grain-boundary phonon scattering dominates and is the only factor to restrict the lattice thermal conductivity. The maximum Z T is anticipated when a is reduced to l E,bulk for majority charge carriers, called the small-grain-size (SGS) limit for nanograined bulk materials. 32 Within nanograins, phonons' point defects (alloy atoms etc.) scattering rates are proportional to the fourth power of the phonon frequency. Some high-frequency phonons can thus have MFPs below the SGS limit and further lower the lattice thermal conductivity. With more defects on grain boundaries, the thermal conductivity can be lower but it is mostly attributed to reduced phonon transmissivity. 33 With coated grain boundaries, 34 it is possible to further reduce the k L by lowering phonon transmissivity at the grain boundaries 35 and prevent grain growth at high temperatures. With further decreased a, both the electrical and thermal conductivities will be scaled down with a as classical size effects, but S will remain the same with unchanged carrier concentration, resulting in saturated Z T . This simple treatment has been used to compare the power factor S 2 σ of new TE materials. 11,12 However, the overall Z T is not evaluated without calculating k L .
In this work, we propose a Z T formulation of general materials under the SGS limit, which is solely based on phonon dispersions and electronic band structures predicted by first-principles calculations. Ideal crystal structures are used as input to the first-principles calculations. Defects are not considered, but will mostly reduce the thermal conductivity more than the electron mobility, leading to higher Z T values. In calibrations with nanograined bulk Si, this SGS estimation yields Z T close to the prediction by accurate models, which enables high-throughput material search for novel TE materials.

ZT Estimation under the SGS Limit
Model description.-Assuming a constant electron MFP l E = a, S 2 σ/l E or simply S 2 σ/a has been used for the power factor comparison in material search. 12 The expressions of S and σ are given by solving the Boltzmann transport equation: is the Fermi-Dirac distribution function and k B is the Boltzmann constant. The function σ( ) is related to the relaxation time τ i,k and group velocity v i,k for an electron state with band index i and momentum k within the electronic band structure of a bulk material. Because MFPs are limited by the grain size a, we have τ i,k = a/|v i,k | and [3] whereê is the unit vector of the transport direction and the integration is across the first Brillouin zone (BZ). In addition, the electronic thermal conductivity k E is given as [4] As the effective values within an isotropic nanograined bulk material, above electrical properties will be averaged over x, y, and z directions forê. The phonon dispersion for a bulk material can be written as ω i,q for a state with branch index i and momentum q. For an individual grain with size a, the phonon MFP within a grain is denoted as l G and it satisfies l −1 G = l −1 P,bulk + a −1 , and the phonon lifetime can be written as τ P i,q = l G /|v P i,q |, where the superscript P indicates quantities for phonon. The lattice thermal conductivity of each grain is denoted as k G and it can be calculated by averaging the kinetic relationship over all threeê directions, where N is total number of q, is volume of unit cell, and C is the volumetric phonon specific heat.
In nanograined materials, all grains are joined by grain boundaries with a thermal resistance R K . Based on an effective medium formulation, 33,36 the effective lattice thermal conductivity for nanograined material k eff is given as , [6] in which for a grain boundary. Equivalently, the effective phonon MFP of nanograined materials is given as [7] At the SGS limit when a l P,bulk , we have l G ≈ a, and l P,eff ≈ 3a/7. Eq. 6 simply gives k eff = 3k G /7 for more accurate Z T predictions.
Dividing σ, k E , and k eff with the common factor a, Z T of a material under the SGS limit can be computed as [8] which can be further optimized at a given temperature by tuning the Fermi level E F . As a more conservative Z T estimation, optical phonons are equally treated as acoustic phonons though some highfrequency optical phonons may still have MFPs much shorter than the grain size and can be overestimated for k eff . The employed electrical property calculations are carried out by modifying the package BOLTZTRAP, 16 where a constant electron MFP is assumed for charge carriers instead of a constant relaxation time in the original code. The full-band k eff calculation is similar to the electrical property calculations. Assigning a unity phonon group velocity in Eq. 5, the phonon code is calibrated with specific heat computations of bulk Si and ZnO, both of which agree within 2% with the experimental data from 50 to 300 K. The electronic band structures are given by the Vienna Ab-initio Simulation Package 37,38 (VASP) based on the density functional theory using GGA approximation. 39 The phonon dispersions are calculated using PHONOPY package. 40

Calibration with nanograined bulk silicon for its SGS limit.-
To evaluate the inaccuracy introduced by the SGS approximation, calibrations have been carried out on n-type nanograined bulk Si using energy-dependent electron and phonon MFPs given for bulk Si. 36,41 An electron concentration of 2 × 10 20 cm −3 is used in the calculations. Only longitudinal acoustic (LA) and transverse acoustic (TA) phonons are considered for k eff calculations because no detailed expression of optical phonon MFPs is given in previous publications. In Eq. 3 for electrons and Eq. 5 for phonons, the grain size a is now replaced by the effective MFP l eff modified from the bulk value l bulk , given as l eff = (l −1 bulk + a −1 ) −1 . All employed parameters and detailed methods are given in Appendices A and B, respectively for electron and phonon calculations for this calibration. Figure 1 shows the grain-size-dependent S 2 σ, k E , k eff , and Z T at 300 K and 1000 K. Below 10 nm, Z T would saturate around the value for the SGS limit (∼0.123 at 300 K, ∼1.08 at 1000 K). The SGS limit k E /a is higher at 1000 K due to thermally activated charge carriers. Despite different parameters used in the analysis, the trend of temperature-dependent TE properties is consistent with previous analysis for nanograined bulk Si. 6 In real samples, charge carriers are also scattered by the grain-boundary energy barrier that is formed as a result of charges trapped by a grain boundary and is associated with a depletion region. 42 In detailed analysis, the scattering rate of this potential field can be calculated 43 and then added to l E using the Matthiessen's rule, with the unknown energy barrier height as a fitting parameter. 6 To evaluate the additional contribution by optical phonons, k eff predicted by Eq. 6 is evaluated with and without considering the optical branches (Fig. 2). At 300 K, 24% of lattice thermal conductivity k L is contributed by the optical phonons under the SGS limit. In comparison, ∼20% k L contribution by optical phonons is suggested within a Si nanowire that has a 10-nm diameter to restrict all phonon MFPs. 44 In the nanowire calculation, l eff based on first-principlescomputed bulk phonon MFPs is used to replace a in Eq. 5. The 4% difference is anticipated because there are still some optical phonons with MFPs remarkably smaller than the 10-nm wire diameter and k L contribution of these phonons is overestimated under the SGS limit. However, such a divergence will not affect the Z T evaluation of a material.

Results and Discussion
Within the scope of this work, the focus is on novel oxides as high-temperature TE materials, which can recover high-quality waste heat from various resources such as industrial furnaces, airplane jet engines, and power plants. In addition to their high-temperature stability, all computed oxides possess a wide bandgap (> 2.1 eV) to prevent the detrimental bipolar conduction that is caused by minority charge carriers thermally excited at elevated temperatures. 1  minority charge carriers cancel out S mediated by majority charge carriers (lower S 2 σ) and simultaneously increase k E (and thus k), resulting in decreased Z T above a threshold temperature. In this aspect, a sufficiently large bandgap of oxides is crucial to maintaining a high Z T at high temperatures.
ZT of standard materials under SGS limit.-Two n-type (wurtzite ZnO, SnO 2 ) and two p-type oxides (ZrOS, Ca 4 P 2 O) as wide-bandgap electrode materials are first investigated at 1500 K. With a high melting point T m > 2000 K, these materials and their nanostructures can be thermally stable during long-term operations and thus reach beyond the state-of-the-art SiGe alloys that are stable only up to ∼1200 K. The selection of electrode materials for TE applications is somewhat anticipated because all existing TE materials tend to have a high roomtemperature σ (∼1 × 10 5 S/m), as suggested by the database of existing TE materials. 45 For device fabrication, such electrode materials can also form superior electrical contacts to minimize Joule heating on junctions to avoid energy loss. At 1500 K, the computed k E , σ/a, S, S 2 σ/a, and Z T are plotted as functions of E F in Figs. 3a-3e. If effective p-type doping can be achieved, Z T 1500K around 3 can be obtained in Ca 4 P 2 O.
The carrier concentration associated with the optimized E F is further used to predict their temperature-dependent Z T s (Fig. 3f) across the whole temperature range, which monotonically increases with elevated temperatures for wide-bandgap materials. More accurate analysis should further consider the dopant activation at varied temperatures, in which the impurity energy level of the selected dopant and its possible impact on the electronic band structure can also be predicted by first principles. 46 for p-type ZnO and SnO 2 , only their n-type Z T s are considered here due to the long-term challenge in their p-type doping. 48,49 Under the SGS limit, analysis based on a single parabolic band suggests that a large effective mass will always benefit S 2 σ and heavy holes are thus better than light electrons. 12 In addition, a high effective mass also leads to a lower σ and thus reduced k E = LσT to benefit the Z T , where the Lorenz number L is roughly 2.4 × 10 −8 W K −2 for heavily doped samples. 1 The reduced k E becomes more important to nanostructured materials with largely suppressed k L . Among these four materials, the slightly higher S 2 σ/a of n-type Ca 4 P 2 O and ZrOS are due to additional electron valleys close to the conduction band edge (Figs. 4a and 4b), which function as extra "electron pockets" to increase σ for the same E F . Such band degeneracy has been used to achieve a high power factor and thus Z T in PbTe 1−x Se x and Mg 2 Si 1−x Sn x alloys. 50,51 In nanograined Ca 4 P 2 O and ZrOS, S 2 σ/a of n-type samples benefits from these additional electron valleys but the maximum Z T s of p-type samples are still higher due to their lower k E (Fig. 3a).
Detailed information for optimized E F and other properties are summarized in Table I. In general, maximized Z T is obtained at |S| from 189 to 309 μV/K, which are higher than those to maximize S 2 σ (100-170 μV/K). The latter is comparable to the 130-187 μV/K range suggested for optimum S 2 σ of conventional materials. 52 Under the SGS limit, k E becomes comparable to k L so that the S optimization is an intermediate situation between two extreme cases, i.e. k E k L and k E k L . The S 2 σ optimization applies to the Z T optimization under the condition k E k L . As an opposite case, k E k L leads to Z T ≈ S 2 σT /k E = S 2 σT /LσT = S 2 /L and a larger |S| always benefits the Z T . In practice, the increased optimum |S| under the SGS limit indicates a lower doping level that is less challenging in experiments.
With four atoms per primitive unit cell, wurtzite ZnO is predicted to reach Z T 1300K = 0.23 under the SGS limit. This value can be treated as the lower Z T bound because the nine optical branches may be largely overestimated for their k eff contribution, particularly the upper six optical branches with significantly higher frequencies than the rest phonon branches (Fig. 5). When only three acoustic branches are considered for k eff , Z T 1300K increases to 0. 66 53 which further benefits from the electronic band structure variation of doped ZnO (Ref. 47) and strong k L reduction due to the alloy scattering of phonons. In particular, the transport of optical phonons are anticipated to be largely suppressed by such point-defect scattering. Instead of completely neglecting these optical branches, their contribution to k eff can also be estimated as the theoretical minimum for amorphous solids, in which the phonon MFPs become half of the phonon wavelength. [54][55][56] An even higher Z T is anticipated in nanograined bulk ZnO alloys, with alloy atoms to scatter short-wavelength phonons and grain boundaries to scatter mid-to long-wavelength phonons. In general, Eq. 8 can be used for fast-screening of novel materials and the obtained Z T can be improved by multi-length scale structures to suppress k L across the whole phonon spectrum.
ZT of general oxides under SGS limit.-As more general cases, other oxides selected from the Inorganic Crystal Structure Database (ICSD) are optimized for their Z T s. Without the exact melting points for some oxides, all Z T optimizations are thus carried out at 800 K that is generally safe for oxides. In particular, PbTiO 3 has ferroelectricparaelectric phase transitions at 763 K and experiences negative thermal expansion. 57 The predicted Z T at 800 K may improve due to possible lower thermal conductivity after phase transition. The obtained S 2 σ/a, S, k eff /a, and k E /a are shown in Figs. 6a and 6b. As the major p-type electrode materials, 58 AlCuO 2 can achieve Z T 800K > 2 under the SGS limit and similarly high Z T s are also predicted for pand n-type HfOS (Fig. 6c).
As general trends, the optimized Z T 800K values are plotted as a function of the bandgap E g (Fig. 7a) and atom number n within a primitive cell (Fig. 7b). Figure 7a shows that the optimized Z T s may still increase when E g k B T , here T is operation temperature. In bulk materials, E g of 6 − 10k B T is suggested for best TE materials with indirect bandgap, 59,60 while direct bandgap materials with special scattering mechanism of charge carriers may have E g > 10k B T for optimum Z T . 61  but the actual Z T for a wide-bandgap material often decreases significantly due to increased k L . This rule becomes invalid with the dramatic k L reduction by nanostructures and alloying atoms, as shown in the high-Z T nanostructured ZnO alloys 53,62 and GaN alloys. 63 Secondly, Z T monotonically increases with n (Fig. 7b). As suggested earlier, S 2 σ benefits from larger n. 12 For phonons, there are 3 acoustic and 3(n − 1) optical branches within a material. The fraction of optical phonons is increased for larger n and their contribution to k eff is usually negligible due to their small group velocities and MFPs. By neglecting optical branch contribution, k L ∼ n −2/3 was proposed for bulk materials by Slack. 64 Under the SGS limit, however, optical phonons become more important because their usually short MFPs are now equal to those for acoustic phonons (i.e. grain size) in Eq. 5. Similar argument can be found for Si nanowires when phonon MFPs are mostly restricted by the wire diameter. 44 As the lower and upper bounds of k eff , considering only acoustic branches leads to k eff ∼ n −1.3 at 800 K (dashed line in Fig. 7c), whereas considering all branches gives k eff ∼ n −0.4 (solid line). The k eff ∼ n −1.3 dependence is close to the k eff ∼ n −1.0 trend estimated for acoustic modes under strong boundary scattering. 56

Conclusions
In summary, a first-principles Z T evaluation of novel materials under the SGS limit is proposed in this work and demonstrated in representative oxides. This enables broad search of next-generation TE materials with low materials cost and environment beneficial. Unrestricted to oxides, the 173473 crystal structures within the ICSD can be re-evaluated for TE applications. The materials search can be further extended to those that do not exist in nature but can be thermodynamically stable based according to their first-principlescomputed formation energy. 13,14,65 A more general approach of nonexisting material predictions can be found for batteries, 66,67 which can be extended to TE materials in the future.

Acknowledgments
We thank Prof. Zi-Kui Liu for discussions on PbTiO 3 . An allocation of computer time from the UA Research Computing High Performance Computing (HPC) and High Throughput Computing (HTC) at the University of Arizona is gratefully acknowledged.

Appendix A: Detailed Charge Carrier Analysis
For the purpose of calibration of electronic properties of Si with the established calculation, we choose the Vining model 41 and adopt its original parameters. A general form of Si 1−x Ge x was considered by Vining. 41 For pure Si, the Ge fraction x is set as zero. A simple two-band analytical electronic energy dispersion is assumed. The effective mass for both electron and hole are set to m * = 1.40m e , where m e is the mass of electron.
The electrical conductivity for each band is given as where z is reduced charge carrier energy z = E/k B T , with temperature T and Boltzmann The subscript i can be either e indicating conduction band, or h for valence band. The total electrical conductivity is simply σ = σ e + σ h . The energy dependent charge carrier relaxation time τ(z) is discussed at the end of this Appendix. The Seebeck coefficient for each band is given as with '−' for conduction band and '+' for valence band. The overall Seebeck coefficient is obtained by weighting each band's contribution by its normalized electrical conductivity The electrical thermal conductivity k E is calculated using Wiedemann-Franz law by k E = L 0 σT , where the Lorenz number L 0 is the sum of three terms: the contribution from conduction band L e , from valence band L h , and bipolar contribution L b . The contribution from each band is given as The bipolar Lorenz number is given as For our nanograined Si, we consider three charge carrier scattering mechanism, so that according to Matthiessen's rule we have where E is energy of charge carriers, τ AC , τ I , and τ N G are relaxation time of scattering with intravalley acoustic phonon, ionized impurity, and nanograins, respectively. The last scattering is what differs from the Vining model. 41 Scattering with acoustic phonon is described with deformation potentials. Following Vining, deformation potential for both electrons and holes are identical as E d = 2.94 eV. 41 We then have where ρ is density, v s is sound velocity, and electronic density of states is D(E) = The screening length r s of coulomb interaction is given as and relaxation time is given as ) unless CC License in place (see abstract where N D is number density of ionized impurity, and ξ = (2kr s ) 2 .
Finally, scattering with nanograins can be calculated using the size of nanograin a and the parabolic band [A9]

Appendix B: Detailed Phonon Analysis
To simplify calculation, isotropic phonon dispersion is assumed here and the phonon dispersion along the (001) direction for bulk silicon 68 is employed. In this case, Eq. 6 is simplified as an integration with respect to ω and summation over branch i: .
Only acoustic phonons are considered here without an explicit expression of energydependent optical phonon MFPs. The phonon MFP l P,i (ω) is modified from the bulk l Bulk P,i (ω) as 1/l P,i (ω) = 1/a + 1/l Bulk P,i (ω). For l Bulk P,i (ω) computations, the bulk phonon lifetime, τ Bulk P,i (ω) = l Bulk P,i (ω)/v P i (ω), is given by the Matthiessen's rule: Here τ I (ω), τ N ,i (ω), τ U,i (ω), and τ E (ω) are the phonon lifetime associated with impurity scattering, normal (N) process, Umklapp (U) process, and electron scattering, respectively. Based on first-principles calculations on bulk Si, the relaxation times of the momentum-conserved N process and momentum-non-conserved U process are calculated as in Refs. 18 and 69, in which T is the absolute temperature,h is the Planck constant divided by 2π, D is the Debye temperature, A N /U,i is material-dependent coefficient further modified for isotropic phonon dispersion of Si (Table BI). 36 For natural Si with isotopes, the impurity-phonon scattering τ I (ω) is given as in Ref. 70, where V 0 is the averaged volume per atom (2.0 × 10 −29 m 3 for Si), g is the mass variance determined by the composition, D(ω) is the phonon density of states at ω. Here g = j f j (1 − m j /m) 2 , with m j as the atomic mass for the j-th isotope, f j the molar percentage of individual isotope atoms, and m = j f j m j the averaged atomic mass. In heavily doped n-type silicon (n > 1.0 × 10 18 cm −3 ), shallow impurity levels within the bandgap start to merge with the conduction band so that the dopants are always completely ionized. 71 For a carrier concentration of 2 × 10 20 cm −3 , there is 0.4 mol% phosphorus atoms in Si and g is estimated to be 2.41 × 10 −4 , which is slightly higher than g = 2.01 × 10 −4 for natural Si. 72 For heavily doped TE materials, the electron scattering of phonons plays an important role in its k L reduction and τ E (ω) is expressed as 163921 1170 507 645 in which k B is the Boltzmann constant, E d = 2.938 eV is the acoustic deformation potential, m * = 1.40m e is the density-of-states effective mass, m e is the free electron mass, d = 2327 kg/m 3 is the density, v g = 5880 m/s is the averaged phonon group velocity, 41 and E F (0.082 eV at 300 K, 0.0036 eV at 1000 K) is the computed Fermi energy referred to the conduction band edge.