Limitations of zT as a Figure of Merit for Nanostructured Thermoelectric Materials

Thermoelectric properties of nanocomposites are numerically studied as a function of average grain size or nanoparticle density by simulating the measurements as they would be done experimentally. In accordance with previous theoretical and experimental results, we find that the Seebeck coefficient, power factor and figure of merit, zT, can be increased by nanostructuring when energy barriers exist around the grain boundaries or embedded nanoparticles. When we simulate the performance of a thermoelectric cooler with the same material, however, we find that the maximum temperature difference, max T ∆ , is much less than expected from the given zT. This occurs because the measurements are done in a way that minimizes Joule heating, but the Joule heating that occurs in operating devices has a large effect for these kinds of materials. The same nanocomposite but without energy barriers at the grain boundaries has a lower measured zT but a higher max T ∆ . The physical reason for these results is explained. The results illustrate the limitations of zT as a figure of merit for nanocomposites with electrically active grain boundaries.


I. INTRODUCTION
To address the challenge of increasing the thermoelectric materials' figure of merit, zT, Hicks and Dresselhaus suggested in 1993 that zT might be enhanced in nanostructures. 1,2 Following this strategy, researchers have been able to substantially and steadily increase zT. 3 Today, there are many reports of zT . 1, and all of them make use of nanostructuring of one kind or another. 4 This progress has primarily been achieved by nanostructuring materials to reduce the lattice thermal conductivity without substantially degrading the electrical conductivity, [5][6][7][8] but benefits to the electronic performance have also been reported. These benefits include an increased Seebeck coefficient (Refs. 9-21) and/or a reduction of bipolar effects. 6 Enhanced Seebeck coefficients are thought to be due to energy filtering near nanoscale inclusions or grain boundaries; 22 theoretical treatments 23 and models 24 have been reported. If nanostructured materials could be engineered to both lower the lattice thermal conductivity and improve the electronic performance, then substantial additional increases in zT would be possible. 4 The physics of interfaces due to nanoscale inclusions or grain boundaries has been studied, and various models have been proposed. In Ref. 25, the authors modeled the grain boundary as a parabolic band material with potential barriers for carriers and found that the model could readily explain experimental results. More detailed models have also been reported. [26][27][28][29][30] In Refs. 31 and 32, the authors showed that recent experimental results in rather complicated thermoelectric materials could be explained using simple models without considering the energy-dependence of the carrier scattering. While these models differ in detail, they all show that band bending at interfaces is the origin of differences observed in the Seebeck coefficients of bulk crystals and nanocomposites. Based on these prior studies, we model the interfaces as electrostatic potential barriers using full, numerical simulations of random structures. The details are discussed in Sec. II.
Two different types of nanocomposite thermoelectric materials are examined. The first type is a p-type polycrystalline Bi 2 Te 3 with grain boundaries that produce energy barriers. The second type is a p-type crystalline Bi 2 Te 3 with embedded silver nanoparticles that produce energy barriers around the nanoparticles. We use these as model systems for which experimental data are available, 17,33 but the conclusions should be more general. The thermoelectric properties as a function of average grain size or nanoparticle density are examined by numerically simulating the measurements of electrical and thermal conductivity and Seebeck coefficient as they would be done experimentally. In addition, we simulate the performance of an ideal thermoelectric cooler (i.e., no interface resistances or shunt conductances) and compare the zT extracted from ΔT max to the zT obtained from the simulated measurements of the thermoelectric parameters.
We focus in this paper on polycrystalline nanocomposites, but similar results (discussed in the supplementary material) are observed for nanocomposites with embedded Ag nanoparticles. For the case where the grain boundary or nanoparticle acts to deplete the grain, energy barriers for majority carriers result, and we find in accordance with experiments that the Seebeck coefficient, power factor (PF), and figure of merit, zT, can be increased by nanostructuring, but when we simulate the performance of a thermoelectric cooler, we find that the maximum temperature difference, ΔT max , is much less that what would be expected from the given zT. An important finding is that conventional measurements of the thermoelectric parameters of nanostructured thermoelectric materials may provide overly optimistic predictions of device performance. This occurs when the grain boundaries are electrically active because Joule heating, which is negligible in the zT measurement, plays a strong role under device operating conditions. The same nanocomposite but with no potential barriers at the grains has a lower zT but produces a higher ΔT max in a TE cooler. The physical reason for these results will be explained.
The paper is organized as follows. Section II describes the numerical techniques for the two-dimensional electrothermal simulations. Also described in Sec. II are the material parameters used in the simulations and the treatment of grain boundaries. Section III presents the results in terms of zT vs grain size where zT is extracted in three ways: (1) by simulating the measurement of the individual thermoelectric parameters, (2) by simulating Harman method measurements, 34 and (3) by deducing zT from ΔT max of a simulated cooler. We discuss the results in Sec. IV and summarize the conclusions in Sec. V.

II. APPROACH
In this study, we use the semiconductor device simulation program, Sentaurus, which solves the coupled partial differential equations that describe electrothermal transport in semiconductors. 35 The following equations 36 were numerically solved: is Poisson's equation, whereD is the displacement vector and ρ is the space charge density, which includes contributions from mobile electrons and holes and ionized dopants. Equations (1b) and (1c) are the electron and hole continuity equations whereJ n andJ p are the current densities that include contributions from gradients in the electrochemical potential and temperature (q is the magnitude of the elementary charge). The term R is the carrier recombination rate, which is only important in the presence of bipolar effects. The final equation (1d) is an energy balance equation for the temperature, T. Local thermodynamic equilibrium is assumed so that the electron, hole, and lattice temperatures are identical. The total heat capacity, c tot , and the total thermal conductivity, κ tot , include contributions from the lattice as well as the electrons and holes. Finally, the heat source term, H, includes contributions from Joule heat, recombination heat, and Thompson heat. 36 Transient simulations are performed only to simulate the Harman method-all other simulations are steady-state.

A. Material-level simulation and material parameters
In this work, Bi 2 Te 3 was chosen as a model thermoelectric material, because it has a simple rhombohedral crystal structure, its material properties are relatively well-known, 37 and experimental results are available. 17,34 The methods and results obtained in this work should, however, be applicable to other thermoelectric materials as well. The overall simulation workflow is summarized in Fig. 1.
Because thermoelectric materials have highly nonparabolic and complex band structures, the first step was to use density functional theory (DFT) calculations with Quantum Espresso 38 to obtain the Bi 2 Te 3 band structure. For both Bi and Te, we used full relativistic PAW pseudopotentials with a 40 Ry plane wave energy cutoff and spin-orbit coupling. An 8 × 8 × 8 and 20 × 20 × 20 Monkhorst-Pack k-point mesh was used for self-consistent and nonself-consistent calculations, respectively.
The second step was to process the full band structure with LanTraP, a program that extracts thermoelectric transport parameters from a given band structure. 9 We assumed that the scattering rate for carriers was proportional to the density-of-states and determined the electron-phonon coupling parameter by matching the resistivity to experimental values. 39 Transport parameters were extracted along the cross-plane direction. We did not include the transport anisotropy of Bi 2 Te 3 because our focus was on the effects of grain boundaries, and the random structures should average out anisotropies.
The simulated electrical resistivity and Seebeck coefficient vs hole density of p-type, crystalline Bi 2 Te 3 are shown in Fig. 2. These parameters are similar to reported values. [40][41][42] The maximum power factor (PF) occurs at a hole concentration of 4 Â 10 19 cm À3 . For the subsequent simulations of polycrystalline materials, each grain was doped at this optimal point to ensure that any performance increase observed in the polycrystalline material was due to grain boundary effects.

B. Electrical simulation of nanocomposite materials
To simulate nanocomposites, a random distribution of grains was first generated, and then the grain boundaries were defined. The algorithm used to generate this random geometry is described in the supplementary material. The grain boundaries were 100 nm thick. A simple thermal interface resistance was added in the middle of each grain boundary region to account for the additional thermal resistance. The value of this thermal interface resistance was set to agree with experimentally reported trends in thermal resistivity vs grain density. 17,43 The doping density of the grain boundaries was set at N GB A ¼ 3 Â 10 18 cm À3 to produce results consistent with experiments. The resulting band bending is about one-half of the bandgap. As shown in Fig. 2(b), the result was that the Seebeck coefficient of the grain boundaries was substantially higher than that of the grains.
The generated polycrystalline geometry was used as an input to the Sentaurus Device simulator. It was critical to resolve the band bending within and near the grain boundaries. For polycrystalline materials, a fine mesh on both sides of the grain boundary was needed. The added burden of a dense mesh is the main reason why we used 2D simulations instead of 3D. A full 3D simulation of a geometry that is large enough to be statistically meaningful would have taken days. As discussed in Sec. IV, simulations of 1D nanocomposites show results that are very similar to those obtained by 2D simulations; it is likely, therefore, that the effects found in the 2D studies would not change in 3D. (See the supplementary material for a description of how the 1D simulations were done.) The Sentaurus simulations used tables of transport parameters vs carrier density for all thermoelectric parameters; these tables were generated using LanTraP with the aforementioned DFT calculations. Phenomena such as the Thomson effect and mobility variations with doping density were, therefore, automatically included. For each random sample, we computed the conductivity by applying a small voltage (0.1 mV), computing the resulting current, and converting the result to a conductivity for the sample. To compute the Seebeck coefficient and thermal conductivity, we applied a  9 Lower right: Grain structure and numerical grid for one randomly generated sample. The dark areas show the high density numerical grid. The inset shows an expanded view of how the grid is refined near and inside the grain boundaries. The average grain size here is 1 μm, and the simulation domain is 1 × 1 mm 2 . Lower left: Simulated performance of a TE cooler. For a typical case, 100 random samples were generated and simulated, and the results were averaged.

FIG. 2. (a) Electrical resistivity and (b)
Seebeck coefficient vs hole concentration for the Bi 2 Te 3 baseline material used in this work. The bulk is doped at the optimal PF point, N G A ¼ 4 Â 10 19 cm À3 . To produce depleted grain boundary regions, N GB A ¼ 3 Â 10 18 cm À3 was used. For these calculations, the orientation of Bi 2 Te 3 was along c-axis.
small temperature difference (1 K), and simulated the resulting open-circuit voltage, which yields the Seebeck coefficient and heat flux, from which we determine the thermal conductivity. From the three parameters, we could compute a material figure of merit, zT, for the sample. We also determined zT by simulating a Harman measurement 34 for each sample. Finally, we extracted zT by simulating ΔT max of an ideal 1-leg TE cooler. The results of 100 such simulations for each average grain size were averaged and plotted.

III. RESULTS
The simulated Seebeck coefficient for polycrystalline Bi 2 Te 3 as a function of average grain size is shown in Fig. 3(a). Because the simulation domain (1 × 1 mm 2 ) is large enough to contain thousands of grains, significant grain size averaging occurs for each randomly generated sample. Still, for each average grain size, 100 randomly generated samples were simulated and averaged. The error bars shown in Fig. 3(a) show that the sample-to-sample variation is small. As intuitively expected, with decreasing average grain size, the Seebeck coefficient increases because the grain boundary regions have higher Seebeck coefficients, and as the average grain size decreases, there are more grain boundaries. Still, the size of the grain boundaries is small, so the grain boundary regions occupy a small fraction of the polycrystalline sample. If we repeat the simulation but assume no reduction in lattice thermal conductivity, then as Fig. 3(a) shows, there is a much smaller dependence of the Seebeck coefficient on grain size.
As noted above, the Seebeck coefficient of the sample increases substantially when the Seebeck coefficient of the grain boundary is high and the thermal conductivity of the grain boundary is low. 27,33 This can be understood with a simple 1D model, 33 where S G (S GB ) is the Seebeck coefficient of the grain (grain boundary), t GB is the thickness of the grain boundary, and κ G (κ GB ) is the thermal conductivity of the grain (grain boundary). Figure 2 shows that S GB % 2S G , but the grain boundary is thin, that it has only a small effect on the overall Seebeck coefficient unless the thermal conductivity of the grain boundary is also very low. This expectation is confirmed by the results shown in Fig. 3(a)-when the thermal conductivity of the grain boundaries is identical to the thermal conductivity of the grains, the Seebeck coefficients of the polycrystalline samples are close to that of the bulk, crystalline reference. Figure 3(b) shows the electrical resistivity vs average grain size. The electrical resistivity is more sensitive to grain size than is the Seebeck coefficient, and it is independent of whether or not the lattice thermal conductivity is decreased. This result can also be understood with a simple, 1D model, 33 Grain boundaries are thin, but as Fig. 2(b) shows, the resistivity of the grain boundaries is more than an order of magnitude larger than that of the grain. This occurs because band bending in the grain boundary exponentially decreases the carrier concentration (it linearly increases the Seebeck coefficient). The smaller the grains are, the more of these grain boundary regions exist and thus the higher the overall resistivity is. As a result, electrical resistivity is inversely proportional to the average grain size.
The thermal conductivity shown in Fig. 4(a) can also be explained in the same way: the more grain boundaries there are, the more interface scattering for phonons occurs, and thus the lower the thermal conductivity is. The power factor calculated using the aforementioned Seebeck coefficient, electrical resistivity, and thermal conductivity is shown in Fig. 4(b). For the particular set of input parameters chosen, PF exhibits a peak value at an average grain size of 3 μm. For smaller grain sizes, the increase in resistivity begins to dominate, and PF decreases. This echoes reports in the literature that show optimal grain sizes exist for various thermoelectric nanocomposites. 42 Finally, we examine the thermoelectric figure of merit of the simulated nanocomposites. The figure of merit was obtained in three different ways. The first was from the simulated measurements of the electrical conductivity, σ, the Seebeck coefficient, S, and the total thermal conductivity, κ, according to the definition where T is the temperature. This is the method most commonly used to report the measured zT of nanocomposite materials. The second way zT determined was from simulated Harman measurements, which we label as zT (Harman). The third way zT determined was from the maximum cooling temperature, ΔT max , of a simulated thermoelectric cooler, Because the simulated cooler is ideal with no contact resistances or thermal shunts, we expect values for zT(ΔT max ) that are similar to those obtained with the first two methods. (A small difference might be expected because of the temperature dependence of the thermoelectric parameters.) The key question in this paper is whether zT(σ, S, κ) predicted from Eq. (4), zT (Harman), and zT(ΔT max ) deduced from Eq. (5) are the same. Figure 5 provides the answer.
As shown in Fig. 5, for all average grain sizes, measurement methods (1) and (2) zT(σ, S, κ) and zT(Harman) are identical but significantly larger than for measurement method (3) zT(ΔT max ). In other words, ΔT max of the cooler built from this polycrystalline material would not be as high as would be expected. As discussed in the supplementary material, the conclusion for embedded nanoparticles is the same. The key point is that when potential barriers and lower thermal resistance occur, these effects should be expected.
Also shown as filled triangles in Fig. 5 are results for the case in which the electrical properties of the grain boundaries are identical to those of the grains and the crystalline reference. For this case, the only difference between the polycrystalline and crystalline samples is the lattice thermal conductivity. Figure 5 shows that for this case, zT is greater than for the crystalline reference and that all three measurement methods for obtaining zT give identical results.

IV. DISCUSSION
The difference between zT predicted from the individually measured electrical and thermal conductivities and Seebeck coefficient (or from a Harman measurement) and zT deduced from the cooler operating at ΔT max originates from the nature of the measurements. Seebeck coefficients are measured by creating a temperature gradient across a sample and then measuring the open-circuit voltage. To minimize the effects of temperature dependent material parameters, the temperature gradient applied across the sample is typically very small-a few Kelvin. 45,46 Harman measurements are typically done at low currents to minimize Joule heating. A thermoelectric cooler operates at high currents where Joule heating is significant.
One-dimensional simulations of the structure shown in Fig. 6 provide some insights into these results. Three different cases were considered: (1) A uniform material with the properties of the grains in the 2D simulations, (2) a 1D nanocomposite with reduced thermal conductivity in the grain boundaries, but with electrical properties that are idential to those of the grain, and (3) a 1D nanocomposite with a reduced thermal conductivity in the grain boundaries and with the lighter doping that produces potential barriers for carriers. Case 1 is the bulk reference. Case 2 is a nanocomposite that benefits only from a reduction of the lattice thermal conductivity, and case 3 is a nanocomposite that also provides an enhanced power factor. For case 3, 1D simulations give zT(σ, S, κ) ¼ zT(Harman) ¼ 1:35 and zT(ΔT max ) ¼ 1:02, which, for an average grain size of 1.3 μm, agrees well with the 2D results shown in Fig. 5. Figure 7 shows the ΔT vs I comparison for the three cases shown in Fig. 6. As expected, case 1, the bulk crystal case, has the lowest performance. For case 2, when the grains only have the added benefits of increased thermal resistance while their electrical properties are identical to those of grains, the cooling performance is significantly enhanced. Finally, note that case 3, for which the grains have an enhanced S and a lower thermal conductivity, gives a lower performance than case 2. This occurs in spite of the fact that the conventionally measured zT for case 3 is zT(σ, S, κ) ¼ 1:35 while that for case 2 is zT(σ, S, κ) ¼ 1:05. Figure 7 also shows ΔT vs I that would be predicted from the conventionally measured zT (1.35) for the case 3 nanocomposite. These results demonstrate that conventional measurements of nanocomposites with electrically active grain boundaries (case 3) substantially overpredict the performance of thermoelectric coolers (compare the solid and dashed lines in Fig. 7). They also show that while potential barriers at grain boundaries can increase the measured zT, they decrease the performance of a thermoelectric cooler in comparison to a nanocomposite with no potential barriers at the grain boundaries (i.e., case 3 vs case 2).
One way to understand these effects is in terms of the compatibility factor discussed by Snyder et al. 47 The grains are optimally doped to have the highest power factor according to Fig. 2, but the grain boundaries are not optimally doped. Intuitively, one might expect this effect to be small because the grain boundaries occupy a small portion of the overall geometry. The difference shown in Fig. 7 between case 2 and case 3 is, however, not insignificant. The additional performance degradation can be understood as due to the incompatibility between the grain and grain boundary elements. The essential point-of-view of the compatibility factor is to analyze the 1D thermoelectric device as infinitesimally small segments cascaded in series. Next, we adopt this view and compare the thermoelectric device under Seebeck measurement conditions and at ΔT max of a thermoelectric cooler.
As shown in Fig. 8 for the Seebeck measurement, the thermoelectric device is subject to a very small temperature difference, and the device is open-circuited electrically. The temperature profile shows the expected linear profile. The details of grain boundaries are visible in the insets of Fig. 8. Since grain boundary regions have high thermal resistance, more temperature drops across them, and the profile shows a zig-zag pattern. As discussed earlier in Eq. (2), this causes the overall Seebeck measurement to favor the Seebeck value of the grain boundaries, which is higher than that of the grains.
Under cooling operation, the electrical current is high, and the temperature variations are steeper than in the Seebeck measurement. There is also a significant amount of Joule heating, which cannot be ignored. The temperature profile at ΔT max is shown in Fig. 9. Overall, the profile behaves as expected for a homogenous device. The details, however, are important. At the cold side, the profile looks like the one seen for the Seebeck measurement, albeit with a higher temperature variation. Both the grain and grain  boundary regions operate as coolers, and the cooling adds up. At the hot side, the profile shows a sawtooth pattern. More importantly, the grain regions no longer act as coolers but as heaters instead. This occurs because the high resistance of the grain boundaries makes them significant Joule heating sources and the low thermal conductance makes it hard for the heat to flow away. In other words, the overall cooling device consists of many small segments of coolers in series with several heaters. This, of course, lowers the cooling performance.
To understand quantitatively why the temperature profile shown in Fig. 8 enhances the sample Seebeck coefficient and why the temperature profile in Fig. 9 does not, consider a one-dimensional example for which the hole current density is From the open-circuit voltage for a sample of length, L, we can deduce the Seebeck coefficient for the polycrystalline sample from Equation (7) shows that the Seebeck coefficient of the sample is strongly influenced by the Seebeck coefficient of regions where the temperature gradient large. When the Seebeck coefficient is measured under open-circuit conditions (Fig. 8), the temperature gradient in the grain boundaries is larger than in the grains, so the overall S increases. On the other hand, at ΔT max for the TEC, Joule heating produces a temperature gradient that changes sign in the grain boundary. According to Eq. (7), this change in sign  of the temperature gradient essentially eliminates the higher Seebeck coefficient of the grain boundary.

V. CONCLUSIONS
In this work, we examined the electronic performance of nanostructured thermoelectric materials using first-principles informed two-dimensional numerical simulation. It is now wellestablished that phonon scattering from grain boundaries or imbedded nanoparticles lowers the lattice thermal conductivity, which increases the material figure of merit, zT. Our focus was on the electronic performance of these materials-specifically on materials in which the depletion of carriers near grain boundaries (or around embedded nanoparticles) increases the Seebeck coefficient. Simulations of these materials with electrically active grain boundaries confirm that when the Seebeck coefficient of the grain boundaries is larger than that of the grains, the overall Seebeck coefficient of the sample only increases significantly when the lattice thermal conductivity is low where the Seebeck coefficient is high. More importantly, we showed that when zT of the nanocomposite is determined from standard measurements of σ, κ, and S, it provides an overly optimistic estimate of the performance of a thermoelectric cooler. When we compare a nanocomposite with reduced lattice thermal conductivity due to grain boundary scattering to another nanocomposite that also has potential barriers at the grain boundaries, we found that the one with potential barriers had a higher measured zT, but it produced a lower ΔT max in a TE cooler.
The parameter space for thermoelectric nanocomposites is large. We have explored only the case for which the grains are optimally doped, the grain boundaries deplete the grains and bend the bands by about one-half of the bandgap, and the thermal resistance of the grain boundary was reduced. (A similar case for embedded nanoparticles is discussed in the supplementary material.) The conclusions of this study are, however, expected to apply more broadlypotential barriers can increase the measured zT of a nanocomposite, but they can reduce the performance of a thermoelectric cooler made from the same nanocomposite. This occurs because electrically active grain boundaries of the type examined here (or of similarly electrically active embedded nanoparticles) are particularly sensitive to Joule heating, which occurs during device operation but not during the measurement of zT. Although we only examined thermoelectric coolers, the same conclusion may also apply to electrical power generation from thermoelectric devices.

SUPPLEMENTARY MATERIAL
See the supplementary material for extraction of thermoelectric material parameters from DFT calculation, the procedure for a random grain geometry generation and coupling with Sentaurus Device simulator, the results with Ag nanoparticles, and comparison between 2D and 1D simulations.