Assessing the Thermoelectric Properties of Sintered Compounds via High-Throughput Ab-Initio Calculations

Shidong Wang, Zhao Wang, Wahyu Setyawan, Natalio Mingo, and Stefano Curtarolo* Department of Mechanical Engineering and Materials Science, Duke University, Durham, North Carolina 27708, USA LITEN, CEA-Grenoble, 17 rue des Martyrs, 38054 Grenoble Cedex 9, France Department of Mechanical Engineering and Materials Science and Department of Physics, Duke University, Durham, North Carolina 27708, USA (Received 6 July 2011; published 29 November 2011)


I. INTRODUCTION
Thermoelectric materials, converting heat into electric power, are one of the subjects of growing interest in the field of energy conversion [1][2][3].Currently, energy harvesting through thermoelectrics is limited by the low efficiencies of present commercial solutions [2].To tackle the problem, several nanostructured materials have been proposed since the late 1990s [4].Among them, the sintered powder nanocomposites stand out as cost-effective solutions with improved thermoelectric efficiency [5][6][7][8][9].In these materials, the presence of grain boundaries alters the electron-and phonon-scattering mechanisms.In particular, the shorter phonon mean free path leads to a thermal conductivity, , smaller than the bulk value.A low is a required ingredient in the recipe for a good thermoelectric system.The conversion capability of a system is captured by the dimensionless figure of merit, ZT ¼ P Á T=.Here, T is the temperature and P ¼ S 2 is the power factor, with being the electrical conductivity and S the Seebeck coefficient.In addition to a small , a large P is another ingredient for efficient materials.While experimental trial-and-error approaches are expensive and time-consuming, theoretical investigations, especially within the high-throughput approach [10,11], could be both cost-and time-effective in pinpointing good candidates.
There are two different approaches for assessing thermoelectric properties from first principles.In the lowthroughput approach, the properties are addressed as accurately as possible by using a minimum possible number of assumptions and adjustable parameters.A few articles have been published in this direction, which solve the Boltzmann transport equation [12], and aim at computing all the scattering rates ab initio (e.g., Ref. [13] for SiGe, and Refs.[14,15] or [16] for mobility or of Si, respectively).The approach is accurate, but it is limited to perfect crystals and simple alloys: it is computationally very expensive and usually difficult to transfer to other systems.In the high-throughput (HT) approach, rather than high accuracy, the goal is to compute as many compounds as possible using highly automated techniques.Results are then saved in large material databases for further processing (e.g., the databases in Ref. [17]).The HT approach has been used to find new materials and trends with data-mining algorithms [10,11,[18][19][20][21].As an example, Madsen reported an ab-initio search for thermoelectrics containing Zn and Sb [22].Both approaches have advantages: the HT method generates more candidates than the low-throughput one, but the tradeoff between accuracy and speed/transferability results in occasional outliers in the HT method, which have to be corrected manually.
The constant-relaxation-time approximation, very common in ab-initio transport calculations [22,23], is hardly justifiable within HT frameworks, since carriers in different materials have different relaxation times, r .Therefore, many works have restricted the search to the value of the Seebeck coefficient only, which does not depend on r within the constant-relaxation-time approximation.The Seebeck coefficient alone, however, is not enough for comparing thermoelectrics.There is a necessary practical rule related to S: To have ZT > 1 at room temperature, one needs S > 156 V=K.(This can be easily derived from the Wiedemann-Franz law.)Although not sufficient, the threshold can be used to reduce the size of the candidate pool.
For sintered powders, the situation is radically different.For these systems, a simple and physically sound model adequate for HT calculations can be introduced.Grain boundaries strongly affect carrier scattering.It is reasonable to assume that, for rough interfaces, the scattering becomes diffusive: Carriers lose the memory about their path, and they get scattered in all directions with equal probability.This picture, forming the basis of the ''diffusemismatch model'' of phonon transport, has been employed with success for understanding the thermal conductance of solid-solid interfaces [24][25][26][27].Thus, it is reasonable to assume a similar diffusive-scattering ansatz for the limiting effect of grain boundaries on charge transport: The carrier mean free path (MFP), , becomes of the order of the grain size, L. For grains smaller than the MFP of the bulk, one can approximate ' L. Contrary to the constantrelaxation-time approximation, this assumption, called the constant-mean-free-path approximation, is based on physical grounds.The value of is not fitted from experimental data, but it is concretely linked to a well-defined quantity, the grain size.This allows for fair comparisons between power factors so that good candidates for sintering are possibly selected.
Representing with bulk e the average electron MFP in bulk single crystals, we define the small-grain-size (SGS) limit of a material as the regime where L < bulk e and where the power factor P becomes proportional to the grain size, P ¼ A Á L, where A is a material's specific constant [9].Outside of this regime (i.e., for large grain sizes), the power factor tends to its bulk value.Often, the grain-size reduction affects phonon transport before it changes electron transport.This is because bulk phonon mean free paths bulk ph are often larger than bulk electron ones bulk e .Then, in the cases where bulk e < L < bulk ph , the lattice thermal conductivity will be proportional to the grain size L, and ZT is expected to increase as 1=L.When L < bulk e , the power factor becomes proportional to the grain size, and therefore ZT saturates to a maximum value when the SGS limit is reached.This approximate argument suggests that the best ZT of a sintered material will be achieved within the SGS limit.Assessing the relative small-grain-size capabilities of different materials tells us whether their power factors are more likely to be higher or lower relative to each other in the SGS limit (i.e., when ZT is maximized).The characterization does not conclude which materials have the highest absolute ZT, as this would require knowing the lattice thermal conductivity.Nevertheless, a high power factor is a crucial property for a thermoelectric material, since the maximum power that can be generated by a device is directly related to its power factor P dev as W ¼ P dev ÁT 2 =4L dev [1].
In this article, we develop the model for sintered powders, focusing on the ordering of the power factors within the SGS limit.We present the results for more than 2500 compounds extracted from the online repository of electronic-structure calculations AFLOWLIB (the first database in Ref. [17]), based on compounds from the Inorganic Crystal Structure Database (ICSD) [28,29].We then introduce regression analysis to reveal correlations between the power factor and different physical quantities such as carrier mass and band structures.

II. COMPUTATIONAL SCHEME FOR THE POWER FACTOR OF SMALL-GRAINED COMPOUNDS
A. Transport-property calculations The power factor, P ¼ S 2 , is calculated from the definitions of the electrical conductivity and the Seebeck coefficient, S: where f 0 1=½e ðÀÞ=k B T À 1 is the Fermi-Dirac distribution function, is the chemical potential, T is the temperature, is the mean free path, and AEðÞ is an energy-dependent function defined below.As usual, we find the maximum power factor by adjusting the chemical potentials for both n-and p-doped compounds [14].
Especially at high concentrations, dopants can significantly alter the electronic structures of the host compounds and might invalidate the rigid-band approximation (RBA) used in our calculations.It has been shown that the RBA fails when the Fermi level is $0:1 eV into the band [30].Certainly, more accurate results of P could be obtained by characterizing the whole electronic structure at different doping levels.However, due to the size of the materials library, this task is not computationally feasible within the current HT approach.Furthermore, the optimized power factors typically occur for the Fermi level about 0.02 eV into the band, in a regime in which the RBA is expected to be applicable [14,15,31,32].Hence, we consider the RBA to be a reasonable compromise between efficiency and accuracy.
The function AEðÞ is related to the scattering and the electronic structure.For small-grained compounds, boundary scattering dominates.As discussed before, we assume a constant-mean-free path equal to the average grain size.The lifetime for a given state of momentum k, band index , and velocity v ;k is ;k ¼ =jv ;k j.The function AEðÞ is defined as Here, v z is the electron group velocity in the transport direction, and the integral is carried out over the first Brillouin zone (BZ).We have performed transport calculations in all of the three principal axes.To obtain the power factor of the sintered powders, one needs to consider the random distribution of the grains.Although effectivemedium theories have been developed to model random systems, the calculation requires the direction-dependent thermal conductivities of the grains [33].Without this information, the safest guess is that the power factor will fall between the minimum and maximum values calculated along the principal axes.(Values are listed in the supplemental material; see Ref. [34].)To integrate Eq. ( 3), we use a uniform grid in the first BZ and calculate energies and velocities at the various points: f ;n ; v ;n g.We choose a mesh of at least 10 000 points, obtained by interpolating the HT ab-initio values and by using the modified-quadratic-Shepard method [35].The function AEðÞ reduces to the summation where g a ðxÞ ¼ e Àjx=aj 2 =a ffiffiffiffi p is a Gaussian function, N grid is the number of grid points, and V 0 is the volume per atom.An adjustable Gaussian smearing parameter a is used to control the smoothness of AEðÞ.This parameter, related to the energy spacing of the underlying grids, is optimized to obtain acceptable integration smoothness without losing information from the electronic-structure data: a is chosen as the ''k-points-average'' energy difference between the eigen-energies close to the gap edges (within 1 eV).
Density functional theory is used with augmented-wave pseudopotentials [42,43], and with the generalized gradient approximation for exchange-correlation functions (Perdew-Burke-Ernzerhof [44]).Calculations are performed without zero-point motion and at zero temperature.All structures are fully relaxed several times.Convergence tolerance of $1 meV=atom is expected using dense grids of 3000-4000 k points per reciprocal atom (KPPRA).At the beginning of the structural relaxation, a spin-polarized calculation is performed for all structures.If the magnetization is smaller than 0:025 B =atom, the spin degree of freedom is further neglected to enhance speed.A much denser grid with 10 000 KPPRA is used to produce accurate charge densities and density of states.In all the grid generations, the Monkhorst-Pack scheme [45,46] is used except for face-centered cubic, hexagonal, and rhombohedral systems, in which a À-centered mesh is automatically adopted for faster convergence.In addition, the effective density of k points in each reciprocal direction is chosen self-consistently so that it forms a Bravais lattice compatible with the reciprocal lattice.After this step, a mesh of at least 10 000 true k points is interpolated [35] for integrating the power factor [47]. On-site Coulomb interactions of localized electrons are corrected with the DFT þ U technique [48] for 3d-orbitals in transition metals and 4f-orbitals in lanthanides.Dudarev's [49] implementation of DFT þ U is used.It relies on U eff U À J, the difference between the Hubbard-like Coulomb interaction U and the screened Stoner exchange parameter J.The used values of U eff are listed in Table I and in Ref. [10].For all the possible Brillouin-zone shapes, the band structures' integration paths are taken from the AFLOW standard [11].

III. ANALYSIS AND DISCUSSION OF THE RESULTS
The 20 n-and p-doped powders with the highest SGS power factor normalized to grain size (P=L) are shown in Fig. 1.(The properties of all compounds included in this study are listed in the supplemental material; see Ref. [34].) P is direction dependent.Panels 1(a) and 1(b) show the best 20 systems with respect to the simple direction-averaged value, hPi=L.In polycrystalline samples, this average does not necessarily correspond to the actual P=L; therefore, in panels 1(c)-1(f), we show the compounds with the highest maximum [1(c) and 1(d)] and highest minimum [1(e) and 1(f)] P i =L values along the principal axes i.
To estimate the magnitude of L for each material in the SGS limit, one needs to know the value of the mean free path which can be extracted from the experimental mobility.Unfortunately, the mobility does not appear to have been reported for a large number of the compounds in the list.Usually, nonionic systems have longer carrier MFPs than highly ionic compounds.For the latter, the classical Â is typically several hundred K. (For example, Â $ 670 K for TiO 2 [51].)This yields the shortest MFP in the order of 1 nm.In the case of titanium oxide, for example, the bulk MFP is indeed close to 1 nm, as reported by Hendry et al. [52].In contrast, for ZnO, the use of the simple model in Ref. [53] and the largest measured bulk mobility available in the literature of 200 cm 2 =Vs [54] yields a mean free path of 7 nm.
The thermoelectric properties of micrograined polycrystalline titanium oxides have been measured, giving P $ 2 W=cmK 2 [55].Using the value of bulk e estimated above, our calculations for TiO 2 yield P values in the same  T i   order of magnitude.However, the Magne ´li phase Ti 5 O 9 stands out (Ti 10 O 18 , Pearson symbol aP28, ICSD #31399).
Even with a 1 nm grain size, the P values obtained for Ti 10 O 18 (and also for B 2 Cs 3 Li 2 NaO 6 ) are surprisingly large.For example, the largest experimentally reported power factor corresponds to FeSb 2 , with P$2300 W=cmK 2 at a temperature of 12 K [56], and measured values for strongly correlated systems YbAl 3 (P $ 340 W=cmK 2 at 80 K [57]) or YbAgCu 4 (P $ 235 W=cmK 2 at 20 K [58][59][60]) are regarded as exceptionally high.In addition to having larger effective masses, Ti 10 O 18 has P larger than other titanium-oxide compositions by more than an order of magnitude.Even assuming that the least-optimal directions dominated the overall P of the sintered system, this value is surprising.Since the experimental samples were not doped, the measured P do not correspond to the doped optimized values, which might explain the difference.Doping optimization may not always be simple, especially for large-band-gaps systems.Too-large dopant concentrations may also reduce mobility and affect the electronic structure of the compound.All these precautions should be considered when comparing our HT results to experiments.
The following considerations are relevant.First, for HT calculations, there is no expectation of accuracy beyond the estimate of the order of magnitude.This is sufficient to establish trends and correlations, as the P values span 6 orders of magnitude throughout the list.Also, the properties of a sintered powder critically depend on the fabrication process.Different techniques, such as spark-plasma sintering or hot-isostatic pressing, result in different-quality grain boundaries, which also depend on the particular chemical compounds.Therefore, a high position in the list indicates only a good potential towards high power factor, but it does not necessarily predict a final good thermoelectric material.Factors such as charge trapping at the interfaces may overshadow the good electronic-structure properties captured by the HT-model calculations.
Directions for searching good sintered-powder thermoelectrics are extracted from correlations between P and other properties.Between the variables fx; yg, the Pearson correlation [61] is defined as x;y covðx; yÞ= x y , where covðx; yÞ is the covariance between the variables, and x and y are the standard deviations, respectively.
In Table II, we show the positive correlation between P=L and the band gap, the charge-carrier effective masses, and the ratios between the maximum and minimum chargecarrier effective masses.P=L is more weakly correlated to the hole effective mass (m h ) than to the electron effective mass (m e ), and it is correlated only to the ratio of maximum and minimum electron effective masses.For p-doping, it is not correlated with the ratio of the masses.The reason behind this asymmetric behavior is unclear and currently under investigation.
We find also that compounds with a large number of atoms per primitive cell tend to have the highest P.This trend is depicted in Fig. 2 where the median of P=L is shown to increase with increasing atoms per primitive cell.At room temperature, the state-of-the-art Bi 2 Te 3 compound has been reported to yield power factor as high as P ?¼ 66:1 W=cmK 2 [62].Although most compounds in our list have P lower than this value when L is a few nm, there are still some with P comparable to Bi 2 Te 3 .These are good candidates for sintering.
The correlation found between the power factor and the effective mass can be qualitatively understood by considering a simple model with a parabolic dispersion.Following Mahan and Sofo [63], we can rewrite the thermoelectric properties of this model as with 0 ¼ e 2 =@a 0 , Bohr's radius a 0 , and the integrals Here sðxÞ ¼ @a 0 AEð þ xk B TÞ and the function AEðEÞ is given in Eq. (3).For parabolic dispersion, and under the assumption of a constant-mean-free path , the functions sðxÞ and AEðEÞ reduce to with carrier effective mass m Ã .Then the power factor becomes and its maximum is given by where m 0 is the electron mass and At room temperature, for ¼ 5 nm, the maximum power factor is P max ' 8ð m Ã m e Þ W=cm K 2 , indicating the relationship between P and the carrier effective mass.
The correlations between large power factors and large band gaps or effective masses are phenomenological rules that have been previously understood in the terms of simple models like the one above [64][65][66][67].The approach in this paper is radically different.No a priori assumptions have been made on the electronic structures, and indeed correlations between different properties are found by mathematical algorithms.The validation of the known empirical rules by the high-throughput approach corroborates the validity of the rules.This is in spite of the complexities of the actual band structures, which are much more intricate than the simple models employed to justify these correlations [64][65][66][67].In turn, this also cross-validates the power of the HT approach in discovering correlations: If the relationships P , bandgap and P , effective mass had not been rediscovered, the whole HT algorithmic apparatus might be questioned.
In addition to positively correlating with high P, a large primitive-cell size also benefits lattice thermal conductivity.Generally, the latter is expected to decrease when the number of symmetrically inequivalent atoms in the primitive cell increases.This is mainly due to the fact that the opticalphonon branches usually carry little heat.For a material with n atoms per primitive cell, there are 3 acoustic and 3ðn À 1Þ optical branches.The larger the n, the smaller the fraction of phonons that can contribute efficiently to the thermal conductivity, .This argument was put to practice in an insightful analysis by G. Slack [68].He showed that the actual dependence of on n is $ n À2=3 .Then high ZT may be expected in complex primitive-cell compounds, because of the two simultaneous effects: smaller and larger P. Whereas the thermal conductivity has a À2=3 exponent, we find that the power factor displays a linear dependence, P $ n, roughly.Therefore the power-factor effect may become relatively more important than the thermal-conductivity effect when n becomes very large.HT investigations of could describe its actual dependence on n, and its correlations with other properties, beyond Slack's simple approach.We plan to undertake this computationally demanding investigation in the future.
Note that we recently became aware of the work of Parker and Singh about the potential thermoelectric performance of hole-doped Bi 2 Se 3 [69].This material is included in our database.(Both Bi 2 Se 3 and Bi 2 Te 3 have been calculated with spin-orbit coupling.)Within our nanopowder framework, p-doped Bi 2 Te 3 has about $2:5 times larger averaged power factor hP p i=L than p-type Bi 2 Se 3 .(See page 16 of the supplemental material in Ref. [34].)

IV. CONCLUSIONS
By using HT ab-initio calculations, we have calculated the thermoelectric properties of more than 2500 sintered compounds extracted from the online repository of electronic-structure calculations AFLOWLIB (the first database in Ref. [17]), based on compounds from the ICSD.We have assumed (i) diffusive electron scattering at grain boundaries, (ii) grain sizes smaller than the mean free path of the bulk compounds, and (iii) constant-mean-free path equal to the grain size of the compound.Taking the same grain size for all compounds, we have calculated all the power factors and identified the best candidates for sintering.Guiding rules for searching for better materials are drawn by the correlations between the power factor and other physical properties.It is found that sintered thermoelectric compounds with expected large power factors tend to have large band gaps, heavy carrier effective masses, and many atoms per primitive cell.

FIG. 1 .
FIG. 1. Materials selection.The best 20 compounds for (a) n-and (b) p-doped average normalized power factor hPi=L by grain size L; (c) n-and (d) p-doped principal-axis maximum P i =L; and (e) n-and (f) p-doped principal-axis minimum P i =L.The corresponding Seebeck coefficients are also shown.

TABLE I .
[10]es of U eff parameters (in eV) for the Dudarev's GGA þ U approach implemented in AFLOW.Data taken from Ref.[10]and references therein.

TABLE II .
Calculated Pearson correlations between the normalized power factor and other properties: band gap (E g ), carriers' effective mass (m e , m h ), and ratios between the maximum and minimum effective masses.Power factor is positively correlated to most of these properties.
FIG. 2. Median of normalized power factors by grain sizeversus number of atoms per primitive cell for n-doped compounds (red squares) and p-doped compounds (green dots).High P correlates with the number of atoms per cell.